Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/coupling/point_association.cc
c73f13fec450features/coupler_compon…

Last Change
Raw File

point_association.cc
View Options

/**
* @file point_association.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Thu Jul 24 14:21:58 2014
*
* @brief Coupling tool to associate point-like DOFs
*
* @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/>.
*
*/
//#define TIMER
/* -------------------------------------------------------------------------- */
#include "point_association.hh"
#include "filter_geometry.hh"
#include "geometry_manager.hh"
#include "lib_bridging.hh"
#include "lib_continuum.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "reference_manager_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
PointAssociation::PointAssociation(const std::string &name)
: LMObject(name), DofAssociation(name), pointsAlist("pointsAlist:" + name),
pointsBlist("pointsBlist:" + name), unmatchedA("unmatchedA:" + name),
unmatchedB("unmatchedB:" + name), assoc_found(0), tolerance(1e-5),
path(""), weight_path(""), weights("weights:" + this->getID()),
weightsA("weightsA:" + this->getID()),
weightsB("weightsB:" + this->getID()),
true_disp("ComputeTrueDisplacement:" + this->getID()) {}
/* -------------------------------------------------------------------------- */
PointAssociation::~PointAssociation() {}
/* -------------------------------------------------------------------------- */
template <typename ContA, typename ContB>
void PointAssociation::init(ContA &contA, ContB &contB) {
if (this->in_group_A && contA.hasRefManager())
contA.getRefManager()->addSubSet(this->pointsAlist);
if (this->in_group_B && contB.hasRefManager())
contB.getRefManager()->addSubSet(this->pointsBlist);
/* registering computes for outer world */
LM_TOIMPLEMENT;
// FilterManager::getManager().addObject(&this->pointsAlist, false);
// FilterManager::getManager().addObject(&this->pointsBlist, false);
constexpr UInt Dim = ContA::Dim;
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(contA, this->unmatchedA);
LM_TOIMPLEMENT;
// this->local_points = this->unmatchedA.size();
// this->buildPositions(this->unmatchedA.getOutput());
}
if (this->in_group_B)
this->buildPointList(contB, this->unmatchedB);
DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP);
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : exchange coarse geometries", DBG_DETAIL);
this->exchangeGeometries(this->unmatchedA, this->unmatchedB);
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : exchange points positions", DBG_DETAIL);
this->exchangePositions(Dim);
MPI_Barrier(MPI_COMM_WORLD);
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();
}
MPI_Barrier(MPI_COMM_WORLD);
DUMP(this->getID() << " : ExchangeAssociation", DBG_DETAIL);
ContainerArray<UInt> managed;
this->exchangeAssociationInformation(managed);
if (!this->in_group_A || !this->in_group_B) {
if (this->in_group_A) {
LM_ASSERT(managed.size() == this->local_points * this->nb_zone_B,
"internal error " << managed.size() << " != "
<< this->local_points * this->nb_zone_B);
std::vector<std::vector<UInt>> unassociated;
this->createDuoVectorA("A", managed, unassociated);
for (UInt i = 0; i < this->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");
}
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->filterAssoc();
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 && contA.hasRefManager()) {
try {
contA.getRefManager()->attachObject(this->getDuoVector(),
this->pointsAlist);
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
} catch (LibMultiScaleException &e) {
}
}
if (this->in_group_B && contB.hasRefManager()) {
try {
contB.getRefManager()->attachObject(this->getDuoVector(),
this->pointsBlist);
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
} catch (LibMultiScaleException &e) {
}
}
// 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(contA, contB);
}
}
/* -------------------------------------------------------------------------- */
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.build<dispatch>(pointsAlist);
this->weights.setParam("FILENAME", this->weight_path);
this->weights.init();
this->weights.build<dispatch>(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.resize(this->pointsAlist.nbElem());
// weightsB.resize(this->pointsAlist.nbElem());
weightsA.assign(pointsAlist.size(), 0.);
weightsB.assign(pointsAlist.size(), 0.);
// compute nominator and denominator
UInt index = 0;
for (auto &&point : pointsAlist) {
Real mass = point.mass();
weightsA[index] = weights.getOutputAsArray()[index] * mass;
++index;
}
index = 0;
for (auto &&point : pointsBlist) {
Real mass = point.mass();
weightsB[index] = (1. - weights.getOutputAsArray()[index]) * mass;
++index;
}
for (UInt index = 0; index < weights.getOutputAsArray().size(); ++index) {
Real denom = weightsA[index] + weightsB[index];
weightsA[index] /= denom;
weightsB[index] /= denom;
}
LM_TOIMPLEMENT;
// FilterManager::getManager().addObject(&this->weights, false);
// FilterManager::getManager().addObject(&this->weightsA, false);
// FilterManager::getManager().addObject(&this->weightsB, false);
}
/* -------------------------------------------------------------------------- */
template <typename container, typename PointList>
void PointAssociation::buildPointList(container &cont, PointList &pointList) {
casted_component<FilterGeometry>::cast<typename container::ContainerSubset>
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_local_master);
filter_geom.init();
filter_geom.buildManual(cont);
DUMP("constructed filtergeom with " << filter_geom->size() << " points ",
DBG_INFO_STARTUP);
pointList.setParam("PATH", this->path);
pointList.setParam("TOL", this->tolerance);
pointList.init();
pointList.buildManual(filter_geom);
DUMP("constructed filter_point with " << pointList->size() << " points ",
DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
void PointAssociation::fillGrid() {
LM_TOIMPLEMENT;
// constexpr UInt Dim = ContainerA::Dim;
// create my grid
// this->grid.reset();
// Geometry &geom = *GeometryManager::getManager().getGeometry(this->geomID);
// Cube cube = geom.getBoundingBox();
// LM_TOIMPLEMENT;
// // this->grid = new SpatialGridPoint<UInt, Dim>(cube, this->grid_division);
// UInt index = 0;
// for (auto &&point : this->unmatchedB) {
// auto position = point.position0();
// // I'ts important to use addAtom here instead of addElement, because here
// // we're sensitive to atoms on the edge of grid_boxes
// this->grid->addAtom(index, position);
// ++index;
// }
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename V1, typename V2>
inline bool coincidence(V1 &pointA, V2 &pointB, const Real &tolerance) {
Real distance = 0;
for (UInt dir = 0; dir < Dim; ++dir) {
distance += fabs(pointA[dir] - pointB[dir]);
}
return distance < tolerance;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename V1, typename V2>
inline Real temp_coincidence(V1 &pointA, V2 &pointB) {
Real distance = 0;
for (UInt dir = 0; dir < Dim; ++dir) {
distance += fabs(pointA[dir] - pointB[dir]);
}
return distance;
}
/* -------------------------------------------------------------------------- */
void PointAssociation::associate() {
LM_TOIMPLEMENT;
// this->pointToElement.assign(this->local_points, lm_uint_max);
// constexpr UInt Dim = ContainerA::Dim;
// LM_ASSERT(this->positions.size() == this->local_points * Dim,
// "inconsistency " << this->positions.size()
// << " != " << this->local_points * Dim);
// // Loop over all positions of A world
// bool found;
// Real min_dist;
// UInt not_found_counter = 0;
// Real *posA_ptr = &this->positions[0];
// for (UInt indexA = 0; indexA < this->local_points;
// ++indexA, posA_ptr += Dim) {
// // extract the points of pointsBlist which are candidates for association
// VectorView<Dim> positionA(posA_ptr);
// LM_TOIMPLEMENT;
// // auto &&subsetPointsB = this->grid->findSet(positionA);
// found = false;
// min_dist = lm_real_max;
// // loop the candidates
// LM_TOIMPLEMENT;
// // for (auto &&indexB : subsetPointsB) {
// // auto &&positionB = unmatchedB.get(indexB).position0();
// // Real dist = temp_coincidence<Dim>(positionA, positionB);
// // if (dist < min_dist)
// // min_dist = dist;
// // if (coincidence<Dim>(positionA, positionB, this->tolerance)) {
// // this->assoc[indexA] = indexB;
// // ++this->assoc_found;
// // found = true;
// // if (indexA < 10)
// // DUMP("associate point " << positionA << ", " << indexB
// // << " with position " << indexA,
// // DBG_DETAIL);
// // break;
// // }
// // }
// if (!found) {
// ++not_found_counter;
// {
// std::stringstream sstr;
// sstr << not_found_counter << ": couldn't associate point (";
// for (UInt i = 0; i < Dim - 1; ++i)
// sstr << positionA[i] << ", ";
// sstr << positionA[Dim - 1] << "), "
// << "the minimal distance was " << min_dist;
// DUMP(sstr.str(), DBG_INFO);
// }
// }
// }
}
/* -------------------------------------------------------------------------- */
void PointAssociation::clearAll() {
this->assoc_found = 0;
this->pointToElement.clear();
LM_TOIMPLEMENT;
// this->grid.reset();
// pointsAlist->clear();
// pointsBlist->clear();
}
/* -------------------------------------------------------------------------- */
void PointAssociation::filterPointListAForUnmatched() {
LM_TOIMPLEMENT;
// this->pointsAlist->clear();
// for (UInt i = 0; i < this->pointToElement.size(); ++i) {
// if (this->pointToElement[i] != UINT_MAX) {
// DUMP("for container, atom " << i << " becomes "
// << this->pointsAlist->size(),
// DBG_ALL);
// this->pointsAlist->push_back(this->unmatchedA->get(i));
// } else {
// DUMP("atom filtered " << i << " " << this->unmatchedA->get(i),
// DBG_DETAIL);
// }
// }
// this->pointsAlist->copyContainerInfo(unmatchedA);
LM_TOIMPLEMENT;
// this->pointsAlist->copyReleaseInfo(unmatchedA);
}
/* -------------------------------------------------------------------------- */
void PointAssociation::filterPointListBForUnmatched() {
LM_TOIMPLEMENT;
// this->pointsBlist->clear();
// UInt counter = 0;
// for (UInt i = 0; i < this->pointToElement.size(); ++i) {
// if (this->pointToElement[i] != UINT_MAX) {
// DUMP("for container, atom " << i << " becomes "
// << this->pointsBlist->size(),
// DBG_ALL);
// this->pointsBlist->push_back(unmatchedB->get(this->pointToElement[i]));
// this->pointToElement[i] = counter;
// ++counter;
// } else {
// DUMP("atom filtered " << i, DBG_DETAIL);
// }
// }
// this->pointsBlist->getBoundingBox() = unmatchedB->getBoundingBox();
}
/* -------------------------------------------------------------------------- */
void PointAssociation::checkPositionCoherency() {
LM_TOIMPLEMENT;
// try {
// std::stringstream sstr;
// sstr << this->getID() << "-duovector-" << current_step << "-" <<
// current_stage;
// this->getDuoVector().print(sstr.str());
// }
// catch (LibMultiScaleException & e){
// }
// extract the position of local atoms
// if (this->in_group_A)
// this->buildPositions(this->pointsAlist.getOutput());
// constexpr UInt Dim = ContainerA::Dim;
// // transport the positions to B side
// this->distributeVectorA2B("temp_positions", this->positions);
// if (this->in_group_B) {
// UInt point_idx = 0;
// for (auto &&point : this->pointsBlist) {
// auto &&local_pos = point.position0();
// std::stringstream sstr;
// for (UInt i = 0; i < Dim; ++i)
// sstr << this->positions[Dim * point_idx + i] << " ";
// sstr << " and ";
// for (UInt i = 0; i < Dim; ++i)
// sstr << local_pos[i] << " ";
// std::string _tmp = sstr.str();
// DUMP("for point " << point_idx << " compare " << _tmp, DBG_DETAIL);
// for (UInt k = 0; k < Dim; ++k) {
// if (fabs(this->positions[Dim * point_idx + k] - local_pos[k]) > 1e-5)
// {
// LM_FATAL(this->getID()
// << ": inconsistency in the communication scheme "
// << this->positions[Dim * point_idx + k]
// << " != " << local_pos[k] << " " << point_idx);
// }
// }
// ++point_idx;
// }
// }
// DUMP("position coherency OK", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename ContA, typename ContB>
void PointAssociation::updateForMigration(ContA &contA, ContB &contB) {
// constexpr UInt Dim = ContainerA::Dim;
DUMP("update migration for " << this->getID(), DBG_INFO);
DUMP("update local_points " << this->getID(), DBG_INFO);
STARTTIMER("syncMigration resize");
// update local number of points
if (this->in_group_A)
this->local_points = this->pointsAlist.size();
// update local number of points
if (this->in_group_B)
this->local_points = this->pointsBlist.size();
this->positions.resize(spatial_dimension * this->local_points);
STOPTIMER("syncMigration resize");
DUMP("synchMigration " << this->getID(), DBG_INFO);
STARTTIMER("syncMigration duo");
if ((this->in_group_A && contA.hasRefManager()) ||
(this->in_group_B && !contB.hasRefManager()))
this->synchronizeMigration(this->comm_group_A, this->comm_group_B);
if ((this->in_group_B && contB.hasRefManager()) ||
(this->in_group_A && !contA.hasRefManager()))
this->synchronizeMigration(this->comm_group_B, this->comm_group_A);
STOPTIMER("syncMigration duo");
DUMP("update release " << this->getID(), DBG_INFO);
LM_TOIMPLEMENT; // to review
// pointsAlist.setRelease(this->contA.getRelease());
// pointsBlist.setRelease(this->contB.getRelease());
if (check_coherency) {
DUMP("check position coherency " << this->getID(), DBG_INFO);
STARTTIMER("check position coherency");
this->checkPositionCoherency();
STOPTIMER("check position coherency");
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC PointAssociation
This coupler computes the association between point degrees of freedom of
different domains
*/
/* LMHERITATE ComputeTrueDisplacement
*/
void PointAssociation::declareParams() {
DofAssociation::declareParams();
this->addSubParsableObject(true_disp);
/* LMKEYWORD TOL
max distance between two points to consider them coincident
*/
this->parseKeyword("TOL", tolerance);
/* LMKEYWORD PATH
file path to the filteratom file
*/
this->parseKeyword("PATH", path);
/* LMKEYWORD PBC
Specify if periodic boundary conditions should be considered for each
direction
*/
this->parseVectorKeyword("PBC", spatial_dimension, this->pbc);
/* LMKEYWORD CHECK_COHERENCY
Perform a systematic check of the communication scheme.
\textbf{Be careful, it is extremely computationally expensive}
*/
this->parseKeyword("CHECK_COHERENCY", check_coherency, false);
/* LMKEYWORD WEIGHT_PATH
file path to the python script to generate the weights
*/
this->parseKeyword("WEIGHT_PATH", weight_path);
}
/* -------------------------------------------------------------------------- */
// DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(PointAssociation)
// DECLARE_BRIDGING_CONTINUUM_ATOMIC_TEMPLATE(PointAssociation)
__END_LIBMULTISCALE__
c4science · Help