Page MenuHomec4science
Files F65176772

bridging.cc
No OneTemporary
Actions

File Metadata

Created
Sat, Jun 1, 13:14

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.
*
*/
/* -------------------------------------------------------------------------- */
#include "bridging.hh"
#include "compute_extract.hh"
#include "factory_multiscale.hh"
#include "filter_geometry.hh"
#include "geometry_manager.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),
pointList(createOutput("pointList")), meshList(createOutput("meshList")),
contMesh(createOutput("contMesh")), grid(createOutput("grid")) {
assoc_found = 0;
unmatchedPointList.setCommGroup(this->comm_group_A);
unmatchedMeshList.setCommGroup(this->comm_group_B);
// contains the elements selected by geometry
this->createOutput("unmatchedMeshList") = unmatchedMeshList;
// contains the points selected by geometry
this->createOutput("unmatchedPointList") = unmatchedPointList;
is_in_continuum = this->comm_group_B.amIinGroup();
is_in_atomic = this->comm_group_A.amIinGroup();
}
/* -------------------------------------------------------------------------- */
Bridging::~Bridging() {}
/* -------------------------------------------------------------------------- */
void Bridging::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;
}
}
/* -------------------------------------------------------------------------- */
// void Bridging::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);
// }
// }
/* -------------------------------------------------------------------------- */
void Bridging::mesh2Point(const ContainerArray<Real> &data_node,
ContainerArray<Real> &data_atom) {
if (this->local_points == 0) {
DUMP("Warning : atomic container is empty:"
<< " cannot interpolate mesh fields on control points (check "
"geometries ?!)",
DBG_WARNING);
return;
}
data_atom = smatrix.transpose() * data_node.matrix();
}
/* -------------------------------------------------------------------------- */
void Bridging::mesh2Point(FieldType field_type,
ContainerArray<Real> &data_point) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field_type);
extracted_field.compute(this->meshList);
this->mesh2Point(extracted_field.evalArrayOutput(), data_point);
}
/* -------------------------------------------------------------------------- */
void Bridging::leastSquarePointData(ContainerArray<Real> &,
ContainerArray<Real> &, UInt) {
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
void Bridging::copySlaveValues(ContainerArray<Real> &v) {
for (auto &&pair : local_pbc_pairs) {
UInt i1 = pair.first;
UInt i2 = pair.second;
v(i1) = v(i2);
}
}
/* -------------------------------------------------------------------------- */
/* Parallel Methods */
/* -------------------------------------------------------------------------- */
void Bridging::commPoint2Continuum(ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->distributeVectorA2B("communicate buffer", buf);
}
}
/* -------------------------------------------------------------------------- */
void Bridging::commContinuum2Point(ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->distributeVectorB2A("communicate buffer", buf);
}
}
/* -------------------------------------------------------------------------- */
void Bridging::synchSumBuffer(ContainerArray<Real> &buf) {
if (this->local_points != 0) {
this->synchronizeVectorBySum("synchsumvector", buf);
}
}
/* -------------------------------------------------------------------------- */
void Bridging::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);
LM_TOIMPLEMENT;
// 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");
LM_TOIMPLEMENT;
// this->smatrix->changeAtomIndirection(this->pointToElement[at],
// current_indexes[at],
// indirection);
++indirection;
}
LM_TOIMPLEMENT;
// 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];
// }
}
}
/* -------------------------------------------------------------------------- */
void Bridging::setPBCPairs(std::vector<std::pair<UInt, UInt>> &pairs) {
if (this->is_in_continuum)
pbc_pairs = pairs;
}
/* -------------------------------------------------------------------------- */
UInt Bridging::getNumberPoints() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
UInt Bridging::getNumberElems() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
UInt Bridging::getNumberNodes() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
void Bridging::attachVector(ContainerArray<Real> &tab) {
this->pointList.get<ContainerInterface>().attachObject(tab);
}
/* -------------------------------------------------------------------------- */
void Bridging::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->getOutput<ContainerInterface>("pointList").size();
this->positions.resize(this->local_points, spatial_dimension);
}
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 */
void Bridging::declareParams() {
this->addSubParsableObject(unmatchedMeshList);
DofAssociation::declareParams();
}
/* -------------------------------------------------------------------------- */
void Bridging::setPointField(FieldType field, ContainerArray<Real> &buffer) {
StimulationField impose_field("impose_point_field");
impose_field.setParam("FIELD", field);
impose_field.compute("dofs"_input = this->pointList, "field"_input = buffer);
}
/* -------------------------------------------------------------------------- */
void Bridging::setPointField(FieldType field, Real val, UInt dim) {
ContainerArray<Real> data_point(pointList.get<ContainerInterface>().size(),
dim);
for (auto &&v : data_point.rowwise()) {
v = val;
}
setPointField(field, data_point);
}
/* -------------------------------------------------------------------------- */
void Bridging::addPointField(FieldType field, ContainerArray<Real> &buffer) {
StimulationField impose_field("impose_point_field");
impose_field.setParam("FIELD", field);
impose_field.setParam("ADDITIVE", true);
impose_field.compute("dofs"_input = this->pointList, "field"_input = buffer);
}
/* -------------------------------------------------------------------------- */
void Bridging::setMeshField(FieldType field, ContainerArray<Real> &buffer) {
StimulationField impose_field("impose_mesh_field");
impose_field.setParam("FIELD", field);
impose_field.compute("dofs"_input = this->meshList, "field"_input = buffer);
}
/* -------------------------------------------------------------------------- */
void Bridging::addMeshField(FieldType field, ContainerArray<Real> &buffer) {
StimulationField impose_field("impose_mesh_field");
impose_field.setParam("FIELD", field);
impose_field.setParam("ADDITIVE", true);
impose_field.compute("dofs"_input = this->meshList, "field"_input = buffer);
}
/* -------------------------------------------------------------------------- */
void Bridging::setPointField(FieldType field, Eigen::VectorXd val) {
ContainerArray<Real> data_point(pointList.get<ContainerInterface>().size(),
val.rows());
for (auto &&v : data_point.rowwise()) {
v = val;
}
setPointField(field, data_point);
}
/* -------------------------------------------------------------------------- */
void Bridging::setMeshField(FieldType field, Eigen::VectorXd val) {
ContainerArray<Real> data_mesh(pointList.get<ContainerInterface>().size(),
val.rows());
data_mesh = val;
setPointField(field, data_mesh);
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__

Event Timeline

c4science · Help