Homec4science
Diffusion LibMultiScale (features/pybind) Browse

libmultiscale/src/coupling/arlequin_template.cc
b05d5f2aba78features/pybind

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 <fstream>
#include "lm_common.hh"
#include "bridging.hh"
#include "bridging_atomic_continuum.hh"
#include "arlequin_template.hh"
#include "weighting.hh"
#include "lib_bridging.hh"
#include "filter_manager.hh"
/* -------------------------------------------------------------------------- */
__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,1),
weightMD("weight-md-" + name,1),
lambdasC("lambdas-fe-" + name,1),
lambdasA("lambdas-md-" + name,1)
{
quality = 1e-3;
multi_time_step = 1;
size_constraint = 0;
/* registering computes for outer world */
FilterManager::getManager().addObject(&this->weightFE,false);
FilterManager::getManager().addObject(&this->lambdasC,false);
FilterManager::getManager().addObject(&this->weightMD,false);
FilterManager::getManager().addObject(&this->lambdasA,false);
}
/* -------------------------------------------------------------------------- */
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.nbElem();
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.nbElem()
<< " elements",DBG_INFO);
typename Parent::IteratorNodesSubset itrec = bridging_zone.meshList.getIterator();
typename Parent::RefNode n;
UInt j = 0;
Real pos0[3] = {0.0,0.0,0.0};
for (n = itrec.getFirst(); !itrec.end() ;
n = itrec.getNext() , ++j) {
n.getPositions0(pos0);
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
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void ArlequinTemplate<DomainA,DomainC,Dim>::correctContinuumWeights(){
typename Parent::IteratorNodesSubset itrec =
bridging_zone.meshList.getContainerNodes().getIterator();
lambdasC.assign(weightFE.nbElem(),0);
typename Parent::RefNode n;
UInt j = 0;
for (n = itrec.getFirst() , j = 0; !itrec.end() ; n = itrec.getNext() , ++j ) {
lambdasC[j] = weightFE[j]*n.mass();
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void ArlequinTemplate<DomainA,DomainC,Dim>::computeAtomWeights(){
Weighting<ATOMEFLAG,Dim,LINEAR> P(bridging_geom);
typename Parent::ContainerAtomsSubset & pointList = bridging_zone.pointList;
UInt nb = pointList.nbElem();
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);
typename Parent::IteratorAtomsSubset itrec = pointList.getIterator();
typename Parent::RefAtom at = itrec.getFirst();
UInt at_index = 0;
Real pos0[3] = {0.0,0.0,0.0};
for (at = itrec.getFirst(),at_index=0;!itrec.end();
at = itrec.getNext(),++at_index) {
LM_ASSERT(at_index < weightMD.nbElem(),
"overflow detected " << at_index << " >= "
<< weightMD.nbElem());
at.getPositions0(pos0);
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
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void ArlequinTemplate<DomainA,DomainC,Dim>::correctAtomWeights(){
typename Parent::ContainerAtomsSubset &
pointList = bridging_zone.pointList;
lambdasA.assign(weightMD.nbElem(),0);
bridging_zone.attachVector(lambdasA);
typename Parent::IteratorAtomsSubset itrec = pointList.getIterator();
typename Parent::RefAtom at = itrec.getFirst();
UInt at_index = 0;
for (at = itrec.getFirst(),at_index=0;!itrec.end();
at = itrec.getNext(),++at_index) {
LM_ASSERT(at_index < weightMD.nbElem(),
"overflow detected " << at_index
<< " >= " << weightMD.nbElem());
DUMP("correcting lambdasA[" << at_index << "] = " << weightMD[at_index]
<< " *= " << at.mass(),DBG_DETAIL);
lambdasA[at_index] = weightMD[at_index]*at.mass();
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void ArlequinTemplate<DomainA,DomainC,Dim>::buildAtomsConstraintMatrix(){
typename Parent::IteratorAtomsSubset itrec = bridging_zone.pointList.getIterator();
typename Parent::RefAtom at;
UInt i = 0;
for (itrec.getFirst(); !itrec.end() ; at = itrec.getNext(), ++i) {
LM_ASSERT(i < lambdasA.nbElem(),
"overflow detected " << i << " >= " << lambdasA.nbElem());
LM_ASSERT(i < A.size(),
"overflow detected " << i << " >= " << A.size());
LM_ASSERT(!isinf(A[i]) && !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(){
typename Parent::IteratorAtomsSubset itrec = bridging_zone.pointList.getIterator();
typename Parent::RefAtom at = itrec.getFirst();
UInt at_index=0;
for (at = itrec.getFirst(),at_index=0;!itrec.end();
at = itrec.getNext(),++at_index) {
LM_ASSERT(at_index*Dim < rhs.size(),
"overflow detected " << at_index*Dim
<< " >= " << rhs.size());
for (UInt k = 0; k < Dim; ++k) {
DUMP("rhs[" << at_index << "]=" << rhs[Dim*at_index+k]
<< " before correction",DBG_ALL);
LM_ASSERT(!isnan(at.velocity(k)) && !isinf(at.velocity(k)),
"problem with this atom : " << at
<< "velocity(" << k << ") is not correct : "
<< at.velocity(k)
<< "at_index = " << at_index << " ");
rhs[Dim*at_index+k] -= at.velocity(k);
DUMP(at_index << " is adding " << -1.0*at.velocity(k) << " to rhs["
<< at_index << "] = " << rhs[Dim*at_index+k],DBG_ALL);
}
}
}
/* -------------------------------------------------------------------------- */
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(!isnan(rhs[Dim*i]) && !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(!isnan(rhs[Dim*i]) && !isnan(rhs[Dim*i+1])
&& !isinf(rhs[Dim*i]) && !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(!isnan(rhs[Dim*i]) && !isnan(rhs[Dim*i+1])
&& !isnan(rhs[Dim*i+2]) && !isinf(rhs[Dim*i])
&& !isinf(rhs[Dim*i+1]) && !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(){
typename Parent::_Vec_ & 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(){
typename Parent::IteratorAtomsSubset it = bridging_zone.pointList.getIterator();
typename Parent::RefAtom at = it.getFirst();
UInt i = 0;
LM_ASSERT(size_constraint == bridging_zone.pointList.nbElem(),
"something bad is happening "
<< size_constraint << " "
<< bridging_zone.pointList.nbElem());
for (at = it.getFirst(); !it.end();at = it.getNext(),++i)
{
LM_ASSERT(i < lambdasA.nbElem(),
"overflow detected " << i << " >= " << lambdasA.nbElem());
LM_ASSERT(i*Dim < rhs.size(),
"overflowdetected " << i*Dim << " >= " << rhs.size());
DUMP("correcting atom " << i << "/(" << size_constraint << ") on "
<< at.velocity(0) << " by adding "
<< rhs[i*Dim+0]/lambdasA[i],DBG_DETAIL);
LM_ASSERT(lambdasA[i],
"weight associated with atom " << i << " is zero : abort");
at.velocity(0) += rhs[i*Dim]/lambdasA[i];
if (Dim > 1)
at.velocity(1) += rhs[i*Dim+1]/lambdasA[i];
if (Dim == 3)
at.velocity(2) += rhs[i*Dim+2]/lambdasA[i];
// Srinivasa: Fixed minor bug by introducing following 3 if conditions
if(Dim == 1){
LM_ASSERT(!isnan(at.velocity(0)) && !isinf(at.velocity(0)),
"problem with this atom : " << at
<< "velocity value is not correct : " << at.velocity(0) << " "
<< "applied correction was : " << rhs[i*Dim] << " "
<< "weight was : " << lambdasA[i]
<< "at_index = " << i << " ");
}
if(Dim == 2){
LM_ASSERT(!isnan(at.velocity(0)) && !isnan(at.velocity(1))
&& !isinf(at.velocity(0)) && !isinf(at.velocity(1)),
"problem with this atom : " << at
<< "velocity value is not correct : " << at.velocity(0) << " "
<< at.velocity(1) << " "
<< "applied correction was : " << rhs[i*Dim] << " "
<< rhs[i*Dim+1] << " "
<< "weight was : " << lambdasA[i]
<< "at_index = " << i << " ");
}
if(Dim == 3){
LM_ASSERT(!isnan(at.velocity(0)) && !isnan(at.velocity(1))
&& !isnan(at.velocity(2)) && !isinf(at.velocity(0))
&& !isinf(at.velocity(1)) && !isinf(at.velocity(2)),
"problem with this atom : " << at
<< "velocity value is not correct : " << at.velocity(0) << " "
<< at.velocity(1) << " " << at.velocity(2)
<< "applied correction was : " << rhs[i*Dim] << " "
<< rhs[i*Dim+1] << " " << rhs[i*Dim+2]
<< "weight was : " << lambdasA[i]
<< "at_index = " << 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);
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(ArlequinTemplate)
__END_LIBMULTISCALE__
c4science · Help