25#include "communication.hpp"
27#include "system/System.hpp"
31#include <boost/mpi/collectives/all_reduce.hpp>
32#include <boost/range/algorithm/find_if.hpp>
42 if (VolumeInitDone && !BoundariesFound) {
46 calc_volume_force(cs);
51 auto const &bonded_ias = *
get_system().bonded_ias;
54 if (not BoundariesFound) {
55 BoundariesFound = std::ranges::any_of(bonded_ias, [](
auto const &kv) {
56 return (boost::get<IBMVolCons>(&(*kv.second)) !=
nullptr);
60 if (!VolumeInitDone && BoundariesFound) {
66 for (
auto &kv : bonded_ias) {
67 if (
auto *v = boost::get<IBMVolCons>(kv.second.get())) {
71 BoundariesFound =
true;
72 if (v->volRef == 0.) {
73 v->volRef = VolumesCurrent[v->softID];
78 VolumeInitDone =
true;
85 auto const it = boost::find_if(p1.
bonds(), [&](
auto const &bond) ->
bool {
86 return boost::get<IBMVolCons>(bonded_ias.at(bond.bond_id()).get());
89 return (it != p1.
bonds().end())
90 ? boost::get<IBMVolCons>(bonded_ias.
at(it->bond_id()).get())
103 auto const &bonded_ias = *
get_system().bonded_ias;
106 std::vector<double> tempVol(VolumesCurrent.size());
109 cs.
bond_loop([&tempVol, &box_geo, &bonded_ias](
110 Particle &p1,
int bond_id, std::span<Particle *> partners) {
113 if (vol_cons_params &&
114 boost::get<IBMTriel>(bonded_ias.at(bond_id).get()) !=
nullptr) {
123 auto const x1 = box_geo.unfolded_position(p1.
pos(), p1.
image_box());
124 auto const x2 = x1 + box_geo.get_mi_vector(p2.
pos(), x1);
125 auto const x3 = x1 + box_geo.get_mi_vector(p3.
pos(), x1);
137 const double v321 = x3[0] * x2[1] * x1[2];
138 const double v231 = x2[0] * x3[1] * x1[2];
139 const double v312 = x3[0] * x1[1] * x2[2];
140 const double v132 = x1[0] * x3[1] * x2[2];
141 const double v213 = x2[0] * x1[1] * x3[2];
142 const double v123 = x1[0] * x2[1] * x3[2];
144 tempVol[vol_cons_params->softID] +=
145 1.0 / 6.0 * (-v321 + v231 + v312 - v132 - v213 + v123);
151 boost::mpi::all_reduce(
comm_cart, tempVol.data(),
152 static_cast<int>(tempVol.size()),
153 VolumesCurrent.data(), std::plus<double>());
157void ImmersedBoundaries::calc_volume_force(
CellStructure &cs) {
158 if (!BoundariesFound)
162 auto const &bonded_ias = *
get_system().bonded_ias;
165 std::span<Particle *> partners) {
166 if (boost::get<IBMTriel>(bonded_ias.at(bond_id).get()) !=
nullptr) {
171 auto const vol_cons_params = vol_cons_parameters(bonded_ias, p1);
172 if (not vol_cons_params)
175 auto const current_volume =
176 VolumesCurrent[static_cast<unsigned int>(vol_cons_params->softID)];
180 Particle &p2 = *partners[0];
181 Particle &p3 = *partners[1];
186 auto const x1 = box_geo.unfolded_position(p1.pos(), p1.image_box());
190 auto const a12 = box_geo.get_mi_vector(p2.pos(), x1);
191 auto const a13 = box_geo.get_mi_vector(p3.pos(), x1);
195 auto const n = vector_product(a12, a13);
196 const double ln = n.norm();
197 const double A = 0.5 * ln;
198 const double fact = vol_cons_params->kappaV *
199 (current_volume - vol_cons_params->volRef) /
202 auto const nHat = n / ln;
203 auto const force = -fact * A * nHat;
214 auto const new_size = bond.
softID + 1u;
215 if (new_size > VolumesCurrent.size()) {
216 VolumesCurrent.resize(new_size);
static IBMVolCons const * vol_cons_parameters(BondedInteractionsMap const &bonded_ias, Particle const &p1)
Data structures for bonded interactions.
container 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.
void register_softID(IBMVolCons &bond)
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.
void set_volumes_view(std::vector< double > const &volumes)
unsigned int softID
ID of the large soft particle to which this node belongs.
Struct holding all information for one particle.
auto const & image_box() const
auto const & bonds() const