Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88780000
bridging.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 20, 15:46
Size
31 KB
Mime Type
text/x-c++
Expires
Tue, Oct 22, 15:46 (2 d)
Engine
blob
Format
Raw Data
Handle
21819668
Attached To
rLIBMULTISCALE LibMultiScale
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__
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
Bridging<ContainerPoints, ContainerMesh, Dim>::Bridging(const std::string &name,
ContainerPoints &cP,
ContainerMesh &cM)
: LMObject(name), DofAssociation<ContainerPoints, ContainerMesh, Dim>(
name, cP, cM),
unmatchedPointList("unmatched-point-" + name),
pointList("matched-point-" + name),
unmatchedMeshList("unmatched-fe-" + name), meshList("matched-fe-" + name),
node_shape("node_shape-" + name), contPoints(cP), contMesh(cM) {
assoc_found = 0;
smatrix = NULL;
/* registering computes for outer world */
FilterManager::getManager().addObject(unmatchedMeshList);
FilterManager::getManager().addObject(unmatchedPointList);
FilterManager::getManager().addObject(pointList);
FilterManager::getManager().addObject(meshList);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
Bridging<ContainerPoints, ContainerMesh, Dim>::~Bridging() {
if (smatrix)
delete smatrix;
smatrix = NULL;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::init() {
DUMPBYPROC("selecting DOFs in bridging zone", DBG_INFO_STARTUP, 0);
if (this->in_group_A) {
this->buildPointList();
this->buildPositions(unmatchedPointList);
}
if (this->in_group_B)
this->buildContinuumDOFsList();
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();
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);
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, Dim);
this->filterAssoc();
this->buildNodeShape();
}
this->deleteGrid();
DUMP(this->getID() << " : all done for scheme generation", DBG_MESSAGE);
MPI_Barrier(MPI_COMM_WORLD);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildContinuumDOFsList() {
unmatchedMeshList.setParam("GEOMETRY", this->geomID);
unmatchedMeshList.buildManual(contMesh);
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);
if (this->grid)
this->deleteGrid();
this->grid = new SpatialGridElem<UInt, 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) {
RefNode 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 < 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 ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildNodeList() {
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 ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildPointList() {
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 ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildShapeMatrix(
UInt nb_atoms) {
this->pointToElement.assign(nb_atoms, UINT_MAX);
if (unmatchedMeshList->getContainerElems().size() == 0) {
DUMP("elem_rec is empty!", DBG_WARNING);
// delete grid;
return;
}
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<Dim> x;
for (UInt k = 0; k < Dim; ++k)
x[k] = this->positions(i, k);
std::vector<UInt> &subset_elts = this->grid->findSet(x);
std::vector<UInt>::iterator it = subset_elts.begin();
for (it = subset_elts.begin(); it != subset_elts.end(); ++it) {
RefElem 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();
smatrix = new ShapeMatrix<Dim>(meshList.getContainerElems().size());
UInt i = 0;
for (auto &&el : meshList.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;
for (UInt k = 0; k < Dim; ++k)
X[k] = this->positions(i, k);
auto el = meshList.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] = meshList.subIndex2Index(nodes[nd]);
smatrix->fill(atom_index, this->pointToElement[i], nodes, local_nodes, shapes);
++atom_index;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildLocalPBCPairs() {
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 ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildNodeShape() {
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, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::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]];
}
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh,
Dim>::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, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::filterArray(
ContainerArray<Real> &array, UInt stride) {
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);
for (UInt k = 0; k < stride; ++k)
array[stride * index + k] = array[stride * i + k];
++index;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::allocateBuffer(
UInt size, UInt stride) {
if (size > buffer.size()) {
buffer.assign(stride * size, 0);
DUMP("allocating buffer of size " << size * stride, DBG_INFO);
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::interpolatePointData(
std::vector<Real> &data_atom, std::vector<Real> &data_node, UInt stride) {
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, stride);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::cumulPBC(
ContainerArray<Real> &data, UInt stride) {
UInt npairs = local_pbc_pairs.size();
DUMP("Detected " << npairs << " local pairs", DBG_INFO);
for (UInt pair = 0; pair < npairs; ++pair) {
UInt i1 = local_pbc_pairs[pair].first;
UInt i2 = local_pbc_pairs[pair].second;
for (UInt k = 0; k < stride; ++k) {
data[i2 * stride + k] += data[stride * i1 + k];
DUMP("cumul pbc : [" << i2 * stride + k << "] = " << data[i2 * stride + k]
<< " += [" << i2 * stride + k
<< "] = " << data[stride * i1 + k] << " := "
<< data[stride * i1 + k] + data[stride * i2 + k],
DBG_DETAIL);
}
}
for (UInt pair = 0; pair < npairs; ++pair) {
UInt i1 = local_pbc_pairs[pair].first;
UInt i2 = local_pbc_pairs[pair].second;
for (UInt k = 0; k < stride; ++k)
data[stride * i1 + k] = data[stride * i2 + k];
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::leastSquarePointData(
std::vector<Real> &, std::vector<Real> &, UInt, UInt) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::solveLeastSquare(
ContainerArray<Real> &mesh_data, ContainerArray<Real> &atomic_data,
UInt stride) {
// build right hand side of leastsquare system
mesh_data.assign(stride * meshList.size(), 0);
smatrix->buildLeastSquareRHS(mesh_data, atomic_data, stride);
cumulPBC(mesh_data, stride);
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)
for (UInt k = 0; k < stride; ++k) {
DUMP(mesh_data[i * stride + k]
<< " / " << node_shape[i] << " = "
<< mesh_data[i * stride + k] / node_shape[i],
DBG_DETAIL);
mesh_data[i * stride + k] /= node_shape[i];
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::averageOnElements(
ContainerArray<Real> &data_atomic, ContainerArray<Real> &data_mesh,
UInt stride) {
data_mesh.assign(stride * meshList.getContainerElems().size(), 0);
smatrix->averageOnElements(data_atomic, data_mesh, stride);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::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, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::buildLeastSquareMatrix(
math::Matrix &mat) {
smatrix->buildLeastSquareMatrix(&mat);
cumulPBC(mat);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::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, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::copySlaveValues(
std::vector<Real> &v, UInt stride) {
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 k = 0; k < stride; ++k) {
v[i1 * stride + k] = v[i2 * stride + k];
}
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints, ContainerMesh, Dim>::extractShapeMatrix(
math::Matrix &mat) {
smatrix->extractShapeMatrix(&mat);
}
/* -------------------------------------------------------------------------- */
/* Parallel Methods */
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Bridging<DomainA, DomainC, Dim>::communicateBufferAtom2Continuum(
UInt stride, ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->distributeVectorA2B("communicate buffer", buf, stride);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Bridging<DomainA, DomainC, Dim>::communicateBufferContinuum2Atom(
UInt stride, ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->distributeVectorB2A("communicate buffer", buf, stride);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Bridging<DomainA, DomainC, Dim>::synchSumBuffer(
UInt stride, ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->synchronizeVectorBySum("synchsumvector", buf, stride);
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void Bridging<ContainerA, ContainerB, Dim>::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, UInt Dim>
void Bridging<DomainA, DomainC, Dim>::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(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, UInt Dim>
void Bridging<ContainerA, ContainerB, Dim>::declareParams() {
this->addSubParsableObject(unmatchedMeshList);
DofAssociation<ContainerA, ContainerB, Dim>::declareParams();
}
/* -------------------------------------------------------------------------- */
DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(Bridging)
DECLARE_BRIDGING_DD_CONTINUUM_TEMPLATE(Bridging)
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment