26#include "communication.hpp"
28#include "system/System.hpp"
33#include <boost/mpi/collectives/all_reduce.hpp>
34#include <boost/serialization/utility.hpp>
49 cs.
bond_loop([&area, &volume, &box_geo, &bonded_ias, molType](
50 Particle &p1,
int bond_id, std::span<Particle *> partners) {
51 if (p1.
mol_id() != molType)
54 if (boost::get<OifGlobalForcesBond>(bonded_ias.
at(bond_id).get())) {
56 auto const p22 = p11 + box_geo.
get_mi_vector(partners[0]->pos(), p11);
57 auto const p33 = p11 + box_geo.
get_mi_vector(partners[1]->pos(), p11);
63 auto const VOL_dn = VOL_norm.norm();
64 auto const VOL_hz = (1.0 / 3.0) * (p11[2] + p22[2] + p33[2]);
65 volume -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz;
79 cs.
bond_loop([&box_geo, &bonded_ias, area, volume, molType](
80 Particle &p1,
int bond_id, std::span<Particle *> partners) {
81 if (p1.
mol_id() != molType)
84 auto const *bond_ptr = bonded_ias.
at(bond_id).get();
85 if (
auto const *bond = boost::get<OifGlobalForcesBond>(bond_ptr)) {
87 auto const p22 = p11 + box_geo.
get_mi_vector(partners[0]->pos(), p11);
88 auto const p33 = p11 + box_geo.
get_mi_vector(partners[1]->pos(), p11);
94 auto const VOL_vv = (volume - bond->V0) / bond->V0;
96 auto const VOL_force = (1. / 3.) * bond->kv * VOL_vv * VOL_A * VOL_norm;
98 auto const h = (1. / 3.) * (p11 + p22 + p33);
100 auto const deltaA = (area - bond->A0_g) / bond->A0_g;
102 auto const m1 = h - p11;
103 auto const m2 = h - p22;
104 auto const m3 = h - p33;
106 auto const m1_length = m1.norm();
107 auto const m2_length = m2.norm();
108 auto const m3_length = m3.norm();
110 auto const fac = bond->ka_g * VOL_A * deltaA /
111 (m1_length * m1_length + m2_length * m2_length +
112 m3_length * m3_length);
114 p1.
force() += fac * m1 + VOL_force;
115 partners[0]->force() += fac * m2 + VOL_force;
116 partners[1]->force() += fac * m3 + VOL_force;
123void OifGlobal::run_force_kernel() {
125 auto &box_geo = *system.box_geo;
126 auto &bonded_ias = *system.bonded_ias;
127 auto &cell_structure = *system.cell_structure;
131 auto const local =
calc_oif_mesh(i, box_geo, cell_structure, bonded_ias);
132 auto const global = boost::mpi::all_reduce(
comm_cart, local, std::plus());
133 auto const area = std::abs(global[0]);
134 auto const volume = std::abs(global[1]);
135 if (area < 1e-100 and volume < 1e-100) {
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
container for bonded interactions.
mapped_type at(key_type const &key) const
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
boost::mpi::communicator comm_cart
The communicator.
double area_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the area of triangle between vectors P1,P2,P3, by computing the cross product P1P2 x P1P3 an...
Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the normal vector of a triangle.
static void add_oif_global_forces(double area, double volume, int molType, BoxGeometry const &box_geo, CellStructure &cs, BondedInteractionsMap const &bonded_ias)
Distribute the OIF global forces to all particles in the mesh.
static auto calc_oif_mesh(int molType, BoxGeometry const &box_geo, CellStructure &cs, BondedInteractionsMap const &bonded_ias)
Calculate the mesh volume and area.
Routines to calculate the OIF global forces for a particle triple (triangle from mesh).
Describes a cell structure / cell system.
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
Struct holding all information for one particle.
auto const & image_box() const
auto const & mol_id() const
auto const & force() const