Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/coupling/point_association_inline_impl.hh
2864081480a0lammps_2020

Last Change
Raw File

point_association_inline_impl.hh
View Options

/**
* @file point_association_inline_impl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Mon Jun 16 17:13:26 2014
*
* @brief computes the association between point degrees of freedom of
* different domains
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "filter_geometry.hh"
#include "point_association.hh"
#include "stimulation_zero.hh"
/* -------------------------------------------------------------------------- */
#include <string>
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA, typename ContB>
inline void PointAssociation::mixField(ContA &pointsAlist, ContB &pointsBlist) {
if (this->in_group_A() && !this->in_group_B())
LM_FATAL("not a parallel routine yet");
if (pointsAlist.size() != pointsBlist.size())
LM_FATAL("inconsistency here");
constexpr UInt Dim = ContA::Dim;
this->field_buffer.resize(pointsAlist.size(), Dim);
// make the average in the field_buffer thing
UInt index = 0;
UInt index_weight = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] = weightsA[index_weight] * _field[dir];
}
++index_weight;
}
index = 0;
index_weight = 0;
for (auto &&point : pointsBlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] += weightsB[index_weight] * _field[dir];
}
++index_weight;
}
// apply the result to the fields
index = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
_field[dir] = this->field_buffer[index];
}
}
index = 0;
for (auto &&point : pointsBlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
_field[dir] = this->field_buffer[index];
}
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA>
inline void PointAssociation::prepareSynchronisation(ContA &pointsAlist) {
if (ftype == _displacement) {
this->prepareSynchronisationDisplacement(pointsAlist);
return;
}
constexpr UInt Dim = ContA::Dim;
if (this->in_group_A()) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
UInt index = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] = _field[dir];
}
}
} else if (this->in_group_B()) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
}
}
/* -------------------------------------------------------------------------- */
// specialize this for ComputeTrueDisplacement
template <typename ContA>
inline void
PointAssociation::prepareSynchronisationDisplacement(ContA &pointsAlist) {
constexpr UInt Dim = ContA::Dim;
if (this->in_group_A()) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
DUMP(this->getID() << ": bounding box of " << pointsAlist.getID() << " is "
<< pointsAlist->getBoundingBox(),
DBG_INFO);
this->true_disp.build<dispatch>(this->pointsAlist);
UInt index = 0;
for (auto &&disp : this->true_disp.getOutputAsArray()) {
this->field_buffer[index] = disp;
DUMP(this->getID() << ": prepareSyncDisp.template field<" << _displacement
<< ">(" << index % Dim << ") = " << disp,
DBG_DETAIL);
++index;
}
} else if (this->in_group_B()) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype> inline void PointAssociation::synchronise() {
std::stringstream sstr;
sstr << "Sync " << ftype << " for object " << this->getID();
this->field_buffer.resize(this->local_points, spatial_dimension);
this->distributeVectorA2B(sstr.str(), this->field_buffer);
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA, typename ContB>
inline void PointAssociation::applyBuffer(ContA &pointsAlist,
ContB &pointsBlist) {
if (ftype == _position) {
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0)
this->template applyBufferPosition<false, false, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 0)
this->template applyBufferPosition<true, false, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 0)
this->template applyBufferPosition<true, true, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 1)
this->template applyBufferPosition<true, true, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 1)
this->template applyBufferPosition<true, false, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 0)
this->template applyBufferPosition<false, true, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 1)
this->template applyBufferPosition<false, true, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 1)
this->template applyBufferPosition<false, false, true>(pointsAlist,
pointsBlist);
return;
}
// constexpr UInt Dim = ContainerA::Dim;
for (UInt point_idx = 0; point_idx < pointsBlist->size(); ++point_idx) {
auto point = pointsBlist->get(point_idx);
DUMP(this->getID() << ": apply buffer for point " << point, DBG_DETAIL);
auto &&field_vals = this->field_buffer(point_idx);
auto &&_field = point.template field<ftype>();
_field = field_vals;
// for (UInt dir = 0; dir < Dim; ++dir) {
// _field[dir] = field_vals[dir];
// DUMP(this->getID() << ": point.template field<" << ftype << ">(" << dir
// << ") = " << field_vals[dir],
// DBG_DETAIL);
// }
}
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, UInt Dim, bool pbc_x,
bool pbc_y, bool pbc_z>
inline void fix_pbc(Vec1 &pos, Vec2 &xmin, Vec3 &d_size) {
for (UInt i = 0; i < Dim; ++i) {
if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) {
int offset = floor((pos[i] - xmin[i]) / d_size[i]);
pos[i] -= offset * d_size[i];
}
}
}
/* -------------------------------------------------------------------------- */
template <bool pbc_x, bool pbc_y, bool pbc_z, typename ContA, typename ContB>
inline void PointAssociation::applyBufferPosition(ContA &pointsAlist,
ContB &pointsBlist) {
DUMP(this->getID() << ": bounding box of " << this->pointsBlist.getID()
<< " is " << pointsBlist->getBoundingBox(),
DBG_INFO);
constexpr UInt Dim = ContA::Dim;
Cube &bbox = pointsBlist->getBoundingBox();
auto &&d_size = bbox.getSize();
auto &&xmin = bbox.getXmin<Dim>();
for (UInt point_idx = 0; point_idx < pointsBlist->size(); ++point_idx) {
auto point = pointsBlist->get(point_idx);
Vector<Dim> field_vals = this->field_buffer.row(point_idx);
DUMP(this->getID() << ": apply buffer position for point " << point,
DBG_DETAIL);
fix_pbc<Vector<Dim>, Vector<Dim>, Vector<3, Real>, Dim, pbc_x, pbc_y,
pbc_z>(field_vals, xmin, d_size);
auto &&_field = point.template field<_position>();
for (UInt dir = 0; dir < Dim; ++dir) {
_field[dir] = field_vals[dir];
DUMP(this->getID() << ": point.template field<" << _position << ">("
<< dir << ") = " << field_vals[dir],
DBG_DETAIL);
}
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype> inline void PointAssociation::enforceField() {
if (this->in_group_A()) {
STARTTIMER("Point Assoc. Synch.");
LM_TOIMPLEMENT;
// this->prepareSynchronisation<ftype>();
STOPTIMER("Point Assoc. Synch.");
}
this->synchronise<ftype>();
if (this->in_group_B()) {
STARTTIMER("Apply Assoc. Synch.");
LM_TOIMPLEMENT;
// this->applyBuffer<ftype>();
STOPTIMER("Apply Assoc. Synch.");
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype> inline void PointAssociation::zeroFields() {
StimulationZero stimulation_zero("zero_field:" + this->getID());
stimulation_zero.setParam("FIELD", ftype);
stimulation_zero.compute(pointsAlist);
stimulation_zero.compute(pointsBlist);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainB>
void PointAssociation::init(DomainA &domA, DomainB &domB) {
auto &all_group = Communicator::getCommunicator().getGroup("all");
auto &contA = domA.getContainer();
auto &contB = domB.getContainer();
UInt nb_zone_B = comm_group_B.size();
auto &pointsAlist =
this->pointsAlist.evalOutput<typename DomainA::ContainerSubset>();
auto &pointsBlist =
this->pointsBlist.evalOutput<typename DomainB::ContainerSubset>();
contA.addSubSet(pointsAlist);
contB.addSubSet(pointsBlist);
constexpr UInt Dim = std::decay_t<decltype(contA)>::Dim;
LM_TOIMPLEMENT;
// if (contA.hasRefManager() && contB.hasRefManager())
// LM_FATAL("Not sure what to do if both sides are able to migrate
// particles");
// Build point lists for both domains
if (this->in_group_A()) {
this->buildPointList(domA.getContainer(), this->unmatchedA);
LM_TOIMPLEMENT;
// this->local_points = this->unmatchedA.size();
// this->buildPositions(this->unmatchedA.getOutput());
}
if (this->in_group_B())
this->buildPointList(domB.getContainer(), this->unmatchedB);
DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP);
all_group.synchronize();
DUMP(this->getID() << " : exchange coarse geometries", DBG_DETAIL);
if (this->in_group_A()) {
this->exchangeGeometries<dispatch>(this->unmatchedA.evalOutput());
}
if (this->in_group_B()) {
this->exchangeGeometries<dispatch>(this->unmatchedB.evalOutput());
}
all_group.synchronize();
DUMP(this->getID() << " : exchange points positions", DBG_DETAIL);
this->exchangePositions(Dim);
all_group.synchronize();
if (this->in_group_B()) {
// set the size and assign zero (for help in debugging)
LM_TOIMPLEMENT;
// this->field_buffer.assign(Dim * (pointsAlist.size()), 0);
this->fillGrid();
this->associate();
}
all_group.synchronize();
DUMP(this->getID() << " : ExchangeAssociation", DBG_DETAIL);
ContainerArray<UInt> managed;
this->exchangeAssociationInformation(managed, point2Point);
if (!this->in_group_A() || !this->in_group_B()) {
if (this->in_group_A()) {
LM_ASSERT(managed.size() == this->local_points * nb_zone_B,
"internal error " << managed.size()
<< " != " << this->local_points * nb_zone_B);
std::vector<std::vector<UInt>> unassociated;
this->createDuoVectorA("A", managed, point2Point, unassociated);
for (UInt i = 0; i < nb_zone_B; ++i) {
LM_ASSERT(unassociated[i].size() == 0,
"internal error : "
<< " for proc " << i << " " << unassociated[i].size()
<< " points were already associated");
}
}
if (this->in_group_B())
this->createDuoVectorB("B", point2Point);
}
if (this->in_group_A()) {
this->filterPointListAForUnmatched();
LM_ASSERT(this->unmatchedA.size() == this->pointsAlist.size(),
"internal error " << this->unmatchedA.size()
<< " != " << this->pointsAlist.size());
LM_TOIMPLEMENT;
UInt total_points_from_file;
// UInt total_points_from_file = this->unmatchedA->getNbPointsInFile();
UInt total_points = this->local_points;
this->comm_group_A.allReduce(&total_points, 1, "verif association", OP_SUM);
if (total_points != total_points_from_file) {
if (total_points > total_points_from_file) {
LM_FATAL("This should never happen, not even with a faulty mesh. If "
<< "you get this message, you're looking at a serious bug.");
}
// in this case, need to figure out which points are missing
LM_FATAL("problem in A association " << total_points
<< " != " << total_points_from_file);
}
}
if (this->in_group_B()) {
this->filterPointListBForUnmatched();
this->local_points = this->pointsBlist.size();
this->filterMapping(point2Point);
LM_ASSERT(this->unmatchedB.size() == this->pointsBlist.size(),
"internal error");
#ifndef LM_OPTIMIZED
LM_TOIMPLEMENT;
UInt total_points_from_file;
// UInt total_points_from_file = this->unmatchedB.getNbPointsInFile();
UInt total_points = this->local_points;
this->comm_group_B.allReduce(&total_points, 1, "verif association", OP_SUM);
LM_ASSERT(total_points == total_points_from_file,
"problem in B association " << total_points
<< " != " << total_points_from_file);
#endif
}
checkPositionCoherency();
if (this->in_group_A()) {
pointsAlist.attachObject(this->getDuoVector());
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
}
if (this->in_group_B()) {
pointsBlist.attachObject(this->getDuoVector());
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
}
// DUMP(this->getID() << ": bounding box of " <<
// this->pointsAlist.getID() << " is "
// << this->pointsAlist.getBoundingBox(),DBG_MESSAGE);
// DUMP(this->getID() << ": bounding box of " <<
// this->pointsBlist.getID() << " is "
// << this->pointsBlist.getBoundingBox(),DBG_MESSAGE);
// generate the weights
if (weight_path != "") {
this->computeWeights(domA.getContainer(), domB.getContainer());
}
}
/* -------------------------------------------------------------------------- */
template <typename ContA, typename ContB>
void PointAssociation::computeWeights(ContA &pointsAlist, ContB &pointsBlist) {
ComputeExtract c_positions("ComputeExtract:" + this->getID());
c_positions.setParam("FIELD", _position0);
c_positions.compute(pointsAlist);
this->weights.setParam("FILENAME", this->weight_path);
this->weights.init();
this->weights.compute(c_positions);
if (this->in_group_A() && !this->in_group_B())
LM_FATAL("not a parallel routine yet");
if (this->pointsAlist.size() != this->pointsBlist.size())
LM_FATAL("inconsistency here");
std::vector<Real> denominator;
weightsA.assign(pointsAlist.size(), 0.);
weightsB.assign(pointsBlist.size(), 0.);
ComputeExtract mass("mass:" + this->getID());
mass.setParam("FIELD", _mass);
mass.compute(pointsAlist);
weightsA = weights.evalArrayOutput() * mass.evalArrayOutput().array();
mass.compute(pointsBlist);
weightsB = (1. - weights.evalArrayOutput()) * mass.evalArrayOutput().array();
auto denom = weightsA + weightsB;
weightsA.array() /= denom;
weightsB.array() /= denom;
}
/* -------------------------------------------------------------------------- */
template <typename container, typename PointList>
void PointAssociation::buildPointList(container &cont, PointList &pointList) {
FilterGeometry filter_geom("FilterGeom:" + this->getID());
filter_geom.setParam("GEOMETRY", this->geomID);
filter_geom.setParam("DONOTRESETBOUNDINGBOX", true);
filter_geom.setParam("DO_NOT_REGISTER", true);
filter_geom.setParam("BASED_ON_NODES", true);
filter_geom.setParam("DOF_TYPE", dt_master);
filter_geom.init();
filter_geom.compute(cont);
DUMP("constructed filtergeom with "
<< filter_geom.evalOutput<ContainerInterface>().size() << " points ",
DBG_INFO_STARTUP);
pointList.setParam("PATH", this->path);
pointList.setParam("TOL", this->tolerance);
pointList.init();
pointList.compute(filter_geom);
DUMP("constructed filter_point with "
<< pointList.template evalOutput<ContainerInterface>().size()
<< " points ",
DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help