Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/coupling/arlequin_template.cc
27ca80ef7f29features/reshape_comput…

Last Change
Raw File

arlequin_template.cc
View Options

/**
* @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.
*
*/
#define ZERO_LIMIT 1e-3
/* -------------------------------------------------------------------------- */
#include "factory_multiscale.hh"
#include "arlequin_template.hh"
#include "bridging.hh"
#include "bridging_atomic_continuum.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
#include "weighting.hh"
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename DomainA, typename DomainC, UInt Dim>
ArlequinTemplate<DomainA, DomainC, Dim>::ArlequinTemplate(
const std::string &name, DomainInterface &domA, DomainInterface &domC)
: CouplingAtomicContinuum<DomainA, DomainC, Dim>(name, domA, domC),
bridging_zone(name + "-bridging", Parent::domA, Parent::domC),
MDboundary_zone(name + "-boundary", Parent::domA, Parent::domC),
weightFE("weight-fe-" + name), weightMD("weight-md-" + name),
lambdasC("lambdas-fe-" + name), lambdasA("lambdas-md-" + name) {
quality = 1e-3;
multi_time_step = 1;
size_constraint = 0;
/* registering computes for outer world */
FilterManager::getManager().addObject(this->weightFE);
FilterManager::getManager().addObject(this->lambdasC);
FilterManager::getManager().addObject(this->weightMD);
FilterManager::getManager().addObject(this->lambdasA);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
ArlequinTemplate<DomainA, DomainC, Dim>::~ArlequinTemplate() {
clearAll();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::clearAll() {
// Parent::domA.getRefManager().DetachVector(A,&this->atomes_rec);
A.clear();
rhs.clear();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::allocate(UInt t) {
size_constraint = t;
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.assign((size_constraint)*Dim, 0);
DUMP("Attach vector A", DBG_INFO_STARTUP);
bridging_zone.attachVector(A);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::buildContinuConstraintMatrix() {
if (size_constraint == 0)
return;
bridging_zone.getShapeMatrix().buildContinuConstraint(
A, bridging_zone.node_shape, lambdasC);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::computeContinuumWeights() {
//! weighting object (used to compute the weights of each dofs)
Weighting<CONTINUFLAG, Dim, LINEAR> P(bridging_geom);
UInt nb = bridging_zone.meshList.size();
weightFE.assign(nb, 0);
if (!nb) {
DUMP("There are no nodes in the bridging", DBG_WARNING);
return;
}
#ifndef LM_OPTIMIZED
typename Parent::ContainerElemsSubset &elemList =
bridging_zone.meshList.getContainerElems();
#endif // LM_OPTIMIZED
DUMP("We found " << nb << " nodes concerned in " << elemList.size()
<< " elements",
DBG_INFO);
UInt j = 0;
for (auto &&n : bridging_zone.meshList) {
auto pos0 = n.position0();
weightFE[j] = P.weight(pos0);
if (weightFE[j] < ZERO_LIMIT)
weightFE[j] = quality;
DUMP("weightFE[" << j << "]=" << weightFE[j], DBG_ALL);
#ifdef LIBMULTISCALE_USE_MECA1D
if (dynamic_cast<DomainMeca1D *>(&(this->domC)) &&
Parent::domC.arlequin() == true) {
n.alpha() = weightFE[j];
}
#endif
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::correctContinuumWeights() {
lambdasC.assign(weightFE.size(), 0);
UInt j = 0;
for (auto &&n : bridging_zone.meshList.getContainerNodes()) {
lambdasC[j] = weightFE[j] * n.mass();
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::computeAtomWeights() {
Weighting<ATOMEFLAG, Dim, LINEAR> P(bridging_geom);
auto & pointList = bridging_zone.pointList;
UInt nb = pointList.size();
weightMD.resize(nb);
bridging_zone.attachVector(weightMD);
if (!nb) {
DUMP("We found no atoms in the bridging zone", DBG_WARNING);
return;
}
DUMP("We found " << nb << " concerned atoms", DBG_INFO);
UInt at_index = 0;
for (auto &&at : pointList) {
LM_ASSERT(at_index < weightMD.size(),
"overflow detected " << at_index << " >= " << weightMD.size());
Vector<Dim> pos0 = at.position0();
weightMD[at_index] = P.weight(pos0);
if (weightMD[at_index] < ZERO_LIMIT)
weightMD[at_index] = quality;
DUMP("weightMD[" << at_index << "]=" << weightMD[at_index], DBG_ALL);
#ifdef LIBMULTISCALE_USE_MD1D
if (dynamic_cast<DomainMD1D *>(&(this->domA)) &&
Parent::domA.arlequin() == true) {
at.alpha() = weightMD[at_index];
}
#endif
++at_index;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::correctAtomWeights() {
auto &pointList = bridging_zone.pointList;
lambdasA.assign(weightMD.size(), 0);
bridging_zone.attachVector(lambdasA);
UInt at_index = 0;
for (auto &&at : pointList) {
LM_ASSERT(at_index < weightMD.size(),
"overflow detected " << at_index << " >= " << weightMD.size());
DUMP("correcting lambdasA[" << at_index << "] = " << weightMD[at_index]
<< " *= " << at.mass(),
DBG_DETAIL);
lambdasA[at_index] = weightMD[at_index] * at.mass();
++at_index;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::buildAtomsConstraintMatrix() {
for (UInt i = 0; i < bridging_zone.pointList.size(); ++i) {
LM_ASSERT(i < lambdasA.size(), "overflow detected "
<< i << " >= " << lambdasA.size());
LM_ASSERT(i < A.size(), "overflow detected " << i << " >= " << A.size());
LM_ASSERT(!std::isinf(A[i]) && !std::isnan(lambdasA[i]),
"problem with this constraUInt index : "
<< i << "constraint value and atom weight are : " << A[i]
<< " " << lambdasA[i] << "at_index = " << i << " ");
DUMP("compute A [" << i << "] = " << A[i] << " += 1/" << lambdasA[i]
<< " := " << 1 / lambdasA[i],
DBG_ALL);
A[i] += 1 / lambdasA[i];
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::cleanRHS() {
rhs.assign(bridging_zone.getLocalPoints() * Dim, 0);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::buildContinuRHS() {
typename Parent::_Vec_ &field = Parent::domC.getV();
bridging_zone.getShapeMatrix().buildRHS(rhs, field);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::buildAtomsRHS() {
UInt at_index = 0;
for (auto &&at : bridging_zone.pointList) {
LM_ASSERT(at_index * Dim < rhs.size(),
"overflow detected " << at_index * Dim << " >= " << rhs.size());
DUMP("rhs[" << at_index << "]=" << VectorView<Dim>(&rhs[Dim * at_index])
<< " before correction",
DBG_ALL);
LM_ASSERT(!isNaN(at.velocity()) && !isInf(at.velocity()),
"problem with this atom : " << at << "velocity is not correct : "
<< at.velocity()
<< "at_index = " << at_index << " ");
VectorView<Dim>(&rhs[Dim * at_index]) -= at.velocity();
DUMP(at_index << " is adding " << -1.0 * at.velocity() << " to rhs["
<< at_index
<< "] = " << VectorView<Dim>(&rhs[Dim * at_index]),
DBG_ALL);
++at_index;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::solveConstraint() {
for (UInt i = 0; i < size_constraint; ++i) {
LM_ASSERT(i < A.size(), "overflow detected " << i << " >= " << A.size());
LM_ASSERT(i * Dim < rhs.size(), "overflow detected "
<< i * Dim << " >= " << rhs.size());
DUMP("multL[" << i << "] = " << rhs[Dim * i + 1] / A[i] << " = "
<< rhs[Dim * i + 1] << "/" << A[i],
DBG_ALL);
rhs[Dim * i] = rhs[Dim * i] / A[i];
if (Dim > 1)
rhs[Dim * i + 1] = rhs[Dim * i + 1] / A[i];
if (Dim == 3)
rhs[Dim * i + 2] = rhs[Dim * i + 2] / A[i];
if (Dim == 1) {
LM_ASSERT(!std::isnan(rhs[Dim * i]) && !std::isinf(rhs[Dim * i]),
"problem with this index : "
<< i << "ConstraUInt matrix value is : " << A[i]
<< "applied correction was : " << rhs[i * Dim]
<< "at_index = " << i << " ");
}
if (Dim == 2) {
LM_ASSERT(!std::isnan(rhs[Dim * i]) && !std::isnan(rhs[Dim * i + 1]) &&
!std::isinf(rhs[Dim * i]) && !std::isinf(rhs[Dim * i + 1]),
"problem with this index : "
<< i << "ConstraUInt matrix value is : " << A[i]
<< "applied correction was : " << rhs[i * Dim] << " "
<< rhs[i * Dim + 1] << " at_index = " << i << " ");
}
if (Dim == 3) {
LM_ASSERT(
!std::isnan(rhs[Dim * i]) && !std::isnan(rhs[Dim * i + 1]) &&
!std::isnan(rhs[Dim * i + 2]) && !std::isinf(rhs[Dim * i]) &&
!std::isinf(rhs[Dim * i + 1]) && !std::isinf(rhs[Dim * i + 2]),
"problem with this index : "
<< i << "ConstraUInt matrix value is : " << A[i]
<< "applied correction was : " << rhs[i * Dim] << " "
<< rhs[i * Dim + 1] << " " << rhs[i * Dim + 2]
<< "at_index = " << i << " ");
}
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::applyContinuCorrection() {
auto &field = Parent::domC.getV();
bridging_zone.getShapeMatrix().applyChanging(field, rhs, lambdasC);
field.close();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::applyAtomsCorrection() {
LM_ASSERT(size_constraint == bridging_zone.pointList.size(),
"something bad is happening " << size_constraint << " "
<< bridging_zone.pointList.size());
UInt i = 0;
for (auto &&at : bridging_zone.pointList) {
LM_ASSERT(i < lambdasA.size(), "overflow detected "
<< i << " >= " << lambdasA.size());
LM_ASSERT(i * Dim < rhs.size(), "overflowdetected "
<< i * Dim << " >= " << rhs.size());
DUMP("correcting atom " << i << "/(" << size_constraint << ") on "
<< at.velocity() << " by adding "
<< VectorView<Dim>(&rhs[i * Dim]) / lambdasA[i],
DBG_DETAIL);
LM_ASSERT(lambdasA[i], "weight associated with atom "
<< i << " is zero : abort");
at.velocity() += VectorView<Dim>(&rhs[i * Dim]) / lambdasA[i];
LM_ASSERT(
!isNaN(at.velocity()) && !isInf(at.velocity()),
"problem with this atom : "
<< at << "velocity value is not correct : " << at.velocity() << " "
<< "applied correction was : " << VectorView<Dim>(&rhs[i * Dim])
<< " "
<< "weight was : " << lambdasA[i] << "at_index = " << i << " ");
++i;
}
}
/* -------------------------------------------------------------------------- */
/* 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 */
template <typename DomainA, typename DomainC, UInt Dim>
void ArlequinTemplate<DomainA, DomainC, Dim>::declareParams() {
this->addSubParsableObject(bridging_zone);
this->addSubParsableObject(MDboundary_zone);
/* LMKEYWORD QUALITY
Because of the strong sense brought in the formulation
of the Lagrange constraints zero weights are prohibited.
Because geometrical situation can lead to a zero weight,
the quality factor defines the replacement value for these zeros.
More details on the impact of the factor can be found in \\
\textit{Ghost force reduction and spectral analysis of the 1D bridging
method}\\
\textbf{Guillaume Anciaux, Olivier Coulaud, Jean Roman, Gilles Zerah}\\
\url{http://hal.inria.fr/inria-00300603/en/}
*/
this->parseKeyword("QUALITY", quality);
/* 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.
\textbf{Be careful, it is extremely computationally expensive}
*/
this->parseTag("CHECK_COHERENCY", check_coherency, false);
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(ArlequinTemplate)
__END_LIBMULTISCALE__
c4science · Help