point_association.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Jun 22, 22:21

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__

	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	PointAssociation<ContainerA, ContainerB>::PointAssociation(
	const std::string &name, ContainerA &contA, ContainerB &contB)
	: LMObject(name), DofAssociation<ContainerA, ContainerB>(name, contA,
	contB),
	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()) {

	if (this->in_group_A && this->contA.hasRefManager())
	this->contA.getRefManager()->addSubSet(*this->pointsAlist);

	if (this->in_group_B && this->contB.hasRefManager())
	this->contB.getRefManager()->addSubSet(*this->pointsBlist);

	/* registering computes for outer world */
	LM_TOIMPLEMENT;
	// FilterManager::getManager().addObject(&this->pointsAlist, false);
	// FilterManager::getManager().addObject(&this->pointsBlist, false);
	}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	PointAssociation<ContainerA, ContainerB>::~PointAssociation() {}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::init() {

	constexpr UInt Dim = ContainerA::Dim;

	if (this->contA.hasRefManager() && this->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(this->contA, this->unmatchedA);
	this->local_points = this->unmatchedA->size();
	this->buildPositions(this->unmatchedA.getOutput());
	}
	if (this->in_group_B)
	this->buildPointList(this->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)
	this->field_buffer.assign(Dim * (this->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 && this->contA.hasRefManager()) {
	try {
	this->contA.getRefManager()->attachObject(this->getDuoVector(),
	*this->pointsAlist);
	DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
	} catch (LibMultiScaleException &e) {
	}
	}
	if (this->in_group_B && this->contB.hasRefManager()) {
	try {
	this->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();
	}
	}

	/* -------------------------------------------------------------------------- */

	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::computeWeights() {

	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(this->pointsAlist->size(), 0.);
	weightsB.assign(this->pointsAlist->size(), 0.);

	// compute nominator and denominator
	UInt index = 0;
	for (auto &&point : this->pointsAlist) {
	Real mass = point.mass();
	weightsA[index] = weights.getArray()[index] * mass;
	++index;
	}
	index = 0;
	for (auto &&point : this->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 ContainerA, typename ContainerB>
	template <typename container, typename PointList>
	void PointAssociation<ContainerA, ContainerB>::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);
	}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::fillGrid() {
	// 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;
	}
	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::associate() {

	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);
	}
	}
	}
	}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::clearAll() {
	this->assoc_found = 0;
	this->pointToElement.clear();
	this->grid.reset();
	pointsAlist->clear();
	pointsBlist->clear();
	}

	/* -------------------------------------------------------------------------- */

	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::filterPointListAForUnmatched() {

	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);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::filterPointListBForUnmatched() {

	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();
	}

	/* -------------------------------------------------------------------------- */

	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::checkPositionCoherency() {

	// 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 ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::updateForMigration() {
	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(Dim * this->local_points);
	STOPTIMER("syncMigration resize");

	DUMP("synchMigration " << this->getID(), DBG_INFO);
	STARTTIMER("syncMigration duo");
	if ((this->in_group_A && this->contA.hasRefManager()) \|\|
	(this->in_group_B && !this->contB.hasRefManager()))
	this->synchronizeMigration(this->comm_group_A, this->comm_group_B);
	if ((this->in_group_B && this->contB.hasRefManager()) \|\|
	(this->in_group_A && !this->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
	*/

	template <typename ContainerA, typename ContainerB>
	void PointAssociation<ContainerA, ContainerB>::declareParams() {
	DofAssociation<ContainerA, ContainerB>::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