bridging.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Jun 1, 13:05

bridging.cc
View Options

	/**
	* @file bridging.cc
	*
	* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
	*
	* @date Fri Jul 11 15:47:44 2014
	*
	* @brief Bridging object between atomistic and finite elements
	*
	* @section LICENSE
	*
	* Copyright INRIA and CEA
	*
	* The LibMultiScale is a C++ parallel framework for the multiscale
	* coupling methods dedicated to material simulations. This framework
	* provides an API which makes it possible to program coupled simulations
	* and integration of already existing codes.
	*
	* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
	* within ScAlApplix team and CEA/DPTA Ile de France.
	* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
	* within the LSMS/ENAC laboratory.
	*
	* This software is governed by the CeCILL-C license under French law and
	* abiding by the rules of distribution of free software. You can use,
	* modify and/ or redistribute the software under the terms of the CeCILL-C
	* license as circulated by CEA, CNRS and INRIA at the following URL
	* "http://www.cecill.info".
	*
	* As a counterpart to the access to the source code and rights to copy,
	* modify and redistribute granted by the license, users are provided only
	* with a limited warranty and the software's author, the holder of the
	* economic rights, and the successive licensors have only limited
	* liability.
	*
	* In this respect, the user's attention is drawn to the risks associated
	* with loading, using, modifying and/or developing or reproducing the
	* software by the user in light of its specific status of free software,
	* that may mean that it is complicated to manipulate, and that also
	* therefore means that it is reserved for developers and experienced
	* professionals having in-depth computer knowledge. Users are therefore
	* encouraged to load and test the software's suitability as regards their
	* requirements in conditions enabling the security of their systems and/or
	* data to be ensured and, more generally, to use and operate it in the
	* same conditions as regards security.
	*
	* The fact that you are presently reading this means that you have had
	* knowledge of the CeCILL-C license and that you accept its terms.
	*
	*/

	//#define TIMER
	/* -------------------------------------------------------------------------- */
	#include "bridging.hh"
	#include "compute_extract.hh"
	#include "factory_multiscale.hh"
	#include "filter_geometry.hh"
	#include "geometry_manager.hh"
	#include "lib_bridging.hh"
	#include "lm_common.hh"
	#include "trace_atom.hh"
	/* -------------------------------------------------------------------------- */
	__BEGIN_LIBMULTISCALE__

	Bridging::Bridging(const std::string &name)
	: LMObject(name), DofAssociation(name),
	unmatchedPointList("unmatched-point-" + name),
	unmatchedMeshList("unmatched-fe-" + name),
	node_shape("node_shape-" + name) {
	assoc_found = 0;

	unmatchedPointList.setCommGroup(this->comm_group_A);
	unmatchedMeshList.setCommGroup(this->comm_group_B);
	/* registering computes for outer world */
	FilterManager::getManager().addObject(unmatchedMeshList);
	FilterManager::getManager().addObject(unmatchedPointList);
	FilterManager::getManager().addObject(pointList);
	FilterManager::getManager().addObject(meshList);

	is_in_continuum = this->comm_group_B.amIinGroup();
	is_in_atomic = this->comm_group_A.amIinGroup();
	}
	/* -------------------------------------------------------------------------- */
	Bridging::~Bridging() {}

	/* -------------------------------------------------------------------------- */
	template <typename types, typename F, typename... Args>
	decltype(auto) decorate(F &f, Args... args) {
	call_compute(*this, [&](auto &... args) { f(args...); }, args...);
	//
	}

	template <typename ContA, typename ContC>
	void Bridging::init(ContA &contA, ContC &contC) {

	pointList = std::make_shared<FilterCompatibility<ContA>>("matched-point-" +
	this->getID());
	pointList->setCommGroup(this->comm_group_A);

	meshList = std::make_shared<FilterCompatibility<ContC>>("matched-fe-" +
	this->getID());
	meshList->setCommGroup(this->comm_group_B);

	DUMPBYPROC("selecting DOFs in bridging zone", DBG_INFO_STARTUP, 0);
	if (this->in_group_A) {
	this->buildPointList(contA);
	this->buildPositions(unmatchedPointList);
	}
	if (this->in_group_B)
	this->buildContinuumDOFsList(contA, contC);

	DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP);
	MPI_Barrier(MPI_COMM_WORLD);
	DUMP(this->getID() << " : exchange coarse geometries", DBG_MESSAGE);
	this->exchangeGeometries(this->unmatchedPointList, this->unmatchedMeshList);
	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : exchange points positions", DBG_MESSAGE);
	this->exchangePositions(ContA::Dim);
	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : build shape matrix", DBG_MESSAGE);
	if (this->in_group_B) {
	DUMP(this->getID() << " : build shape matrix (parsing "
	<< this->local_points << " points by position)",
	DBG_INFO);
	this->buildShapeMatrix(this->local_points, contA, contC);
	this->local_points = this->assoc_found;
	DUMP(this->getID() << " : I will manage " << this->local_points << " atoms",
	DBG_INFO);
	}
	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : ExchangeAssociation", DBG_MESSAGE);
	ContainerArray<UInt> managed;
	this->exchangeAssociationInformation(managed);

	if (this->in_group_A)
	for (UInt i = 0; i < this->nb_zone_B; ++i)
	DUMP("second test Do I communicate with " << i << " : "
	<< this->com_with[i],
	DBG_INFO);

	if (this->in_group_B)
	for (UInt i = 0; i < this->nb_zone_A; ++i)
	DUMP("second test Do I communicate with " << i << " : "
	<< this->com_with[i],
	DBG_INFO);

	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID()
	<< " : generate duo vector and detect double atom assignement",
	DBG_MESSAGE);
	std::vector<std::vector<UInt>> unassociated;
	if (this->in_group_A)
	this->createDuoVectorA("Atom", managed, unassociated);

	DUMP(this->getID() << " : doing exchange rejected", DBG_MESSAGE);
	this->exchangeRejectedContinuumOwners(unassociated);

	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : doing filter rejected", DBG_MESSAGE);
	this->filterRejectedContinuumOwners(unassociated);
	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : filtering position & association vector",
	DBG_MESSAGE);

	if (this->in_group_A) {
	this->filterPointListForUnmatched();
	}

	if (this->in_group_B) {
	this->filterArray(this->positions);
	this->filterAssoc();
	this->buildNodeShape();
	}

	// LM_TOIMPLEMENT
	// this->grid.reset();

	DUMP(this->getID() << " : all done for scheme generation", DBG_MESSAGE);

	MPI_Barrier(MPI_COMM_WORLD);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContA, typename ContC>
	void Bridging::buildContinuumDOFsList(ContA &contA, ContC &contMesh) {

	unmatchedMeshList.setParam("GEOMETRY", this->geomID);
	unmatchedMeshList.buildManual(contMesh);

	LM_TOIMPLEMENT;
	// DUMP("we have found " << unmatchedMeshList.getContainerElems().size()
	// << " concerned elements",
	// DBG_INFO_STARTUP);

	// STARTTIMER("Filling spatial-grid");

	// auto contElems = unmatchedMeshList.getContainerElems();
	// auto contNodes = unmatchedMeshList.getContainerNodes();

	// std::vector<UInt> nodes;
	// Geometry &geom = *GeometryManager::getManager().getGeometry(this->geomID);
	// Cube cube = geom.getBoundingBox();
	// DUMP("geometry of the grid => \n" << cube, DBG_INFO_STARTUP);
	// this->grid.reset();
	LM_TOIMPLEMENT;
	// this->grid = new SpatialGridElem<UInt, ContainerPoints::Dim>(
	// cube, this->grid_division);

	// UInt nb_elem = 0;
	// for (auto &&el : contElems) {
	// std::vector<Real> node_coords;
	// auto &&nodes = el.globalIndexes();
	// UInt nb = nodes.size();

	// for (UInt i = 0; i < nb; ++i) {
	// auto nd = el.getNode(i);
	// #ifndef LM_OPTIMIZED
	// Real mass_node = nd.mass();
	// #endif
	// LM_ASSERT(mass_node > 0, "invalid mass" << nodes[i] << " " <<
	// mass_node);

	// auto X = nd.position0();
	// for (UInt i = 0; i < ContainerPoints::Dim; ++i)
	// node_coords.push_back(X[i]);
	// }

	// this->grid->addFiniteElement(nb_elem, node_coords);
	// ++nb_elem;
	// }
	// STOPTIMER("Filling spatial-grid");

	// DUMP("we have found " << nb_elem << " concerned elements",
	// DBG_INFO_STARTUP); DUMP("in average " <<
	// this->grid->getAverageEltByBlock()
	// << " elements par block",
	// DBG_INFO_STARTUP);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerMesh>
	void Bridging::buildNodeList([[gnu::unused]] ContainerMesh &contMesh) {

	LM_TOIMPLEMENT;
	// meshList->computeAlteredConnectivity(contMesh.getContainerNodes());
	UInt nb = meshList->size();
	if (!nb) {
	DUMP("We found no nodes in bridging zone", DBG_INFO_STARTUP);
	return;
	}

	LM_TOIMPLEMENT;
	// DUMP("We found " << nb << " nodes concerned into "
	// << meshList->getContainerElems().size() << " elements",
	// DBG_INFO_STARTUP);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContA> void Bridging::buildPointList(ContA &contPoints) {

	unmatchedPointList.setParam("GEOMETRY", this->geomID);
	unmatchedPointList.buildManual(contPoints);
	this->local_points = unmatchedPointList.size();

	if (!this->local_points)
	DUMP("We found no atom in the bridging zone", DBG_INFO);

	DUMP("We found " << this->local_points << " atom in the bridging zone",
	DBG_INFO_STARTUP);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContA, typename ContB>
	void Bridging::buildShapeMatrix(UInt nb_atoms, ContA &contA,
	ContB &unmatchedMeshList) {

	this->pointToElement.assign(nb_atoms, UINT_MAX);

	if (unmatchedMeshList.getContainerElems().size() == 0) {
	DUMP("elem_rec is empty!", DBG_WARNING);
	LM_TOIMPLEMENT;
	// need to free the output produced by the component
	// delete grid;
	return;
	}

	constexpr UInt Dim = ContA::Dim;
	auto &grid = this->getCastedOutput<SpatialGridElem<UInt, Dim>>("grid");
	std::vector<UInt> nb_atoms_per_element;
	nb_atoms_per_element.resize(unmatchedMeshList.getContainerElems().size());

	/* construction du array d'assoc atoms elem */

	STARTTIMER("Construction association");

	for (UInt i = 0; i < nb_atoms; ++i) {

	#ifndef TIMER
	if (i % 100000 == 0)
	DUMP("construction de l'association - atome " << i << "/" << nb_atoms,
	DBG_INFO);
	#endif

	Vector<ContA::Dim> x;
	x = this->positions(i);

	std::vector<UInt> &subset_elts = grid->findSet(x);
	std::vector<UInt>::iterator it = subset_elts.begin();

	for (it = subset_elts.begin(); it != subset_elts.end(); ++it) {

	auto el = unmatchedMeshList->getContainerElems().get(*it);
	UInt j = *it;

	DUMP("is atom " << i << " within element " << j << "?", DBG_ALL);

	if (el.contains(x)) {
	LM_ASSERT(this->pointToElement[i] == UINT_MAX,
	"this should not happen. "
	<< "I found twice the same atom while filtering");
	DUMP("associating atom " << i << " and element " << j, DBG_ALL);
	this->pointToElement[i] = j;
	++assoc_found;
	break;
	}
	}

	if (it == subset_elts.end()) {
	this->pointToElement[i] = UINT_MAX;
	}
	}

	STOPTIMER("Construction association");

	/* build number of atoms per element */
	for (UInt i = 0; i < nb_atoms; ++i) {
	if (this->pointToElement[i] != UINT_MAX)
	++nb_atoms_per_element[this->pointToElement[i]];
	}

	/* filter element list and node list for the unassociated node and elements */
	filterContainerElems(nb_atoms_per_element);
	buildNodeList(unmatchedMeshList);

	auto &elements = meshList->getArray<ContB::ContainerElems::Ref>();
	smatrix = std::make_shared<ShapeMatrix>(elements.getContainerElems().size());

	UInt i = 0;
	for (auto &&el : elements.getContainerElems()) {
	DUMP("declare sub matrix number " << i << " which will be of size "
	<< nb_atoms_per_element[i] << "x"
	<< el.nNodes(),
	DBG_ALL);
	LM_ASSERT(nb_atoms_per_element[i] > 0, "this should not happen any more");
	smatrix->setSub(i, nb_atoms_per_element[i], el.nNodes());
	++i;
	}
	// Allocate the shapematrix structure (contiguous)
	smatrix->allocate();

	/* set the indirection (alloc tab) plus the shapes */
	std::vector<Real> shapes;

	UInt atom_index = 0;

	for (i = 0; i < nb_atoms; ++i) {
	if (this->pointToElement[i] == UINT_MAX)
	continue;

	Vector<Dim> X;
	X = this->positions(i);

	auto el = elements.getContainerElems().get(this->pointToElement[i]);
	el.computeShapes(shapes, X);
	auto &&nodes = el.globalIndexes();

	IF_TRACED(X, "checking for atom in trouble");

	std::vector<UInt> local_nodes(shapes.size());
	for (UInt nd = 0; nd < shapes.size(); ++nd)
	local_nodes[nd] = elements.subIndex2Index(nodes[nd]);

	smatrix->fill(atom_index, this->pointToElement[i], nodes, local_nodes,
	shapes);
	++atom_index;
	}
	}
	/* -------------------------------------------------------------------------- */
	template <typename ContC> void Bridging::buildLocalPBCPairs(ContC &meshList) {

	UInt npairs = pbc_pairs.size();
	for (UInt pair = 0; pair < npairs; ++pair) {
	UInt ind1 = pbc_pairs[pair].first;
	UInt ind2 = pbc_pairs[pair].second;
	UInt i1 = meshList.subIndex2Index(ind1);
	UInt i2 = meshList.subIndex2Index(ind2);
	if (i1 == UINT_MAX \|\| i2 == UINT_MAX)
	continue;

	std::pair<UInt, UInt> p((UInt)(i1), (UInt)(i2));
	local_pbc_pairs.push_back(p);
	}
	DUMP("Detected " << local_pbc_pairs.size() << " local pairs", DBG_MESSAGE);
	smatrix->alterMatrixIndexesForPBC(local_pbc_pairs);
	}
	/* -------------------------------------------------------------------------- */
	template <typename ContC> void Bridging::buildNodeShape(ContC &meshList) {

	if (this->local_points == 0)
	return;

	UInt nb_coupled_nodes = meshList.size();
	node_shape.assign(nb_coupled_nodes, 0);

	LM_ASSERT(smatrix, "shapematrix object was not allocated "
	<< " eventhough local atoms were detected");
	smatrix->buildNodeShape(node_shape);

	buildLocalPBCPairs();
	cumulPBC(node_shape);

	#ifndef LM_OPTIMIZED

	for (UInt i = 0; i < nb_coupled_nodes; ++i) {
	if (node_shape[i] <= 1e-1) {
	DUMP(i << " " << node_shape[i], DBG_MESSAGE);
	LM_FATAL("Please check your geometry:"
	<< " one or more nodes has zero nodal shape values."
	<< " Changing the FE element size may resolve this issue!"
	<< std::endl);
	}
	}

	#endif // LM_OPTIMIZED

	FilterManager::getManager().addObject(node_shape);
	}
	/* -------------------------------------------------------------------------- */
	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::filterContainerElems(
	std::vector<UInt> &nb_atome_par_element) {

	typename ContainerMeshSubset::ContainerElems new_t(this->getID() + ":elems");
	std::vector<int> new_index;

	for (UInt i = 0; i < nb_atome_par_element.size(); ++i) {
	DUMP("filtering old index element " << i, DBG_ALL);
	if (nb_atome_par_element[i] > 0) {

	DUMP("for container, elem " << i << " becomes " << new_t.size(), DBG_ALL);
	new_index.push_back(new_t.size());

	DUMP("maintenant nb_atome_par_element["
	<< new_t.size() << "]= " << nb_atome_par_element[i]
	<< " (i=" << i << ")",
	DBG_ALL);
	nb_atome_par_element[new_t.size()] = nb_atome_par_element[i];
	new_t.push_back(unmatchedMeshList->getContainerElems().get(i));
	} else
	new_index.push_back(new_t.size());
	}
	nb_atome_par_element.resize(new_t.size());
	meshList.clear();

	// rebuild the element container
	for (UInt i = 0; i < new_t.size(); ++i)
	meshList.getContainerElems().push_back(new_t.get(i));
	// rebuild the association list
	for (UInt i = 0; i < this->pointToElement.size(); ++i) {
	if (this->pointToElement[i] != UINT_MAX) {
	this->pointToElement[i] = new_index[this->pointToElement[i]];
	}
	}
	}
	/* -------------------------------------------------------------------------- */

	void Bridging::filterPointListForUnmatched() {

	pointList.clear();

	for (UInt i = 0; i < this->pointToElement.size(); ++i)
	if (this->pointToElement[i] != UINT_MAX) {
	DUMP("for container, atom " << i << " becomes " << pointList.size(),
	DBG_ALL);
	pointList.push_back(unmatchedPointList->get(i));
	} else {
	DUMP("atom filtered " << i << " " << unmatchedPointList->get(i),
	DBG_DETAIL);
	}
	pointList.copyContainerInfo(unmatchedPointList);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::filterArray(
	ContainerArray<Real> &array) {
	UInt index = 0;
	for (UInt i = 0; i < this->pointToElement.size(); ++i)
	if (this->pointToElement[i] != UINT_MAX) {
	DUMP("for container, atom " << i << " becomes " << index, DBG_ALL);
	LM_ASSERT(i < array.size() && index < array.size(),
	"overflow detected: nmax = " << array.size() << " , index = "
	<< index << " and i = " << i);

	array(index) = array(i);
	++index;
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::interpolatePointData(
	ContainerArray<Real> &data_atom, ContainerArray<Real> &data_node) {

	if (this->local_points == 0) {
	DUMP("Warning : atomic container is empty:"
	<< " cannot interpolate mesh fields on control points (check "
	"geometries ?!)",
	DBG_WARNING);
	return;
	}
	smatrix->interpolate(data_atom, data_node);
	}

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

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::cumulPBC(
	ContainerArray<Real> &data) {

	UInt npairs [[gnu::unused]] = local_pbc_pairs.size();
	DUMP("Detected " << npairs << " local pairs", DBG_INFO);
	for (auto &&pair : local_pbc_pairs) {
	UInt i1 = pair.first;
	UInt i2 = pair.second;

	data(i2) += data(i1);
	}
	for (auto &&pair : local_pbc_pairs) {
	UInt i1 = pair.first;
	UInt i2 = pair.second;
	data(i1) = data(i2);
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::leastSquarePointData(
	ContainerArray<Real> &, ContainerArray<Real> &, UInt) {
	LM_TOIMPLEMENT;
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::solveLeastSquare(
	ContainerArray<Real> &mesh_data, ContainerArray<Real> &atomic_data) {

	// build right hand side of leastsquare system
	mesh_data.resize(meshList.size(), atomic_data.cols());
	smatrix->buildLeastSquareRHS(mesh_data, atomic_data);
	cumulPBC(mesh_data);

	LM_ASSERT(node_shape.size() != 0,
	"node_shape is empty this should not happen, check geometry");
	// solve least square system
	for (UInt i = 0; i < meshList.size(); ++i) {
	mesh_data(i) /= node_shape[i];
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::averageOnElements(
	ContainerArray<Real> &data_atomic, ContainerArray<Real> &data_mesh) {

	data_mesh.assign(meshList.getContainerElems().rows(),
	meshList.getContainerElems().cols());
	smatrix->averageOnElements(data_atomic, data_mesh);
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::clearAll() {
	if (smatrix)
	delete smatrix;
	smatrix = NULL;
	meshList.getContainerElems().clear();
	// nodeList.clear();
	assoc_found = 0;
	this->pointToElement.clear();
	pointList.clear();
	}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::buildLeastSquareMatrix(
	math::Matrix &mat) {
	smatrix->buildLeastSquareMatrix(&mat);
	cumulPBC(mat);
	}

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

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::cumulPBC(math::Matrix &mat) {

	UInt npairs = local_pbc_pairs.size();

	for (UInt pair = 0; pair < npairs; ++pair) {
	UInt i1 = local_pbc_pairs[pair].first;
	UInt i2 = local_pbc_pairs[pair].second;
	for (UInt j = 0; j < mat.n(); ++j) {
	mat(i2, j) += mat(i1, j);
	}
	}

	for (UInt pair = 0; pair < npairs; ++pair) {
	UInt i1 = local_pbc_pairs[pair].first;
	for (UInt j = 0; j < mat.n(); ++j) {
	if (i1 == j)
	mat(i1, j) = 1.;
	else
	mat(i1, j) = 0.;
	}
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::copySlaveValues(
	ContainerArray<Real> &v) {

	for (auto &&pair : local_pbc_pairs) {
	UInt i1 = pair.first;
	UInt i2 = pair.second;
	v(i1) = v(i2);
	}
	}

	/* -------------------------------------------------------------------------- */
	template <typename ContainerPoints, typename ContainerMesh>
	void Bridging<ContainerPoints, ContainerMesh>::extractShapeMatrix(
	math::Matrix &mat) {
	smatrix->extractShapeMatrix(&mat);
	}

	/* -------------------------------------------------------------------------- */
	/* Parallel Methods */
	/* -------------------------------------------------------------------------- */

	template <typename DomainA, typename DomainC>
	void Bridging<DomainA, DomainC>::communicateBufferAtom2Continuum(
	ContainerArray<Real> &buf) {

	if (this->local_points != 0) {
	this->distributeVectorA2B("communicate buffer", buf);
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename DomainA, typename DomainC>
	void Bridging<DomainA, DomainC>::communicateBufferContinuum2Atom(
	ContainerArray<Real> &buf) {

	if (this->local_points != 0) {
	this->distributeVectorB2A("communicate buffer", buf);
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename DomainA, typename DomainC>
	void Bridging<DomainA, DomainC>::synchSumBuffer(ContainerArray<Real> &buf) {

	if (this->local_points != 0) {
	this->synchronizeVectorBySum("synchsumvector", buf);
	}
	}
	/* -------------------------------------------------------------------------- */

	template <typename ContainerA, typename ContainerB>
	void Bridging<ContainerA, ContainerB>::filterRejectedContinuumOwners(
	std::vector<std::vector<UInt>> &unassociated_atoms) {

	if (this->comm_group_A == this->comm_group_B)
	return;

	if (this->in_group_B) {
	DUMP("local points", DBG_MESSAGE);
	if (this->local_points == 0)
	return;

	UInt offset = 0;
	// recompute indexes of continuum atoms
	std::vector<int> current_indexes(this->pointToElement.size());
	{
	UInt indirection = 0;
	for (UInt at = 0; at < this->pointToElement.size(); ++at) {
	current_indexes[at] = -1;
	if (this->pointToElement[at] == UINT_MAX)
	continue;
	current_indexes[at] = indirection;
	++indirection;
	}
	LM_ASSERT(indirection == this->local_points, "this should not happend "
	<< indirection << " "
	<< this->local_points);
	}
	// remove unassociated atoms
	{
	offset = 0;
	UInt cpt = 0;
	for (UInt i = 0; i < this->nb_zone_A; ++i) {
	if (this->com_with[i]) {
	DUMP("atomic processor " << i << " unassociated "
	<< unassociated_atoms[i].size() << " atoms ",
	DBG_INFO_STARTUP);
	UInt size = unassociated_atoms[i].size();
	for (UInt j = 0; j < size; ++j) {
	// get element that owned the atom
	UInt at_index = offset + unassociated_atoms[i][j];
	LM_ASSERT(this->pointToElement[at_index] != UINT_MAX,
	"this should not happend !! " << at_index);
	LM_ASSERT(current_indexes[at_index] != -1,
	"this should not happend !! " << at_index);

	UInt el = this->pointToElement[at_index];
	this->smatrix->removeAtom(current_indexes[at_index], el);
	this->pointToElement[at_index] = UINT_MAX;
	++cpt;
	DUMP("unassociated atom at position "
	<< at_index << " whereas offset = " << offset
	<< " next offset is "
	<< offset + this->nb_points_per_proc[i],
	DBG_INFO_STARTUP);
	}
	offset += this->nb_points_per_proc[i];
	}
	}
	if (cpt) {
	DUMP("removed " << cpt
	<< " atoms from unassociation. Now local atoms is "
	<< this->local_points,
	DBG_MESSAGE);
	this->local_points -= cpt;
	DUMP("and becomes " << this->local_points, DBG_MESSAGE);
	}
	}
	// little check
	{
	UInt cpt_tmp = 0;
	for (UInt i = 0; i < this->pointToElement.size(); ++i)
	if (this->pointToElement[i] != UINT_MAX)
	cpt_tmp++;

	if (cpt_tmp != this->local_points)
	LM_FATAL("biip this should not happend " << cpt_tmp << " "
	<< this->local_points);
	}
	// now renumber indirection in shapematrix
	{
	UInt indirection = 0;
	for (UInt at = 0; at < this->pointToElement.size(); ++at) {
	if (this->pointToElement[at] == UINT_MAX)
	continue;
	LM_ASSERT(current_indexes[at] != -1, "This should not happend");
	this->smatrix->changeAtomIndirection(this->pointToElement[at],
	current_indexes[at], indirection);
	++indirection;
	}
	this->smatrix->swapAtomIndirections();
	}
	// //now i set the duo vector
	this->createDuoVectorB("FE");
	// std::stringstream sstr;
	// sstr << "duoFE-" << this->getID();
	// DuoDistributedVector & duo =
	// this->createDuoVector(this->local_points,this->total_points);

	// UInt counter = 0;
	// offset = 0;
	// for (UInt i = 0 ; i < this->nb_zone_A ; ++i)
	// if (this->com_with[i]){
	// for (UInt j = 0 ; j < this->nb_points_per_proc[i] ; ++j)
	// if (this->pointToElement[j+offset] != UINT_MAX){
	// duo.setDuoProc(counter,i);
	// ++counter;
	// }
	// DUMP("compteur is " << counter << " after treating atoms from " << i
	// << " offset is " << offset,DBG_INFO_STARTUP);
	// offset += this->nb_points_per_proc[i];
	// }
	}
	}

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

	template <typename DomainA, typename DomainC>
	void Bridging<DomainA, DomainC>::updateForMigration() {

	DUMPBYPROC("update migration for " << this->getID(), DBG_INFO, 0);

	STARTTIMER("syncMigration resize");
	// update local number of atoms in atomic part
	if (this->in_group_A) {
	this->local_points = this->pointList.size();
	this->positions.resize(DomainA::Dim * this->local_points);
	}
	STOPTIMER("syncMigration resize");

	STARTTIMER("syncMigration duo");
	this->synchronizeMigration(this->comm_group_A, this->comm_group_B);
	STOPTIMER("syncMigration duo");

	// pointList.setRelease(this->contA.getRelease());
	// meshList.setRelease(this->contB.getRelease());
	}
	/* -------------------------------------------------------------------------- */

	/* LMDESC Bridging
	This class implements a bridging zone where
	point and elements are associated. \\ \\

	For debugging pruposes, several set of DOf are registered to the central
	system:

	\begin{itemize}
	\item unmatched-fe-\${couplerID} : the set of nodes and elements
	that are contained in the geometry provided by GEOMETRY keyword.
	\item unmatched-point-\${couplerID} : the set of points
	that are contained in the geometry provided by GEOMETRY keyword.
	\item matched-fe-\${couplerID} : the set of nodes and elements
	that are containing at least one point of unmatched-point-\${couplerID}.
	\item matched-point-\${couplerID} : the set of points
	that are contained in at least one element of unmatched-fe-\${couplerID}.
	\end{itemize}
	*/

	/* LMHERITANCE filter_geometry dof_association */

	template <typename ContainerA, typename ContainerB>
	void Bridging<ContainerA, ContainerB>::declareParams() {

	this->addSubParsableObject(unmatchedMeshList);

	DofAssociation<ContainerA, ContainerB>::declareParams();
	}
	/* -------------------------------------------------------------------------- */

	// template <typename DomainA, typename DomainC>
	// void Bridging<DomainA, DomainC>::setAtomicDOFs(FieldType field_code) {
	// if (this->local_points == 0) {
	// DUMP("Warning : atomic container is empty: "
	// << "cannot proceed atomic projection (check geometries ?!)",
	// DBG_WARNING);
	// return;
	// }
	// if (this->buffer.size() == 0) {
	// LM_FATAL("function InterpolateAtomicDOFs should be called first!
	// exiting");
	// }

	// UInt at_index = 0;
	// constexpr UInt Dim = DomainA::Dim;
	// for (auto &&at : this->pointList) {
	// at.force() = Vector<Dim>::Zero();
	// if (field_code == _velocity)
	// at.velocity() = VectorView<Dim>(&this->buffer[at_index * Dim]);
	// if (field_code == _displacement) {
	// at.position() =
	// at.position0() + VectorView<Dim>(&this->buffer[at_index * Dim]);
	// // if (at.position(0) >= xmax) at.position(0) -= xsize;
	// // else if (at.position(0) < xmin) at.position(0) += xsize;
	// //****************************************************************
	// // Guillaume needs to include PBC treatment
	// //****************************************************************
	// at.velocity() = Vector<Dim>::Zero();
	// }
	// ++at_index;
	// }
	// }

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

	DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(Bridging)
	DECLARE_BRIDGING_DD_CONTINUUM_TEMPLATE(Bridging)

	__END_LIBMULTISCALE__

bridging.ccNo OneTemporaryActions

File Metadata

bridging.ccView Options

Event Timeline

bridging.cc
No OneTemporary
Actions

bridging.cc
View Options