25#include "communication.hpp"
27#include "system/System.hpp"
35#include <boost/mpi/collectives/all_reduce.hpp>
36#include <boost/range/algorithm/find_if.hpp>
44 if (VolumeInitDone && !BoundariesFound) {
48 calc_volume_force(cs);
55 if (not BoundariesFound) {
56 BoundariesFound = std::any_of(
58 return (boost::get<IBMVolCons>(&(*kv.second)) != nullptr);
62 if (!VolumeInitDone && BoundariesFound) {
69 if (
auto *v = boost::get<IBMVolCons>(&(*kv.second))) {
73 BoundariesFound =
true;
74 if (v->volRef == 0.) {
75 v->volRef = VolumesCurrent[
static_cast<unsigned int>(v->softID)];
80 VolumeInitDone =
true;
85 auto const it = boost::find_if(p1.
bonds(), [](
auto const &bond) ->
bool {
86 return boost::get<IBMVolCons>(bonded_ia_params.at(bond.bond_id()).get());
89 return (it != p1.
bonds().end())
105 std::vector<double> tempVol(VolumesCurrent.size());
112 if (vol_cons_params &&
122 auto const x1 = box_geo.unfolded_position(p1.
pos(), p1.
image_box());
123 auto const x2 = x1 + box_geo.get_mi_vector(p2.
pos(), x1);
124 auto const x3 = x1 + box_geo.get_mi_vector(p3.
pos(), x1);
136 const double v321 = x3[0] * x2[1] * x1[2];
137 const double v231 = x2[0] * x3[1] * x1[2];
138 const double v312 = x3[0] * x1[1] * x2[2];
139 const double v132 = x1[0] * x3[1] * x2[2];
140 const double v213 = x2[0] * x1[1] * x3[2];
141 const double v123 = x1[0] * x2[1] * x3[2];
143 tempVol[
static_cast<unsigned int>(vol_cons_params->softID)] +=
144 1.0 / 6.0 * (-v321 + v231 + v312 - v132 - v213 + v123);
150 boost::mpi::all_reduce(
comm_cart, tempVol.data(),
151 static_cast<int>(tempVol.size()),
152 VolumesCurrent.data(), std::plus<double>());
156void ImmersedBoundaries::calc_volume_force(
CellStructure &cs) {
157 if (!BoundariesFound)
169 auto const vol_cons_params = vol_cons_parameters(p1);
170 if (not vol_cons_params)
173 auto const current_volume =
174 VolumesCurrent[static_cast<unsigned int>(vol_cons_params->softID)];
178 Particle &p2 = *partners[0];
179 Particle &p3 = *partners[1];
184 auto const x1 = box_geo.unfolded_position(p1.pos(), p1.image_box());
188 auto const a12 = box_geo.get_mi_vector(p2.pos(), x1);
189 auto const a13 = box_geo.get_mi_vector(p3.pos(), x1);
193 auto const n = vector_product(a12, a13);
194 const double ln = n.norm();
195 const double A = 0.5 * ln;
196 const double fact = vol_cons_params->kappaV *
197 (current_volume - vol_cons_params->volRef) /
200 auto const nHat = n / ln;
201 auto const force = -fact * A * nHat;
static const IBMVolCons * vol_cons_parameters(Particle const &p1)
Vector implementation and trait types for boost qvm interoperability.
BondedInteractionsMap bonded_ia_params
Field containing the parameters of the bonded ia types.
Data structures for bonded interactions.
mapped_type at(key_type const &key) const
void volume_conservation(CellStructure &cs)
Calculate volumes, volume force and add it to each virtual particle.
void init_volume_conservation(CellStructure &cs)
Initialize volume conservation.
std::shared_ptr< BoxGeometry > box_geo
A stripped-down version of std::span from C++17.
boost::mpi::communicator comm_cart
The communicator.
Describes a cell structure / cell system.
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
Parameters for IBM volume conservation bond.
Struct holding all information for one particle.
auto const & image_box() const
auto const & bonds() const