Page MenuHomec4science
Files F67605125

xiao.cc
No OneTemporary
Actions

File Metadata

Created
Sun, Jun 23, 08:10

xiao.cc
View Options

/**
* @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()
) {
this->bridging_zone.attachVector(A);
// this->allocate(bridging_zone.getNumberLocalMatchedPoints());
}
// build weights
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__

Event Timeline

c4science · Help