Homec4science
Diffusion LibMultiScale (features/pybind) Browse

libmultiscale/src/coupling/bridging.cc
b05d5f2aba78features/pybind

Last Change
Raw File

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 "lm_common.hh"
#include "bridging.hh"
#include "lib_bridging.hh"
#include "filter_geometry.hh"
#include "compute_extract.hh"
#include "compute_extract.cc"
#include "filter_manager.hh"
#include "geometry_manager.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):
DofAssociation<ContainerPoints,ContainerMesh,Dim>(name,cP,cM),
unmatchedPointList("unmatched-point-" + name,Dim),
pointList("matched-point-" + name,Dim),
unmatchedMeshList("unmatched-fe-" + name,Dim),
meshList("matched-fe-" + name,Dim),
node_shape("node_shape-" + name,1),
contPoints(cP),
contMesh(cM)
{
assoc_found = 0;
smatrix = NULL;
/* registering computes for outer world */
FilterManager::getManager().addObject(&unmatchedMeshList,false);
FilterManager::getManager().addObject(&unmatchedPointList,false);
FilterManager::getManager().addObject(&pointList,false);
FilterManager::getManager().addObject(&meshList,false);
}
/* -------------------------------------------------------------------------- */
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.getOutput());
}
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);
std::vector<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.build(contMesh);
DUMP("we have found " << unmatchedMeshList.getContainerElems().nbElem()
<< " concerned elements",DBG_INFO_STARTUP);
STARTTIMER("Filling spatial-grid");
typename ContainerMeshSubset::ContainerElems contElems =
unmatchedMeshList.getContainerElems();
typename ContainerMeshSubset::ContainerNodes contNodes =
unmatchedMeshList.getContainerNodes();
IteratorElemsSubset it = contElems.getIterator();
UInt nb_elem = 0;
std::vector< UInt > nodes;
Real X[Dim];
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);
for(RefElem el = it.getFirst();!it.end();el = it.getNext(), ++nb_elem) {
std::vector<Real> node_coords;
el.globalIndexes(nodes);
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);
nd.getPositions0(X);
for (UInt i = 0; i < Dim; ++i)
node_coords.push_back(X[i]);
}
this->grid->addFiniteElement(nb_elem,node_coords);
}
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.nbElem();
if (!nb){
DUMP("We found no nodes in bridging zone",DBG_INFO_STARTUP);
return;
}
DUMP("We found " << nb << " nodes concerned into "
<< meshList.getContainerElems().nbElem() << " elements",DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints,ContainerMesh,Dim>::
buildPointList(){
unmatchedPointList.setParam("GEOMETRY", this->geomID);
unmatchedPointList.manualBuild(contPoints);
this->local_points = unmatchedPointList.nbElem();
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->assoc.assign(nb_atoms,UINT_MAX);
if (unmatchedMeshList.getContainerElems().nbElem() == 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().nbElem());
/* 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
LM_ASSERT(i*Dim < this->positions.size(),"overflow detected");
Real x[3] = {0.0,0.0,0.0};
for (UInt k = 0; k < Dim; ++k) x[k] = this->positions[i*Dim+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->assoc[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->assoc[i] = j;
++assoc_found;
break;
}
}
if (it == subset_elts.end()){
this->assoc[i]=UINT_MAX;
}
}
STOPTIMER("Construction association");
/* build number of atoms per element */
for (UInt i = 0 ; i < nb_atoms ; ++i){
if (this->assoc[i] != UINT_MAX)
++nb_atoms_per_element[this->assoc[i]];
}
/* filter element list and node list for the unassociated node and elements */
filterContainerElems(nb_atoms_per_element);
buildNodeList();
IteratorElemsSubset itc = meshList.getContainerElems().getIterator();
RefElem el = itc.getFirst();
smatrix = new ShapeMatrix<Dim>(meshList.getContainerElems().nbElem());
el = itc.getFirst();
UInt i = 0;
while (!itc.end()){
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");
// if (nb_atoms_per_element[i]>0)
smatrix->setSub(i,nb_atoms_per_element[i],el.nNodes());
++i;
el = itc.getNext();
}
//Allocate the shapematrix structure (contiguous)
smatrix->allocate();
/* set the indirection (alloc tab) plus the shapes */
std::vector<Real> shapes;
std::vector<UInt > nodes;
UInt atom_index = 0;
for (i = 0 ; i < nb_atoms ; ++i){
if (this->assoc[i] == UINT_MAX)
continue;
Real X[3] = {0.0,0.0,0.0};
for (UInt k = 0; k < Dim; ++k)
X[k] = this->positions[i*Dim+k];
el = meshList.getContainerElems().get(this->assoc[i]);
el.computeShapes(shapes,X);
el.globalIndexes(nodes);
IF_TRACED(X,"checking for atom in trouble");
for(UInt j=0;j<shapes.size();++j)
smatrix->fill(atom_index,nodes[j],
meshList.subIndex2Index(nodes[j]),
this->assoc[i],shapes[j]);
++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.nbElem();
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,false);
}
/* -------------------------------------------------------------------------- */
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;
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.nbElem(),DBG_ALL);
new_index.push_back(new_t.nbElem());
DUMP("maintenant nb_atome_par_element[" << new_t.nbElem() << "]= "
<< nb_atome_par_element[i] << " (i=" << i << ")",DBG_ALL);
nb_atome_par_element[new_t.nbElem()] = nb_atome_par_element[i];
new_t.add(unmatchedMeshList.getContainerElems().get(i));
}
else
new_index.push_back(new_t.nbElem());
}
nb_atome_par_element.resize(new_t.nbElem());
meshList.empty();
//rebuild the element container
for (UInt i = 0 ; i < new_t.nbElem();++i) meshList.getContainerElems().add(new_t.get(i));
//rebuild the association list
for (UInt i = 0 ; i < this->assoc.size();++i){
if (this->assoc[i] != UINT_MAX){
this->assoc[i] = new_index[this->assoc[i]];
}
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints,ContainerMesh,Dim>::
filterPointListForUnmatched(){
pointList.empty();
for (UInt i = 0 ; i < this->assoc.size();++i)
if(this->assoc[i] != UINT_MAX){
DUMP("for container, atom " << i << " becomes " << pointList.nbElem(),DBG_ALL);
pointList.add(unmatchedPointList.get(i));
}
else{
DUMP("atom filtered " << i << " " << unmatchedPointList.get(i) ,DBG_DETAIL);
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints,ContainerMesh,Dim>::
filterArray(std::vector<Real> & array,UInt stride){
UInt index=0;
for (UInt i = 0 ; i < this->assoc.size();++i)
if(this->assoc[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(std::vector<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> & v_out,std::vector<Real> & data_atom,
UInt stride,UInt nmax){
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerPoints, typename ContainerMesh, UInt Dim>
void Bridging<ContainerPoints,ContainerMesh,Dim>::
solveLeastSquare(std::vector<Real> & mesh_data,std::vector<Real> & atomic_data,UInt stride){
// build right hand side of leastsquare system
mesh_data.assign(stride*meshList.nbElem(),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.nbElem(); ++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(std::vector<Real> & data_atomic,
std::vector<Real> & data_mesh,
UInt stride){
data_mesh.assign(stride*meshList.getContainerElems().nbElem(),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().empty();
//nodeList.empty();
assoc_found = 0;
this->assoc.clear();
pointList.empty();
}
/* -------------------------------------------------------------------------- */
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,
std::vector<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,
std::vector<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,
std::vector<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->assoc.size());
{
UInt indirection = 0;
for (UInt at = 0 ; at < this->assoc.size(); ++at){
current_indexes[at] = -1;
if (this->assoc[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->assoc[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->assoc[at_index];
this->smatrix->removeAtom(current_indexes[at_index],el);
this->assoc[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->assoc.size() ; ++i)
if (this->assoc[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->assoc.size(); ++at){
if (this->assoc[at] == UINT_MAX) continue;
LM_ASSERT(current_indexes[at] != -1,"This should not happend");
this->smatrix->changeAtomIndirection(this->assoc[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->assoc[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.nbElem();
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__
c4science · Help