bridging_inline_impl.hh
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Sep 5, 16:23

bridging_inline_impl.hh
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.
	*
	*/

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

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

	DofAssociation::init(contA, contC);
	const UInt Dim = spatial_dimension;

	contMesh = contC;

	pointList.alloc<typename ContA::ContainerSubset>("matched-point-" +
	this->getID());
	meshList.alloc<typename ContC::ContainerSubset>("matched-fe-" +
	this->getID());

	buffer_for_points.acquireContext(contA);
	buffer_for_nodes.acquireContext(contC);

	DUMPBYPROC("selecting DOFs in bridging zone", DBG_INFO_STARTUP, 0);
	if (this->in_group_A) {
	this->buildPointList(contA);
	this->buildPositions<dispatch>(unmatchedPointList.evalOutput());
	}
	if (this->in_group_B)
	this->buildContinuumDOFsList(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<dispatch>(this->unmatchedPointList.evalOutput(),
	this->unmatchedMeshList.evalOutput());
	MPI_Barrier(MPI_COMM_WORLD);

	DUMP(this->getID() << " : exchange points positions", DBG_MESSAGE);
	this->exchangePositions(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<dispatch>(this->local_points,
	this->unmatchedMeshList.evalOutput());
	this->local_points = this->assoc_found;
	DUMP(this->getID() << " : I will manage " << this->local_points
	<< " points",
	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 point assignement",
	DBG_MESSAGE);
	std::vector<std::vector<UInt>> unassociated;
	if (this->in_group_A)
	this->createDuoVectorA("Point", 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) {
	filterPointListForUnmatched<dispatch>(pointList,
	unmatchedPointList.evalOutput());
	}

	if (this->in_group_B) {
	this->filterArray(this->positions);
	this->filterAssoc();
	// this->buildNodeShape<dispatch>(this->meshList);
	}

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

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

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

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

	constexpr UInt Dim = ContC::ContainerSubset::Dim;

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

	auto &unmatched_mesh =
	unmatchedMeshList.evalOutput<typename ContC::ContainerSubset>();
	DUMP("we have found " << unmatched_mesh.getContainerElems().size()
	<< " concerned elements",
	DBG_INFO_STARTUP);

	STARTTIMER("Filling spatial-grid");

	auto contElems = unmatched_mesh.getContainerElems();
	auto contNodes = unmatched_mesh.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);

	auto &grid =
	this->grid.alloc<SpatialGridElem<UInt, Dim>>(cube, this->grid_division);

	UInt nb_elem = 0;
	for (auto &&el : contElems) {
	std::vector<Vector<Dim>> 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();
	LM_ASSERT(mass_node > 0, "invalid mass" << nodes[i] << " " << mass_node);
	#endif

	node_coords.push_back(nd.position0());
	}

	grid.addFiniteElement(nb_elem, node_coords);
	++nb_elem;
	}
	STOPTIMER("Filling spatial-grid");

	DUMP("we have found " << contElems.size() << " concerned elements",
	DBG_INFO_STARTUP);
	DUMP("in average " << grid.getAverageEltByBlock() << " elements par block",
	DBG_INFO_STARTUP);
	}
	/* -------------------------------------------------------------------------- */

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

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

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

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

	template <typename ContMesh>
	void Bridging::buildShapeMatrix(UInt nb_points, ContMesh &unmatchedMeshList) {

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

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

	constexpr UInt Dim = ContMesh::ContainerSubset::Dim;
	auto &grid = this->grid.get<SpatialGridElem<UInt, Dim> &>();

	//****************************************************************
	/* building the association array points <-> elements */
	//****************************************************************
	STARTTIMER("Construction association");
	for (UInt i = 0; i < nb_points; ++i) {

	#ifndef TIMER
	if (i % 100000 == 0)
	DUMP("building association: point " << i << "/" << nb_points, DBG_INFO);
	#endif

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

	this->pointToElement[i] = UINT_MAX;

	for (auto &&j : grid.findSet(x)) {
	auto el = unmatchedMeshList.getContainerElems().get(j);

	DUMP("is point " << 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 point while filtering");
	DUMP("associating point " << i << " and element " << j, DBG_ALL);
	this->pointToElement[i] = j;
	++assoc_found;
	break;
	}
	}
	}

	STOPTIMER("Construction association");
	std::vector<UInt> nb_points_per_element(
	unmatchedMeshList.getContainerElems().size());

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

	auto &meshList =
	this->getOutput<typename ContMesh::ContainerSubset>("meshList");
	/* filter element&node list for the unassociated nodes&elements */
	filterContainerElems<dispatch>(nb_points_per_element, unmatchedMeshList,
	meshList, this->contMesh);

	auto &elements = meshList.getContainerElems();
	auto &nodes = meshList.getContainerNodes();

	smatrix.resize(nodes.size(), assoc_found);

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

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

	auto &&X = this->positions(i);
	auto el = elements.get(this->pointToElement[i]);
	el.computeShapes(shapes, X);
	auto &&node_indices = el.getAlteredConnectivity();

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

	std::vector<UInt> local_nodes(shapes.size());
	for (UInt nd = 0; nd < shapes.size(); ++nd)
	smatrix.insert(node_indices[nd], atom_index) = shapes[nd];

	++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);
	LM_TOIMPLEMENT;
	// 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<dispatch>(this->meshList);
	// 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 ContC1, typename ContC2, typename ContMesh>
	void Bridging::filterContainerElems(std::vector<UInt> &nb_atom_per_element,
	ContC1 &unmatchedMeshList, ContC2 &meshList,
	ContMesh &contMesh) {

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

	for (UInt el = 0; el < nb_atom_per_element.size(); ++el) {
	if (nb_atom_per_element[el] > 0) {

	DUMP("elem " << el << " becomes " << new_t.size(), DBG_ALL);
	new_index.push_back(new_t.size());

	nb_atom_per_element[new_t.size()] = nb_atom_per_element[el];
	new_t.push_back(unmatchedMeshList.getContainerElems().get(el));
	} else
	new_index.push_back(new_t.size());
	}
	nb_atom_per_element.resize(new_t.size());
	meshList.clear();
	// rebuild the element container
	meshList.getContainerElems() = new_t;
	// 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]];
	}
	}

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

	DUMP("We found " << nb << " nodes concerned into "
	<< meshList.getContainerElems().size() << " elements",
	DBG_INFO_STARTUP);
	}
	/* -------------------------------------------------------------------------- */
	template <typename ContA1, typename ContA2>
	void Bridging::filterPointListForUnmatched(ContA1 &pointList,
	ContA2 &unmatchedPointList) {

	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.acquireContext(unmatchedPointList);
	}
	/* -------------------------------------------------------------------------- */

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

	// build right hand side of leastsquare system
	mesh_data.resize(meshList.size(), atomic_data.cols());
	mesh_data = smatrix * atomic_data.matrix();
	cumulPBC(mesh_data);

	LM_TOIMPLEMENT;
	// // solve least square system
	// for (UInt i = 0; i < meshList.size(); ++i) {
	// mesh_data(i) /= node_shape[i];
	// }
	}
	/* -------------------------------------------------------------------------- */

	// template <typename ContC>
	// void Bridging::averageOnElements(ContainerArray<Real> &data_atomic,
	// ContainerArray<Real> &data_mesh,
	// ContC &meshList [[gnu::unused]]) {

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

	// DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(Bridging)
	// DECLARE_BRIDGING_DD_CONTINUUM_TEMPLATE(Bridging)

	template <typename SparseMatrix>
	void Bridging::buildLeastSquareMatrix(SparseMatrix &mat) {
	mat = smatrix * smatrix.transpose();
	cumulPBC(mat);
	}
	/* -------------------------------------------------------------------------- */

	template <typename Matrix> void Bridging::cumulPBC(Matrix &mat) {

	for (auto &&[i1, i2] : local_pbc_pairs) {
	mat.row(i2) += mat.row(i1);
	}

	for (auto &&[i1, i2] : local_pbc_pairs) {
	mat.row(i1) += mat.row(i2);
	}
	}

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

	__END_LIBMULTISCALE__

bridging_inline_impl.hhNo OneTemporaryActions

File Metadata

bridging_inline_impl.hhView Options

Event Timeline

bridging_inline_impl.hh
No OneTemporary
Actions

bridging_inline_impl.hh
View Options