Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/coupling/bridging_inline_impl.hh
e82f2fefbcf8lammps_2020

Last Change
Raw File

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__
c4science · Help