Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/coupling/arlequin_template.cc
0b9417d37655features/coupler_compon…

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 "arlequin_template.hh"
#include "bridging.hh"
#include "bridging_atomic_continuum.hh"
#include "factory_multiscale.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
ArlequinTemplate::ArlequinTemplate(const std::string &name)
: LMObject(name), CouplingAtomicContinuum(name),
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);
}
/* -------------------------------------------------------------------------- */
ArlequinTemplate::~ArlequinTemplate() { clearAll(); }
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::init(DomainA &domA, DomainC &domC) {
using BridgingType = Bridging<typename DomainA::ContainerPoints,
typename DomainC::ContainerMesh>;
bridging_zone =
std::make_unique<BridgingType>(this->getID() + "-bridging", domA, domC);
MDboundary_zone =
std::make_unique<BridgingType>(this->getID() + "-boundary", domA, domC);
}
/* -------------------------------------------------------------------------- */
void ArlequinTemplate::clearAll() {
// Parent::domA.getRefManager().DetachVector(A,&this->atomes_rec);
A.clear();
rhs.clear();
}
/* -------------------------------------------------------------------------- */
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.to_array().setZero();
DUMP("Attach vector A", DBG_INFO_STARTUP);
auto &bridging_zone = *this->bridging_zone;
bridging_zone.attachVector(A);
}
/* -------------------------------------------------------------------------- */
void ArlequinTemplate::buildContinuumConstraintMatrix() {
if (size_constraint == 0)
return;
this->bridging_zone->buildContinuConstraint(
A, this->bridging_zone->node_shape, lambdasC);
}
/* -------------------------------------------------------------------------- */
template <typename ContC>
void ArlequinTemplate::computeContinuumWeights(ContC &mesh) {
this->weightFE.buildManual(mesh);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::correctContinuumWeights() {
lambdasC.assign(weightFE.getArray().size(), 0);
UInt j = 0;
auto &bridging_zone = *this->bridging_zone;
for (auto &&n :
bridging_zone.cast<DomainA, DomainC>().meshList.getContainerNodes()) {
lambdasC[j] = weightFE.getArray()[j] * n.mass();
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContA>
void ArlequinTemplate::computeAtomWeights(ContA &atoms) {
this->weightMD.buildManual(atoms);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::correctAtomWeights() {
auto &bridging_zone = *this->bridging_zone;
auto &pointList = bridging_zone.cast<DomainA, DomainC>().pointList;
lambdasA.assign(weightMD.getArray().size(), 0);
bridging_zone.attachVector(lambdasA);
UInt at_index = 0;
for (auto &&at : pointList) {
LM_ASSERT(at_index < weightMD.getArray().size(),
"overflow detected " << at_index
<< " >= " << weightMD.getArray().size());
DUMP("correcting lambdasA[" << at_index
<< "] = " << weightMD.getArray()[at_index]
<< " *= " << at.mass(),
DBG_DETAIL);
lambdasA[at_index] = weightMD.getArray()[at_index] * at.mass();
++at_index;
}
}
/* -------------------------------------------------------------------------- */
void ArlequinTemplate::buildAtomsConstraintMatrix() {
A += 1. / lambdasA.getArray().to_array();
}
/* -------------------------------------------------------------------------- */
void ArlequinTemplate::cleanRHS() {
rhs.resize(this->bridging_zone->getNumberLocalMatchedPoints(),
spatial_dimension);
rhs.to_array().setZero();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::buildContinuumRHS() {
auto &domC = this->getDomC<DomainC>();
auto field = domC.getField(_velocity);
auto &bridging_zone = *this->bridging_zone;
bridging_zone.cast<DomainA, DomainC>().getShapeMatrix().buildRHS(rhs, field);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::buildAtomsRHS() {
UInt at_index = 0;
constexpr UInt Dim = DomainA::ContainerPoints::Dim;
auto &bridging_zone = *this->bridging_zone;
for (auto &&at : bridging_zone.cast<DomainA, DomainC>().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>
void ArlequinTemplate::solveConstraint() {
constexpr UInt Dim = DomainA::ContainerPoints::Dim;
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 Field>
void ArlequinTemplate::applyContinuumCorrection(Field &&field) {
this->bridging_zone->applyContinuumConstraint(field, rhs, lambdasC);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void ArlequinTemplate::applyAtomsCorrection() {
auto &bridging_zone = this->bridging_zone->cast<DomainA, DomainC>();
LM_ASSERT(size_constraint == bridging_zone.pointList.size(),
"something bad is happening " << size_constraint << " "
<< bridging_zone.pointList.size());
UInt i = 0;
constexpr UInt Dim = DomainA::ContainerPoints::Dim;
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 */
void ArlequinTemplate::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