Skip to content

Commit

Permalink
bond-refactoring: fix compiler warnings
Browse files Browse the repository at this point in the history
  • Loading branch information
@biermanncarl
biermanncarl committed Mar 17, 2021
1 parent 1fb0dc8 commit 4c80030
Show file tree
Hide file tree
Showing 7 changed files with 130 additions and 135 deletions.
4 changes: 2 additions & 2 deletions src/core/bonded_interactions/bonded_coulomb.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,9 @@ struct Bonded_coulomb_bond_parameters {
Bonded_coulomb_bond_parameters() = default;
Bonded_coulomb_bond_parameters(double prefactor);

boost::optional<Utils::Vector3d> pair_force(double const q1q2,
boost::optional<Utils::Vector3d> pair_force(double q1q2,
Utils::Vector3d const &dx) const;
boost::optional<double> pair_energy(double const q1q2,
boost::optional<double> pair_energy(double q1q2,
Utils::Vector3d const &dx) const;

private:
Expand Down
9 changes: 2 additions & 7 deletions src/core/bonded_interactions/bonded_interaction_data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,13 +46,8 @@ double maximal_cutoff_bonded() {
/* Check if there are dihedrals */
auto const any_dihedrals = std::any_of(
bonded_ia_params.begin(), bonded_ia_params.end(), [](auto const &bond) {
if (boost::get<Dihedral_bond_parameters>(&bond)) {
return true;
} else if (boost::get<Tabulated_dihedral_bond_parameters>(&bond)) {
return true;
} else {
return false;
}
return (boost::get<Dihedral_bond_parameters>(&bond) ||
boost::get<Tabulated_dihedral_bond_parameters>(&bond));
});

/* dihedrals: the central particle is indirectly connected to the fourth
Expand Down
64 changes: 32 additions & 32 deletions src/core/energy_inline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,68 +213,68 @@ calc_bonded_energy(Bonded_ia_parameters const &iaparams, Particle const &p1,
auto const dx = get_mi_vector(p1.r.p, p2->r.p, box_geo);
if (auto const *iap = boost::get<Fene_bond_parameters>(&iaparams)) {
return iap->pair_energy(dx);
} else if (auto const *iap =
boost::get<Harmonic_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Harmonic_bond_parameters>(&iaparams)) {
return iap->pair_energy(dx);
} else if (auto const *iap =
boost::get<Quartic_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Quartic_bond_parameters>(&iaparams)) {
return iap->pair_energy(dx);
}
#ifdef ELECTROSTATICS
else if (auto const *iap =
boost::get<Bonded_coulomb_bond_parameters>(&iaparams)) {
if (auto const *iap =
boost::get<Bonded_coulomb_bond_parameters>(&iaparams)) {
return iap->pair_energy(p1.p.q * p2->p.q, dx);

} else if (auto const *iap =
boost::get<Bonded_coulomb_sr_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Bonded_coulomb_sr_bond_parameters>(&iaparams)) {
return iap->pair_energy(p1, *p2, dx);
}
#endif
#ifdef BOND_CONSTRAINT
else if (auto const *iap = boost::get<Rigid_bond_parameters>(&iaparams)) {
if (auto const *iap = boost::get<Rigid_bond_parameters>(&iaparams)) {
return boost::optional<double>(0);
}
#endif
#ifdef TABULATED
else if (auto const *iap =
boost::get<Tabulated_distance_bond_parameters>(&iaparams)) {
if (auto const *iap =
boost::get<Tabulated_distance_bond_parameters>(&iaparams)) {
return iap->pair_energy(dx);
}
#endif
else if (auto const *iap = boost::get<VirtualBond_Parameters>(&iaparams)) {
if (auto const *iap = boost::get<VirtualBond_Parameters>(&iaparams)) {
return boost::optional<double>(0);
} else {
throw BondUnknownTypeError();
}
throw BondUnknownTypeError();
} // 1 partner
else if (n_partners == 2) {
if (n_partners == 2) {
if (auto const *iap =
boost::get<Angle_harmonic_bond_parameters>(&iaparams)) {
return iap->angle_energy(p1.r.p, p2->r.p, p3->r.p);
} else if (auto const *iap =
boost::get<Angle_cosine_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Angle_cosine_bond_parameters>(&iaparams)) {
return iap->angle_energy(p1.r.p, p2->r.p, p3->r.p);
} else if (auto const *iap =
boost::get<Angle_cossquare_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Angle_cossquare_bond_parameters>(&iaparams)) {
return iap->angle_energy(p1.r.p, p2->r.p, p3->r.p);
} else if (auto const *iap =
boost::get<Tabulated_angle_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Tabulated_angle_bond_parameters>(&iaparams)) {
return iap->angle_energy(p1.r.p, p2->r.p, p3->r.p);
} else {
throw BondUnknownTypeError();
}
} // 2 partner
else if (n_partners == 3) {
throw BondUnknownTypeError();
} // 2 partners
if (n_partners == 3) {
if (auto const *iap = boost::get<Dihedral_bond_parameters>(&iaparams)) {
return iap->dihedral_energy(p2->r.p, p1.r.p, p3->r.p, p4->r.p);
} else if (auto const *iap =
boost::get<Tabulated_dihedral_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Tabulated_dihedral_bond_parameters>(&iaparams)) {
return iap->dihedral_energy(p2->r.p, p1.r.p, p3->r.p, p4->r.p);
} else {
throw BondUnknownTypeError();
}
return 0.;
} else if (n_partners == 0) {
throw BondUnknownTypeError();
} // 3 partners
if (n_partners == 0) {
return 0.;
}

Expand Down
78 changes: 38 additions & 40 deletions src/core/forces_inline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,33 +308,33 @@ calc_bond_pair_force(Particle const &p1, Particle const &p2,
Utils::Vector3d const &dx) {
if (auto const *iap = boost::get<Fene_bond_parameters>(&iaparams)) {
return iap->pair_force(dx);
} else if (auto const *iap =
boost::get<Harmonic_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Harmonic_bond_parameters>(&iaparams)) {
return iap->pair_force(dx);
} else if (auto const *iap = boost::get<Quartic_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Quartic_bond_parameters>(&iaparams)) {
return iap->pair_force(dx);
}
#ifdef ELECTROSTATICS
else if (auto const *iap =
boost::get<Bonded_coulomb_bond_parameters>(&iaparams)) {
if (auto const *iap = boost::get<Bonded_coulomb_bond_parameters>(&iaparams)) {
return iap->pair_force(p1.p.q * p2.p.q, dx);
} else if (auto const *iap =
boost::get<Bonded_coulomb_sr_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Bonded_coulomb_sr_bond_parameters>(&iaparams)) {
return iap->pair_force(dx);
}
#endif
#ifdef TABULATED
else if (auto const *iap =
boost::get<Tabulated_distance_bond_parameters>(&iaparams)) {
if (auto const *iap =
boost::get<Tabulated_distance_bond_parameters>(&iaparams)) {
return iap->pair_force(dx);
}
#endif
else if (boost::get<VirtualBond_Parameters>(&iaparams) ||
boost::get<Rigid_bond_parameters>(&iaparams)) {
if (boost::get<VirtualBond_Parameters>(&iaparams) ||
boost::get<Rigid_bond_parameters>(&iaparams)) {
return Utils::Vector3d{};
} else {
throw BondUnknownTypeError();
}
throw BondUnknownTypeError();
}

inline bool add_bonded_two_body_force(Bonded_ia_parameters const &iaparams,
Expand Down Expand Up @@ -372,24 +372,24 @@ calc_bonded_three_body_force(Bonded_ia_parameters const &iaparams,
Particle const &p3) {
if (auto const *iap = boost::get<Angle_harmonic_bond_parameters>(&iaparams)) {
return iap->angle_force(p1.r.p, p2.r.p, p3.r.p);
} else if (auto const *iap =
boost::get<Angle_cosine_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Angle_cosine_bond_parameters>(&iaparams)) {
return iap->angle_force(p1.r.p, p2.r.p, p3.r.p);
} else if (auto const *iap =
boost::get<Angle_cossquare_bond_parameters>(&iaparams)) {
}
if (auto const *iap =
boost::get<Angle_cossquare_bond_parameters>(&iaparams)) {
return iap->angle_force(p1.r.p, p2.r.p, p3.r.p);
}
#ifdef TABULATED
else if (auto const *iap =
boost::get<Tabulated_angle_bond_parameters>(&iaparams)) {
if (auto const *iap =
boost::get<Tabulated_angle_bond_parameters>(&iaparams)) {
return iap->angle_force(p1.r.p, p2.r.p, p3.r.p);
}
#endif
else if (auto const *iap = boost::get<IBM_Triel_Parameters>(&iaparams)) {
if (auto const *iap = boost::get<IBM_Triel_Parameters>(&iaparams)) {
return iap->calc_forces(p1, p2, p3);
} else {
throw BondUnknownTypeError();
}
throw BondUnknownTypeError();
}

inline bool add_bonded_three_body_force(Bonded_ia_parameters const &iaparams,
Expand All @@ -398,18 +398,17 @@ inline bool add_bonded_three_body_force(Bonded_ia_parameters const &iaparams,
if (auto const *iap =
boost::get<Oif_global_forces_bond_parameters>(&iaparams)) {
return false;
} else {
auto const result = calc_bonded_three_body_force(iaparams, p1, p2, p3);
if (result) {
using std::get;
auto const &forces = result.get();
}
auto const result = calc_bonded_three_body_force(iaparams, p1, p2, p3);
if (result) {
using std::get;
auto const &forces = result.get();

p1.f.f += get<0>(forces);
p2.f.f += get<1>(forces);
p3.f.f += get<2>(forces);
p1.f.f += get<0>(forces);
p2.f.f += get<1>(forces);
p3.f.f += get<2>(forces);

return false;
}
return false;
}
return true;
}
Expand All @@ -422,21 +421,20 @@ calc_bonded_four_body_force(Bonded_ia_parameters const &iaparams,
if (auto const *iap =
boost::get<Oif_local_forces_bond_parameters>(&iaparams)) {
return iap->calc_forces(p1, p2, p3, p4);
} else if (auto const *iap = boost::get<IBM_Tribend_Parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<IBM_Tribend_Parameters>(&iaparams)) {
return iap->calc_forces(p1, p2, p3, p4);
} else if (auto const *iap =
boost::get<Dihedral_bond_parameters>(&iaparams)) {
}
if (auto const *iap = boost::get<Dihedral_bond_parameters>(&iaparams)) {
return iap->dihedral_force(p2.r.p, p1.r.p, p3.r.p, p4.r.p);
}
#ifdef TABULATED
else if (auto const *iap =
boost::get<Tabulated_dihedral_bond_parameters>(&iaparams)) {
if (auto const *iap =
boost::get<Tabulated_dihedral_bond_parameters>(&iaparams)) {
return iap->dihedral_force(p2.r.p, p1.r.p, p3.r.p, p4.r.p);
}
#endif
else {
throw BondUnknownTypeError();
}
throw BondUnknownTypeError();
}

inline bool add_bonded_four_body_force(Bonded_ia_parameters const &iaparams,
Expand Down
100 changes: 51 additions & 49 deletions src/core/immersed_boundary/ImmersedBoundaries.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,8 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) {
if (not BoundariesFound) {
BoundariesFound = boost::algorithm::any_of(
bonded_ia_params, [](auto const &bonded_ia_param) {
return (boost::get<IBM_VolCons_Parameters>(&bonded_ia_param) != NULL);
return (boost::get<IBM_VolCons_Parameters>(&bonded_ia_param) !=
nullptr);
});
}

Expand Down Expand Up @@ -83,7 +84,7 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) {
static const IBM_VolCons_Parameters *vol_cons_parameters(Particle const &p1) {
auto it = boost::find_if(p1.bonds(), [](auto const &bond) {
return (boost::get<IBM_VolCons_Parameters>(
&bonded_ia_params[bond.bond_id()]) != NULL);
&bonded_ia_params[bond.bond_id()]) != nullptr);
});

return (it != p1.bonds().end()) ? boost::get<IBM_VolCons_Parameters>(
Expand All @@ -106,7 +107,7 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) {
auto vol_cons_params = vol_cons_parameters(p1);

if (vol_cons_params && (boost::get<IBM_Triel_Parameters>(
&bonded_ia_params[bond_id]) != NULL)) {
&bonded_ia_params[bond_id]) != nullptr)) {
// Our particle is the leading particle of a triel
// Get second and third particle of the triangle
Particle &p2 = *partners[0];
Expand Down Expand Up @@ -152,51 +153,52 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) {

/** Calculate and add the volume force to each node */
void ImmersedBoundaries::calc_volume_force(CellStructure &cs) {
cs.bond_loop([this](Particle &p1, int bond_id,
Utils::Span<Particle *> partners) {
if (boost::get<IBM_Triel_Parameters>(&bonded_ia_params[bond_id]) != NULL) {
// Check if particle has a IBM_Triel bonded interaction and a
// IBM_VolCons bonded interaction. Basically this loops over all
// triangles, not all particles. First round to check for volume
// conservation.
const IBM_VolCons_Parameters *ibmVolConsParameters =
vol_cons_parameters(p1);
if (not ibmVolConsParameters)
return false;
cs.bond_loop(
[this](Particle &p1, int bond_id, Utils::Span<Particle *> partners) {
if (boost::get<IBM_Triel_Parameters>(&bonded_ia_params[bond_id]) !=
nullptr) {
// Check if particle has a IBM_Triel bonded interaction and a
// IBM_VolCons bonded interaction. Basically this loops over all
// triangles, not all particles. First round to check for volume
// conservation.
const IBM_VolCons_Parameters *ibmVolConsParameters =
vol_cons_parameters(p1);
if (not ibmVolConsParameters)
return false;

auto current_volume = VolumesCurrent[ibmVolConsParameters->softID];

auto current_volume = VolumesCurrent[ibmVolConsParameters->softID];

// Our particle is the leading particle of a triel
// Get second and third particle of the triangle
Particle &p2 = *partners[0];
Particle &p3 = *partners[1];

// Unfold position of first node.
// This is to get a continuous trajectory with no jumps when box
// boundaries are crossed.
auto const x1 = unfolded_position(p1.r.p, p1.l.i, box_geo.length());

// Unfolding seems to work only for the first particle of a triel
// so get the others from relative vectors considering PBC
auto const a12 = get_mi_vector(p2.r.p, x1, box_geo);
auto const a13 = get_mi_vector(p3.r.p, x1, box_geo);

// Now we have the true and good coordinates
// This is eq. (9) in @cite dupin08a.
auto const n = vector_product(a12, a13);
const double ln = n.norm();
const double A = 0.5 * ln;
const double fact = ibmVolConsParameters->kappaV *
(current_volume - ibmVolConsParameters->volRef) /
current_volume;

auto const nHat = n / ln;
auto const force = -fact * A * nHat;

p1.f.f += force;
p2.f.f += force;
p3.f.f += force;
}
return false;
});
// Our particle is the leading particle of a triel
// Get second and third particle of the triangle
Particle &p2 = *partners[0];
Particle &p3 = *partners[1];

// Unfold position of first node.
// This is to get a continuous trajectory with no jumps when box
// boundaries are crossed.
auto const x1 = unfolded_position(p1.r.p, p1.l.i, box_geo.length());

// Unfolding seems to work only for the first particle of a triel
// so get the others from relative vectors considering PBC
auto const a12 = get_mi_vector(p2.r.p, x1, box_geo);
auto const a13 = get_mi_vector(p3.r.p, x1, box_geo);

// Now we have the true and good coordinates
// This is eq. (9) in @cite dupin08a.
auto const n = vector_product(a12, a13);
const double ln = n.norm();
const double A = 0.5 * ln;
const double fact = ibmVolConsParameters->kappaV *
(current_volume - ibmVolConsParameters->volRef) /
current_volume;

auto const nHat = n / ln;
auto const force = -fact * A * nHat;

p1.f.f += force;
p2.f.f += force;
p3.f.f += force;
}
return false;
});
}
Loading

0 comments on commit 4c80030

Please sign in to comment.