diff --git a/src/coupling/arlequin_template.cc b/src/coupling/arlequin_template.cc
index 4239ae6..df8144a 100644
--- a/src/coupling/arlequin_template.cc
+++ b/src/coupling/arlequin_template.cc
@@ -1,182 +1,158 @@
/**
* @file arlequin_template.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 15:05:56 2013
*
* @brief Internal class to factor code for the Arlequin kind methods
*
* @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 "arlequin_template.hh"
#include "bridging.hh"
#include "bridging_atomic_continuum.hh"
#include "factory_multiscale.hh"
#include "lm_common.hh"
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
ArlequinTemplate::ArlequinTemplate(const std::string &name)
: LMObject(name), CouplingAtomicContinuum(name),
bridging_zone(this->getID() + "-bridging"),
boundary_zone(this->getID() + "-boundary"),
weight_mesh("weight-mesh-" + name), weight_point("weight-point-" + name),
lambdas_mesh("lambdas-mesh-" + name),
lambdas_point("lambdas-point-" + name) {
this->weight_mesh.setCommGroup(comm_group_continuum());
this->lambdas_mesh.setCommGroup(comm_group_continuum());
this->weight_point.setCommGroup(comm_group_atomic());
this->lambdas_point.setCommGroup(comm_group_atomic());
this->createOutput("weight-mesh") = this->weight_mesh;
this->createOutput("lambdas-mesh") = this->lambdas_mesh;
this->createOutput("weight-point") = this->weight_point;
this->createOutput("lambdas-point") = this->lambdas_point;
}
/* -------------------------------------------------------------------------- */
-INSTANCIATE_DISPATCH(ArlequinTemplate::computeWeights)
-template <typename ContA, typename ContC>
-void ArlequinTemplate::computeWeights(ContA &pointList, ContC &meshList) {
+INSTANCIATE_DISPATCH(ArlequinTemplate::computeAtomicWeights)
+template <typename ContA>
+void ArlequinTemplate::computeAtomicWeights(ContA &pointList) {
// point weights
- if (this->is_in_atomic()) {
- this->weight_point.compute(pointList);
- bridging_zone.attachVector(this->weight_point.evalArrayOutput());
+ this->weight_point.compute(pointList);
+ bridging_zone.attachVector(this->weight_point.evalArrayOutput());
- lambdas_point.assign(weight_point.evalArrayOutput().size(), 0);
- this->bridging_zone.attachVector(lambdas_point);
+ lambdas_point.assign(weight_point.evalArrayOutput().size(), 0);
+ this->bridging_zone.attachVector(lambdas_point);
- auto &&weights = weight_point.evalArrayOutput();
+ auto &&weights = weight_point.evalArrayOutput();
- for (auto &&[at, lbda, weight] : zip(pointList, lambdas_point, weights)) {
- lbda = weight * at.mass();
- }
+ for (auto &&[at, lbda, weight] : zip(pointList, lambdas_point, weights)) {
+ lbda = weight * at.mass();
}
+}
+/* -------------------------------------------------------------------------- */
+
+INSTANCIATE_DISPATCH(ArlequinTemplate::computeContinuumWeights)
+template <typename ContC>
+void ArlequinTemplate::computeContinuumWeights(ContC &meshList) {
// mesh weights
- if (this->is_in_continuum()) {
- this->weight_mesh.compute(meshList);
- lambdas_mesh.assign(weight_mesh.evalArrayOutput().size(), 0);
+ this->weight_mesh.compute(meshList);
+ lambdas_mesh.assign(weight_mesh.evalArrayOutput().size(), 0);
- for (auto &&[nd, lbda, weight] :
- zip(meshList.getContainerNodes(), lambdas_mesh,
- weight_mesh.evalArrayOutput())) {
- lbda = weight * nd.mass();
- }
+ for (auto &&[nd, lbda, weight] :
+ zip(meshList.getContainerNodes(), lambdas_mesh,
+ weight_mesh.evalArrayOutput())) {
+ lbda = weight * nd.mass();
}
}
/* -------------------------------------------------------------------------- */
-// void ArlequinTemplate::allocate(UInt t) {
-// size_constraint = t;
-// const UInt Dim = spatial_dimension;
-// DUMP("initial number of atoms in rec " << size_constraint,
-// DBG_INFO_STARTUP); DUMP("allocation of " << Dim * (size_constraint) << "
-// Reals",
-// DBG_INFO_STARTUP);
-// A.assign(size_constraint, 0);
-// rhs.resize(size_constraint, Dim);
-// rhs.array().setZero();
-
-// DUMP("Attach vector A", DBG_INFO_STARTUP);
-// bridging_zone.attachVector(A);
-//}
-
-/* -------------------------------------------------------------------------- */
-
-// void ArlequinTemplate::cleanRHS() {
-// rhs.resize(this->bridging_zone.getNumberLocalMatchedPoints(),
-// spatial_dimension);
-// rhs.array().setZero();
-// }
-
-// /* --------------------------------------------------------------------------
-// */
-
/* LMDESC ArlequinTemplate
This class is used internally.
It is to be used
while two zones are declared, one for the coupling
and one for providing a stiff boundary condition to the atoms.
In the coupling a linear weight function is built.
*/
/* LMHERITANCE dof_association compute_arlequin_weight */
void ArlequinTemplate::declareParams() {
this->addSubParsableObject(bridging_zone);
this->addSubParsableObject(boundary_zone);
this->addSubParsableObject(weight_mesh);
this->addSubParsableObject(weight_point);
/* LMKEYWORD GEOMETRY
Set the bridging/overlaping zone where the Lagrange multipliers are
to be computed.
*/
this->parseKeyword("GEOMETRY", bridging_geom);
/* LMKEYWORD BOUNDARY
Set the boundary geometry where the atom velocities are to be fixed from
the interpolated finite element velocities.
*/
this->parseKeyword("BOUNDARY", boundary_geom);
/* LMKEYWORD CHECK_COHERENCY
Perform a systematic check of the communication scheme.
**Careful, it is computationally expensive**
*/
this->parseTag("CHECK_COHERENCY", check_coherency, false);
}
/* ------------------------------------------------------------------------ */
__END_LIBMULTISCALE__
diff --git a/src/coupling/arlequin_template.hh b/src/coupling/arlequin_template.hh
index 6826a6d..3b28026 100644
--- a/src/coupling/arlequin_template.hh
+++ b/src/coupling/arlequin_template.hh
@@ -1,103 +1,105 @@
/**
* @file arlequin_template.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 15:05:56 2013
*
* @brief Internal class to factor code for the Arlequin kind methods
*
* @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.
*
*/
#ifndef __LIBMULTISCALE_ARLEQUIN_TEMPLATE_HH__
#define __LIBMULTISCALE_ARLEQUIN_TEMPLATE_HH__
/* -------------------------------------------------------------------------- */
#include "bridging_atomic_continuum.hh"
#include "compute_arlequin_weight.hh"
#include "coupling_atomic_continuum.hh"
#include "lib_continuum.hh"
#include "lib_md.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class ArlequinTemplate : public CouplingAtomicContinuum {
public:
ArlequinTemplate(const std::string &name);
~ArlequinTemplate() = default;
void declareParams();
//! The bridging zone
BridgingAtomicContinuum bridging_zone;
//! The bridge zone used for surface effects
BridgingAtomicContinuum boundary_zone;
protected:
- DECORATE_FUNCTION_DISPATCH(computeWeights, subMD, subCONTINUUM)
+ DECORATE_FUNCTION_DISPATCH(computeAtomicWeights, subMD)
//! Compute the weighting function of the coupled MD DOFs
- template <typename ContA, typename ContC>
- void computeWeights(ContA &pointList, ContC &meshList);
+ template <typename ContA> void computeAtomicWeights(ContA &pointList);
+ DECORATE_FUNCTION_DISPATCH(computeContinuumWeights, subCONTINUUM)
+ //! Compute the weighting function of the coupled MD DOFs
+ template <typename ContC> void computeContinuumWeights(ContC &meshList);
protected:
//! geometry id for the bridging zone
LMID bridging_geom;
//! geometry id for the boundary zone
LMID boundary_geom;
//! weighting vector for FE DOFs
ComputeArlequinWeight weight_mesh;
//! weighting vector for MD DOFs
ComputeArlequinWeight weight_point;
//! weighting vector for FE DOFs
ContainerArray<Real> lambdas_mesh;
//! weighting vector for MD DOFs
ContainerArray<Real> lambdas_point;
//! flag to know wether we have to chack for coherency
bool check_coherency;
};
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
#endif /* __LIBMULTISCALE_ARLEQUIN_TEMPLATE_HH__ */
diff --git a/src/coupling/bridging.cc b/src/coupling/bridging.cc
index 63da7f0..f9f6622 100644
--- a/src/coupling/bridging.cc
+++ b/src/coupling/bridging.cc
@@ -1,417 +1,416 @@
/**
* @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("unmatchedPointList-" + name),
unmatchedMeshList("unmatchedMeshList-" + 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;
}
/* -------------------------------------------------------------------------- */
Bridging::~Bridging() {}
/* -------------------------------------------------------------------------- */
void Bridging::filterArray(ContainerArray<Real> &array) {
UInt index = 0;
for (UInt i = 0; i < pointToElement.size(); ++i)
if (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 (not this->in_group_B())
return;
DUMP("local points", DBG_MESSAGE);
if (this->local_points == 0)
return;
// loop through the unassociated atoms declared per each processor
// The total list of atom index rejected is constructed
UInt offset = 0;
UInt cpt_unassociated = 0;
std::vector<UInt> atom_to_remove;
for (auto &&[proc, c_w] : enumerate(com_with)) {
if (not c_w)
continue;
DUMP("atomic processor " << proc << " unassociated "
<< unassociated_atoms[proc].size() << " atoms ",
DBG_INFO_STARTUP);
for (auto &&unassociated : unassociated_atoms[proc]) {
// reconstruct the index for atom declared unassociated
UInt unassociated_index = offset + unassociated;
atom_to_remove.push_back(unassociated_index);
pointToElement[unassociated_index] = UINT_MAX;
// auto el_index = this->pointToElement[unassociated_index];
// LM_ASSERT(el_index != UINT_MAX,
// "this should not happend !! " << unassociated_index);
// // remove the column associated to atom in the shapematrix
// this->smatrix.removeAtom(unassociated_index);
// // alter the atom-element mapping
// pointToElement[unassociated_index] = UINT_MAX;
++cpt_unassociated;
DUMP("unassociated atom at position "
<< unassociated_index << " whereas offset = " << offset
<< " next offset is " << offset + this->nb_points_per_proc[proc],
DBG_INFO_STARTUP);
}
offset += this->nb_points_per_proc[proc];
}
smatrix.removeAtoms(atom_to_remove);
// corrects the local_points variable
if (cpt_unassociated) {
DUMP("removed " << cpt_unassociated
<< " atoms from unassociation. Now local atoms is "
<< this->local_points,
DBG_MESSAGE);
this->local_points -= cpt_unassociated;
DUMP("and becomes " << this->local_points, DBG_MESSAGE);
}
// compute contiguous indexes for atoms (compression of indexes)
// the result is a mapping from old indexes to new ones
// output: mapping_old_to_new_indexes
// std::vector<UInt> mapping_old_to_new_indexes;
// mapping_old_to_new_indexes.assign(pointToElement.size(), UINT_MAX);
// UInt indirection = 0;
// for (auto &&[at_index, el_index] : enumerate(pointToElement)) {
// if (el_index == UINT_MAX)
// continue;
// mapping_old_to_new_indexes[at_index] = indirection;
// ++indirection;
// }
// LM_ASSERT(indirection == this->local_points, "this should not happend "
// << indirection << " "
// << this->local_points);
// checks that number of mapped atoms equals local atoms variable
UInt cpt_tmp = 0;
for (auto &&el_index : pointToElement) {
if (el_index != 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
// indirection = 0;
// for (auto &&[at_index, el_index] : enumerate(pointToElement)) {
// if (el_index == UINT_MAX)
// continue;
// UInt new_index = mapping_old_to_new_indexes[at_index];
// LM_ASSERT(new_index != UINT_MAX, "This should not happend");
// smatrix.changeAtomIndirection(at_index, new_index, indirection);
// ++indirection;
// }
// this->smatrix.swapAtomIndirections();
// now i set the duo vector
this->createDuoVectorB("FE", pointToElement);
}
/* -------------------------------------------------------------------------- */
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:
- unmatched-fe-\${couplerID} : the set of nodes and elements
that are contained in the geometry provided by GEOMETRY keyword.
- unmatched-point-\${couplerID} : the set of points
that are contained in the geometry provided by GEOMETRY keyword.
- matched-fe-\${couplerID} : the set of nodes and elements
that are containing at least one point of unmatched-point-\${couplerID}.
- matched-point-\${couplerID} : the set of points
that are contained in at least one element of unmatched-fe-\${couplerID}.
*/
/* LMHERITANCE filter_geometry dof_association */
void Bridging::declareParams() {
this->addSubParsableObject(unmatchedMeshList);
DofAssociation::declareParams();
}
/* -------------------------------------------------------------------------- */
void Bridging::setPointField(FieldType field, ContainerArray<Real> &buffer) {
+ if (not this->is_in_atomic())
+ return;
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);
+ ContainerArray<Real> data_point(this->local_points, 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());
+ ContainerArray<Real> data_point(this->local_points, 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());
+ ContainerArray<Real> data_mesh(this->local_points, val.rows());
data_mesh = val;
setPointField(field, data_mesh);
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
diff --git a/src/coupling/bridging_atomic_continuum.cc b/src/coupling/bridging_atomic_continuum.cc
index 877d458..bc0d408 100644
--- a/src/coupling/bridging_atomic_continuum.cc
+++ b/src/coupling/bridging_atomic_continuum.cc
@@ -1,344 +1,341 @@
/**
* @file bridging_atomic_continuum.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief Bridging object between atomistic and finite elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
//#define TIMER
/* -------------------------------------------------------------------------- */
#include "bridging_atomic_continuum.hh"
#include "compute_extract.hh"
#include "lm_common.hh"
#include "reference_manager.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
BridgingAtomicContinuum::BridgingAtomicContinuum(const std::string &name)
: LMObject(name), Bridging(name), comm_group_continuum(this->comm_group_B),
comm_group_atomic(this->comm_group_A) {}
/* -------------------------------------------------------------------------- */
BridgingAtomicContinuum::~BridgingAtomicContinuum() {}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::updateForMigration() {
Bridging::updateForMigration();
if (this->check_coherency) {
checkPositionCoherency();
}
}
/* -------------------------------------------------------------------------- */
INSTANCIATE_DISPATCH(BridgingAtomicContinuum::init)
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum::init(DomainA &domA, DomainC &domC) {
auto all_group = Communicator::getCommunicator().getGroup("all");
Bridging::init(domA.getContainer(), domC.getContainer());
- auto &pointList = this->pointList.get<typename DomainA::ContainerSubset>();
-
- domA.getContainer().addSubSet(pointList);
-
- DUMP(this->getID() << " : attach duo vector", DBG_MESSAGE);
-
if (in_group_A()) {
+ auto &pointList = this->pointList.get<typename DomainA::ContainerSubset>();
+ domA.getContainer().addSubSet(pointList);
+ DUMP(this->getID() << " : attach duo vector", DBG_MESSAGE);
pointList.attachObject(this->getDuoVector());
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
}
all_group.synchronize();
}
/* -------------------------------------------------------------------------- */
// void BridgingAtomicContinuum::averageDOFsOnElements(
// ContainerArray<Real> &data_mesh, FieldType field) {
// LM_FATAL("to adapt");
// if (comm_group_atomic == comm_group_continuum) {
// ComputeExtract extracted_field("ComputeExtract:" + this->getID());
// extracted_field.setParam("FIELD", field);
// extracted_field.buildManual(*this->pointList);
// this->averageOnElements<dispatch>(extracted_field.getArray(), data_mesh,
// meshList);
// return;
// }
// ComputeExtract extracted_field("ComputeExtract:" + this->getID());
// extracted_field.setParam("FIELD", field);
// if (this->is_in_atomic)
// extracted_field.buildManual(*this->pointList);
// LM_TOIMPLEMENT;
// // extracted_field.resize(Parent::local_points * Dim);
// /* do communication */
// this->communicateBufferAtom2Continuum(extracted_field.getArray());
// if (this->is_in_continuum)
// this->averageOnElements<dispatch>(extracted_field.getArray(), data_mesh,
// meshList);
// }
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::leastSquareAtomicData(
ContainerArray<Real> &data_mesh, ContainerArray<Real> &data_atom) {
if (comm_group_atomic == comm_group_continuum) {
this->solveLeastSquare<dispatch>(data_mesh, data_atom,
this->evalOutput("meshList"));
return;
}
/* do communication */
this->commPoint2Continuum(data_atom);
if (this->is_in_continuum()) {
this->solveLeastSquare<dispatch>(data_mesh, data_atom,
this->getOutput("meshList"));
}
}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::leastSquareAtomicDOFs(
ContainerArray<Real> &data_mesh, FieldType field) {
LM_FATAL("to adapt");
if (comm_group_atomic == comm_group_continuum) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
extracted_field.compute(this->pointList);
this->solveLeastSquare<dispatch>(
data_mesh, extracted_field.evalArrayOutput(), this->meshList);
return;
}
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
if (this->is_in_atomic())
extracted_field.compute(this->pointList);
LM_TOIMPLEMENT;
// extracted_field.resize(Parent::local_points * Dim);
this->leastSquareAtomicData(data_mesh, extracted_field.evalArrayOutput());
}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::buildAtomicDOFsNoiseVector(std::vector<Real> &,
FieldType field) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
if (this->is_in_atomic())
extracted_field.compute(this->pointList);
if (this->is_in_continuum()) {
LM_TOIMPLEMENT;
// auto &field = extracted_field.evalArrayOutput();
// this->interpolatePointData(field, ??);
for (UInt i = 0; i < this->local_points * spatial_dimension; ++i)
LM_TOIMPLEMENT;
// extracted_field[i] *= -1;
}
/* do communication */
this->synchSumBuffer(extracted_field.evalArrayOutput());
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum::checkPositionCoherency(DomainA &domA,
DomainC &domC) {
LM_TOIMPLEMENT;
// constexpr UInt Dim = DomainC::ContainerPoints::Dim;
// if (this->in_group_B) {
// // interpolate atomic sites to check everything is ok from this
// auto &&field = domC.getField(_position0);
// UInt *assoc_tmp = &this->pointToElement[0];
// for (UInt i = 0; i < this->local_points; ++i) {
// Vector<Dim> X;
// Vector<Dim> pos = this->positions(i);
// X = this->interpolate(i, field, assoc_tmp[i]);
// for (UInt k = 0; k < Dim; ++k) {
// if (X[k] != pos[k])
// if (fabs(pos[k] - X[k]) > 1e-5) {
// searchAtomInLocalPool(
// static_cast<typename DomainA::ContainerSubset &>(
// this->pointList),
// pos);
// DUMP("mismatch in positions (X), atom "
// << i << " dist is " << pos[k] - X[k]
// << " : probably a redistribution trouble check "
// "migrations\n"
// << "computed position is " << X[0] << " " << X[1] << " "
// << X[2] << " position was " << pos[0] << " " << pos[1]
// << " " << pos[2],
// DBG_MESSAGE);
// }
// }
// }
// }
// this->distributeVectorB2A("temp_positions", this->positions);
// if (this->in_group_A) {
// UInt i = 0;
// UInt global_flag = 0;
// for (auto &&at : this->pointList) {
// UInt local_flag = 0;
// Vector<Dim> pos = VectorView<Dim>(&this->positions[Dim * i]);
// Vector<Dim> at_pos = at.position0();
// for (UInt k = 0; k < Dim; ++k) {
// if (pos[k] != at_pos[k])
// if (!local_flag && fabs(pos[k] - at_pos[k]) > 1e-5) {
// searchAtomInLocalPool(this->pointList, pos);
// DUMP("mismatch in positions (X), atom "
// << i << " dist is "
// << this->positions[Dim * i] - at.position0()[0]
// << " : probably a redistribution trouble check
// migrations"
// << std::endl
// << "ref is " << at << "\nand send position is " <<
// pos[0]
// << " " << pos[1] << " " << pos[2],
// DBG_MESSAGE);
// local_flag = 1;
// global_flag = 1;
// }
// }
// ++i;
// }
// if (global_flag) {
// this->getDuoVector().print("buggy-redistrib");
// LM_FATAL("abort due to migration problem");
// }
// }
// DUMP("position coherency OK", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename Cont, typename Vec>
void BridgingAtomicContinuum::searchAtomInLocalPool(Cont &container, Vec &x) {
LM_TOIMPLEMENT;
// UInt i = 0;
// for (auto &&at : container) {
// #ifdef TRACE_ATOM
// auto pos = at.position0();
// #endif
// IF_COORD_EQUALS(pos, x) {
// DUMP("actually searched atom is at ref " << at << " bridge index " <<
// i,
// DBG_MESSAGE);
// break;
// }
// ++i;
// }
// if (i == container.size())
// DUMP("atom " << x[0] << " " << x[1] << " " << x[2] << " "
// << " not found locally : probably has migrated"
// << " or is simply lost !"
// << " NOT GOOD",
// DBG_MESSAGE);
// i = 0;
// for (auto &&at : domA.getContainer()) {
// #ifdef TRACE_ATOM
// auto pos = at.position0();
// #endif
// IF_COORD_EQUALS(pos, x) {
// DUMP("actually searched atom is at ref " << at << " bridge index " <<
// i,
// DBG_MESSAGE);
// break;
// }
// ++i;
// }
// if (i == domA.getContainer().size())
// DUMP("atom not found locally : probably has migrated"
// << "or is simply lost !"
// << " NOT GOOD",
// DBG_MESSAGE);
}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::projectAtomicFieldOnMesh(FieldType field_type) {
this->projectAtomicFieldOnMesh(field_type, this->buffer_for_points);
if (is_in_atomic())
this->setPointField(field_type, this->buffer_for_points);
}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::projectAtomicFieldOnMesh(
FieldType field_type, ContainerArray<Real> &buffer) {
const UInt Dim = spatial_dimension;
if (is_in_continuum() && !this->total_points)
DUMP("Warning : atomic container is empty: "
<< "cannot proceed atomic projection (check geometries ?!) "
<< this->getID(),
DBG_WARNING);
if (!this->local_points)
return;
buffer.resize(this->local_points, Dim);
if (is_in_continuum()) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field_type);
extracted_field.compute(this->meshList);
this->mesh2Point(extracted_field.evalArrayOutput(), buffer);
buffer.setCommGroup(comm_group_continuum);
}
this->distributeVectorB2A("correction_fictifs", buffer);
}
/* -------------------------------------------------------------------------- */
void BridgingAtomicContinuum::checkPositionCoherency() { LM_TOIMPLEMENT; }
__END_LIBMULTISCALE__
diff --git a/src/coupling/bridging_inline_impl.hh b/src/coupling/bridging_inline_impl.hh
index 8303829..2a5837e 100644
--- a/src/coupling/bridging_inline_impl.hh
+++ b/src/coupling/bridging_inline_impl.hh
@@ -1,502 +1,507 @@
/**
* @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) {
auto all_group = Communicator::getCommunicator().getGroup("all");
DofAssociation::init(contA, contC);
const UInt Dim = spatial_dimension;
contMesh = contC;
- pointList.alloc<typename ContA::ContainerSubset>("pointList-" +
- this->getID());
- meshList.alloc<typename ContC::ContainerSubset>("meshList-" + this->getID());
-
- pointList->acquireContext(contA);
- meshList->acquireContext(contC);
- unmatchedPointList.acquireContext(contA);
- unmatchedMeshList.acquireContext(contC);
+ if (in_group_A()) {
+ pointList.alloc<typename ContA::ContainerSubset>("pointList-" +
+ this->getID());
+ pointList->acquireContext(contA);
+ unmatchedPointList.acquireContext(contA);
+ }
+
+ if (in_group_B()) {
+ meshList.alloc<typename ContC::ContainerSubset>("meshList-" +
+ this->getID());
+ meshList->acquireContext(contC);
+ unmatchedMeshList.acquireContext(contC);
+ }
buffer_for_points.acquireContext(contA);
buffer_for_nodes.acquireContext(contC);
DUMPBYPROC("selecting DOFs in bridging zone", DBG_INFO_STARTUP, 0);
if (in_group_A()) {
this->buildPointList(contA);
this->buildPositions<dispatch>(unmatchedPointList.evalOutput());
}
if (in_group_B())
this->buildContinuumDOFsList(contC);
DUMP(this->getID() << " : Generating Communication Scheme", DBG_INFO_STARTUP);
all_group.synchronize();
DUMP(this->getID() << " : exchange coarse geometries", DBG_MESSAGE);
if (this->in_group_A()) {
this->exchangeGeometries<dispatch>(this->unmatchedPointList.evalOutput());
}
if (in_group_B()) {
this->exchangeGeometries<dispatch>(this->unmatchedMeshList.evalOutput());
}
all_group.synchronize();
DUMP(this->getID() << " : exchange points positions", DBG_MESSAGE);
this->exchangePositions(Dim);
all_group.synchronize();
DUMP(this->getID() << " : build shape matrix", DBG_MESSAGE);
if (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);
}
all_group.synchronize();
DUMP(this->getID() << " : ExchangeAssociation", DBG_MESSAGE);
ContainerArray<UInt> managed;
this->exchangeAssociationInformation(managed, pointToElement);
for (auto &&[i, c_w] : enumerate(this->com_with))
DUMP("second test Do I communicate with " << i << " : " << c_w, DBG_INFO);
all_group.synchronize();
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, pointToElement, unassociated);
DUMP(this->getID() << " : doing exchange rejected", DBG_MESSAGE);
this->exchangeRejectedContinuumOwners(unassociated);
all_group.synchronize();
DUMP(this->getID() << " : doing filter rejected", DBG_MESSAGE);
if (this->in_group_B())
this->filterRejectedContinuumOwners(unassociated);
all_group.synchronize();
DUMP(this->getID() << " : filtering position & association vector",
DBG_MESSAGE);
if (this->in_group_A()) {
filterPointListForUnmatched<dispatch>(pointList);
}
if (this->in_group_B()) {
this->filterArray(this->positions);
this->filterMapping(pointToElement);
// this->buildNodeShape<dispatch>(this->meshList);
}
DUMP(this->getID() << " : all done for scheme generation", DBG_MESSAGE);
all_group.synchronize();
}
/* -------------------------------------------------------------------------- */
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 per 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) {
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();
this->removeOutput("grid");
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);
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(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);
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 (pointToElement[i] != UINT_MAX)
++nb_points_per_element[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, 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 (pointToElement[i] == UINT_MAX)
continue;
auto &&X = this->positions(i);
auto el = elements.get(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 ContMesh>
void Bridging::filterContainerElems(std::vector<UInt> &nb_atom_per_element,
ContMesh &contMesh) {
using ContainerSubset = typename ContMesh::ContainerSubset;
typename ContainerSubset::ContainerElems new_t(this->getID() + ":elems");
std::vector<int> new_index;
auto &unmatchedMeshList =
this->unmatchedMeshList.evalOutput<ContainerSubset>();
auto &meshList = this->meshList.get<ContainerSubset>();
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 < pointToElement.size(); ++i) {
if (pointToElement[i] != UINT_MAX) {
pointToElement[i] = new_index[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 ContA>
void Bridging::filterPointListForUnmatched(ContA &pointList) {
auto &unmatchedPointList =
this->unmatchedPointList.evalOutput<typename ContA::ContainerSubset>();
pointList.clear();
for (UInt i = 0; i < pointToElement.size(); ++i)
if (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__
diff --git a/src/coupling/kobayashi.cc b/src/coupling/kobayashi.cc
index 29cc80b..489998d 100644
--- a/src/coupling/kobayashi.cc
+++ b/src/coupling/kobayashi.cc
@@ -1,189 +1,189 @@
/**
* @file kobayashi.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief Kobayashi's bridging method
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "kobayashi.hh"
#include "arlequin_template.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
Kobayashi::Kobayashi(const std::string &name)
: LMObject(name), ArlequinTemplate(name) {
// naming the bridging zones
std::stringstream fname;
fname << name << "-bridging";
}
/* -------------------------------------------------------------------------- */
Kobayashi::~Kobayashi() {}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::init(DomainA &domA, DomainC &domC) {
if (lm_world_size > 1)
LM_FATAL("this coupler works only in sequential");
std::stringstream fname;
// initialize the bridging_zone object
bridging_zone.init(domA, domC);
// build weights
- this->computeWeights<dispatch>(this->bridging_zone.pointList,
- this->bridging_zone.meshList);
+ this->computeAtomicWeights<dispatch>(this->bridging_zone.pointList);
+ this->computeContinuumWeights<dispatch>(this->bridging_zone.meshList);
DUMP("Computing initial least square constraint", DBG_INFO);
buildCGMatrix();
// treat now MD boundary zone
boundary_zone.init(domA, domC);
DUMP("Init ... done", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::coupling(CouplingStage stage, DomainA &domA [[gnu::unused]],
DomainC &domC [[gnu::unused]]) {
if (stage != COUPLING_STEP4)
return;
boundary_zone.projectAtomicFieldOnMesh(_velocity);
boundary_zone.projectAtomicFieldOnMesh(_displacement);
correctContinuum(_velocity);
correctContinuum(_displacement);
correctAtoms(_velocity);
correctAtoms(_displacement);
}
/* -------------------------------------------------------------------------- */
void Kobayashi::buildCGMatrix() {
auto &smatrix = bridging_zone.smatrix;
this->A = smatrix * smatrix.transpose();
this->solver.compute(this->A);
}
/* -------------------------------------------------------------------------- */
ContainerArray<Real> &Kobayashi::buildLeastSquareField(FieldType field_type) {
auto &smatrix = bridging_zone.smatrix;
ComputeExtract field(this->getID());
field.setParam("FIELD", field_type);
field.compute(this->bridging_zone.pointList);
Eigen::VectorXd rhs = smatrix * field.evalArrayOutput().matrix();
Eigen::VectorXd x = this->solver.solve(rhs.matrix());
auto &ls_field = this->bridging_zone.buffer_for_nodes;
ls_field.resize(x.rows(), x.cols());
ls_field = x;
return ls_field;
}
/* -------------------------------------------------------------------------- */
void Kobayashi::correctContinuum(FieldType field_type) {
auto &LS_field = this->buildLeastSquareField(field_type);
auto &weight_mesh = this->weight_mesh.evalArrayOutput();
ComputeExtract field(this->getID());
field.setParam("FIELD", field_type);
field.compute(this->bridging_zone.meshList);
for (auto &&[w, ls_f, f] : zip(weight_mesh, LS_field.rowwise(),
field.evalArrayOutput().rowwise())) {
if (w <= .5) {
Real speed_projection = this->quality;
if (w <= 1e-2)
speed_projection = 1;
ls_f = speed_projection * (ls_f - f);
}
}
this->bridging_zone.addMeshField(field_type, LS_field);
}
/* -------------------------------------------------------------------------- */
void Kobayashi::correctAtoms(FieldType field_type) {
ContainerArray<Real> LS_field;
bridging_zone.projectAtomicFieldOnMesh(field_type, LS_field);
auto &weightMD = this->weight_point.evalArrayOutput();
ComputeExtract field(this->getID());
field.setParam("FIELD", field_type);
field.compute(this->bridging_zone.pointList);
for (auto &&[w, ls_f, f] :
zip(weightMD, LS_field.rowwise(), field.evalArrayOutput().rowwise())) {
if (w <= .5) {
Real speed_projection = this->quality;
if (w <= 1e-2)
speed_projection = 1;
ls_f = speed_projection * (ls_f - f);
}
}
this->bridging_zone.addPointField(field_type, LS_field);
}
/* -------------------------------------------------------------------------- */
/* LMDESC KOBAYASHI
This keyword is used to set a coupling method
which will implement the method of Kobayashi and co-workers
*A simple dynamical scale-coupling method for concurrent simulation of
hybridized atomistic/coarse-grained-particle system,*
**Ryo Kobayashi, Takahide Nakamura, Shuji Ogata**,
*International Journal for Numerical Methods in Engineering Volume 83,
Issue 2, pages 249-268, 9 July 2010.*
This coupling component works only in 1D and in sequential.
*/
/* LMHERITANCE arlequin_template */
/* LMEXAMPLE
COUPLING_CODE md fe KOBAYASHI GEOMETRY 3 BOUNDARY 4 GRID_DIVISIONX 10 QUALITY
1e-3 VERBOSE
*/
void Kobayashi::declareParams() { ArlequinTemplate::declareParams(); }
/* -------------------------------------------------------------------------- */
DECLARE_COUPLER_INIT_MAKE_CALL(Kobayashi, domA, domC)
DECLARE_COUPLING_MAKE_CALL(Kobayashi, domA, domC)
__END_LIBMULTISCALE__
diff --git a/src/coupling/xiao.cc b/src/coupling/xiao.cc
index e50ca16..a823ba7 100644
--- a/src/coupling/xiao.cc
+++ b/src/coupling/xiao.cc
@@ -1,348 +1,352 @@
/**
* @file xiao.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jan 15 17:00:43 2014
*
* @brief Bridging Domain/Arlequin method
*
* @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 "xiao.hh"
#include "lm_common.hh"
#include "stimulation_zero.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
Xiao::Xiao(const std::string &name) : LMObject(name), ArlequinTemplate(name) {
this->createOutput("unmatched-point-boundary");
this->createOutput("unmatched-point-bridging");
this->createOutput("unmatched-mesh-boundary");
this->createOutput("unmatched-mesh-bridging");
this->createOutput("matched-point-boundary");
this->createOutput("matched-point-bridging");
this->createOutput("matched-mesh-boundary");
this->createOutput("matched-mesh-bridging");
this->createOutput("rhs");
this->createOutput("A");
}
/* -------------------------------------------------------------------------- */
Xiao::~Xiao() {}
/* -------------------------------------------------------------------------- */
void Xiao::coupling(CouplingStage stage) {
if (!this->is_in_continuum() && !this->is_in_atomic())
return;
if (stage == COUPLING_STEP2) {
DUMP("Zero boundary condition", DBG_INFO);
STARTTIMER("ZERO_BC");
boundary_zone.setPointField(_force, 0., spatial_dimension);
STOPTIMER("ZERO_BC");
}
if (stage != COUPLING_STEP3)
return;
STARTTIMER("BuildRHS");
DUMP("Build RHS", DBG_INFO);
this->buildRHS(_velocity);
STOPTIMER("BuildRHS");
DUMP("Synch with migration bridging", DBG_INFO);
STARTTIMER("updatesForMigrations");
bridging_zone.updateForMigration();
// this->size_constraint = bridging_zone.getNumberLocalMatchedPoints();
DUMP("Synch with migration boundary", DBG_INFO);
boundary_zone.updateForMigration();
STOPTIMER("updatesForMigrations");
DUMP("Correcting surface effect", DBG_INFO);
STARTTIMER("projectAtomicVelocitiesOnMe.hh");
// Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_displacement);
boundary_zone.projectAtomicFieldOnMesh(_velocity);
STOPTIMER("projectAtomicVelocitiesOnMe.hh");
STARTTIMER("BuildRHS");
DUMP("Build RHS", DBG_INFO);
this->buildRHS(_velocity);
STOPTIMER("BuildRHS");
DUMP("solving constraint", DBG_INFO);
STARTTIMER("Solving constraint");
this->solveConstraint();
STOPTIMER("Solving constraint");
DUMP("Applying constraint", DBG_INFO);
STARTTIMER("Correcting");
this->applyCorrection(_velocity);
STOPTIMER("Correcting");
DUMP(this->getID() << " : coupling stage all done", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
void Xiao::solveConstraint() {
for (auto &&[r, a] : zip(rhs.rowwise(), A)) {
r /= a;
}
}
/* -------------------------------------------------------------------------- */
void Xiao::applyCorrection(FieldType field) {
if (this->is_in_continuum()) {
auto &correction = this->bridging_zone.buffer_for_nodes;
correction.resize(rhs.rows(), rhs.cols());
correction = -(this->bridging_zone.smatrix * this->rhs.matrix()).array();
for (auto &&[corr, lbda] : zip(correction.rowwise(), this->lambdas_mesh)) {
corr /= lbda;
}
this->bridging_zone.addMeshField(field, correction);
}
if (this->is_in_atomic()) {
auto &correction = this->bridging_zone.buffer_for_points;
correction.resize(rhs.rows(), rhs.cols());
for (auto &&[cor, r, lbda] :
zip(correction.rowwise(), rhs.rowwise(), lambdas_point)) {
LM_ASSERT(lbda, "weight associated with atom is zero : abort");
cor = r / lbda;
}
this->bridging_zone.addPointField(field, correction);
}
}
/* -------------------------------------------------------------------------- */
void Xiao::buildRHS(FieldType field_type) {
DUMP("Clean RHS", DBG_INFO);
if (this->rhs.rows() == 0)
this->rhs.resize(0, spatial_dimension);
this->rhs.setZero();
if (this->is_in_continuum()) {
this->bridging_zone.mesh2Point(field_type, this->rhs);
}
if (this->is_in_atomic()) {
ComputeExtract field(this->getID());
field.setParam("FIELD", field_type);
field.compute(this->bridging_zone.pointList);
if (not this->is_in_continuum()) {
this->rhs = field.evalArrayOutput().array();
} else {
this->rhs -= field.evalArrayOutput().array();
}
}
bridging_zone.synchronizeVectorBySum("rhs", this->rhs);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Xiao::init(DomainA &domA, DomainC &domC) {
ArlequinTemplate::init(domA, domC);
// allocate the bridging zone (parallel version)
bridging_zone.setParam("GEOMETRY", this->bridging_geom);
bridging_zone.setParam("CHECK_COHERENCY", this->check_coherency);
bridging_zone.setParam("CENTROID", this->centroid_flag);
// set the pbc pairs
bridging_zone.setPBCPairs(domC.getPBCpairs());
// initialize the bridging_zone object
bridging_zone.init(domA, domC);
// allocate the vectors necessary to continue
- if (this->is_in_atomic() or bridging_zone.getNumberLocalMatchedPoints()) {
+ if (this->is_in_atomic() // or bridging_zone.getNumberLocalMatchedPoints()
+ ) {
this->bridging_zone.attachVector(A);
// this->allocate(bridging_zone.getNumberLocalMatchedPoints());
}
// build weights
- this->computeWeights<dispatch>(bridging_zone.pointList,
- bridging_zone.meshList);
+ if (this->is_in_atomic())
+ this->computeAtomicWeights<dispatch>(bridging_zone.pointList);
+
+ if (this->is_in_continuum())
+ this->computeContinuumWeights<dispatch>(bridging_zone.meshList);
/* build constraint matrix */
this->buildConstraintMatrix();
/* now treat the boundary zone */
boundary_zone.setParam("GEOMETRY", this->boundary_geom);
boundary_zone.setParam("CHECK_COHERENCY", this->check_coherency);
boundary_zone.setParam("CENTROID", this->centroid_flag);
// initialize the bridging_zone object
boundary_zone.init(domA, domC);
// this->boundary_zone.projectAtomicDOFsOnMesh(_displacement);
boundary_zone.projectAtomicFieldOnMesh(_velocity);
boundary_zone.setPointField(_force, 0., spatial_dimension);
this->getOutput("unmatched-point-boundary") =
boundary_zone.unmatchedPointList;
this->getOutput("unmatched-point-bridging") =
bridging_zone.unmatchedPointList;
this->getOutput("unmatched-mesh-boundary") = boundary_zone.unmatchedMeshList;
this->getOutput("unmatched-mesh-bridging") = bridging_zone.unmatchedMeshList;
this->getOutput("matched-point-boundary") = boundary_zone.pointList;
this->getOutput("matched-point-bridging") = bridging_zone.pointList;
this->getOutput("matched-mesh-boundary") = boundary_zone.meshList;
this->getOutput("matched-mesh-bridging") = bridging_zone.meshList;
this->getOutput("rhs") = rhs;
this->getOutput("A") = A;
}
/* -------------------------------------------------------------------------- */
void Xiao::buildConstraintMatrix() {
if (this->is_in_continuum()) {
auto &shp = this->bridging_zone.smatrix;
Array lumped_shape(shp.rows(), 1);
for (UInt i = 0; i < shp.rows(); ++i) {
lumped_shape(i) = shp.row(i).sum();
}
A = (shp.transpose() * (1. / lambdas_mesh * lumped_shape).matrix());
}
if (this->is_in_atomic()) {
if (not this->is_in_continuum())
A = Array::Zero(lambdas_point.rows(), lambdas_point.cols());
A += 1. / lambdas_point;
}
/* synch constraint matrix */
bridging_zone.synchronizeVectorBySum("A", this->A);
}
/* -------------------------------------------------------------------------- */
/* LMDESC XIAO
This class implements the BridgingDomain/Arlequin method.
Two zones are declared, one for the coupling (light blue atoms)
one for providing a stiff boundary condition to the atoms (purple atoms).
It defines a coupling area as depicted below:
.. image:: images/bridging-zone-BM.svg
In the bridging part a linear weight function is built the weight the
energies in the Arlequin method flavor.
The detailed algorithm is:
.. math::
\left\{
\begin{array}{l}
\displaystyle\hat{M}_I \ddot{\mathbf{u}}_I = - \mathbf{f}^R_I +
\sum_{k=1}^L \mathbf{\lambda}_k
\frac{\partial \mathbf{g}_k}{\partial \mathbf{u}_I} \\
\displaystyle \hat{m}_i \ddot{\mathbf{d}}_i = \mathbf{f}^R_i + \sum_{k=1}^L
\mathbf{\lambda}_k
\frac{\partial \mathbf{g}_k}{\partial \mathbf{d}_i}
\end{array}
\right.
where :math:`\alpha(X_I)` and :math:`\alpha(X_i)` are the weight associated
with nodes and atoms respectively, where
:math:`\hat{M}_I=\alpha(X_I)M_I` and :math:`\hat{m}_i=(1-\alpha(X_i))m_i`
and where the constraint multipliers are computed by solving
:math:`H \Lambda = \mathbf{g}^\star` with:
.. math::
H_{ik} = \Delta t \left( \sum_J \varphi_J(\mathbf{X}_I) \hat{M_J}^{-1}
\frac{\partial \mathbf{g}_k}{\partial \mathbf{u}_J} -
\hat{m}_I^{-1} \frac{\partial \mathbf{g}_k}{\partial \mathbf{d}_i}
\right),
\mathbf{g}_i^\star = \sum_J \varphi_J(\mathbf{X}_I)
\dot{\mathbf{u}}_J^\star
- \dot{\mathbf{d}}_i^\star,
where :math:`u^\star, d^\star` are the displacement obtained after one time
step
by the integration scheme when the constraints are not applied.
For debugging puposes several computes are registered to the central
system:
- weight-fe-\${couplerID} : it is a compute containing the weights
associated with nodes
inside of the bridging zone (the :math:`\alpha(X_I)`).
- weight-md-\${couplerID} : it is a compute containing the weights
associated with atoms
inside of the bridging zone (the :math:`\alpha(X_i)`).
- lambdas-fe-\${couplerID} : it is a compute containing the
:math:`\hat{M}_I`.
- lambdas-md-\${couplerID} : it is a compute
containing the :math:`\hat{m}_i`.
*/
/* LMEXAMPLE
COUPLING_CODE bdmID md fe XIAO GEOMETRY 3 BOUNDARY 4 GRID_DIVISIONX
10 GRID_DIVISIONY 10 QUALITY 1e-3
*/
/* LMHERITANCE arlequin_template */
void Xiao::declareParams() {
ArlequinTemplate::declareParams();
/* LMKEYWORD CENTROID
Set the selection of the bridging zones based on centroid.
*/
this->parseTag("CENTROID", centroid_flag, false);
}
/* -------------------------------------------------------------------------- */
DECLARE_COUPLER_INIT_MAKE_CALL(Xiao, domA, domC)
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__