point_association.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jun 6, 09:51

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.buildManual(pointsAlist);
	this->weights.setParam("FILENAME", this->weight_path);
	this->weights.init();
	this->weights.buildManual(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.getArray()[index] * mass;
	++index;
	}
	index = 0;
	for (auto &&point : pointsBlist) {
	Real mass = point.mass();
	weightsB[index] = (1. - weights.getArray()[index]) * mass;
	++index;
	}

	for (UInt index = 0; index < weights.getArray().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__

point_association.ccNo OneTemporaryActions

File Metadata

point_association.ccView Options

Event Timeline

point_association.cc
No OneTemporary
Actions

point_association.cc
View Options