Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/filter/compute_arlequin_weight.cc
d8e9d1cc530clammps_2020

Last Change
Raw File

compute_arlequin_weight.cc
View Options

/**
* @file compute_arlequin_weight.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @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 "compute_arlequin_weight.hh"
#include "lib_continuum.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
///* --------------------------------------------------------------------------
///*/
// enum FunctionOrder { LINEAR = 1, THIRD = 3, CUSTOM = 4 };
///* -------------------------------------------------------------------------
///*/
//
// struct Weight1D {
// template <FunctionOrder order> static Real weight(Real x, Real R);
//};
/* ------------------------------------------------------------------------ */
Real ComputeArlequinWeight::computeWeight(Real x, Real R) { return x / R; }
/* ------------------------------------------------------------------------ */
// template <> inline Real Weight1D::weight<THIRD>(Real x, Real R) {
// Real R3 = R * R * R;
// return (-1. * x * x * (2. * x - 3 * R) / R3);
//}
//
///* --------------------------------------------------------------------------
///*/
// template <> inline Real Weight1D::weight<CUSTOM>(Real x, Real R) {
// /* Real r=R,q=8*R/10,epsilon=0.1; */
// /* Real a = -(epsilon*r+q-r)/(q-r)/(q-r)/q; */
// /* Real b = (q*q*epsilon -r*r +q*q +epsilon*r*r)/(q-r)/(q-r)/q; */
// /* Real c = - (q+epsilon*q-r)*r/(q-r)/(q-r); */
// /* if(x < q) */
// /* return (1-epsilon)/q*x; */
// /* return (1-a*x*x - b*x -c); */
// /* ********************** autre experience **************/
// /* Real epsilon = .1; */
// /* // Real q = R/16;//marche pas */
// /* Real q = R/32; */
// /* if (x < q) */
// /* return 0; */
// /* return 1/(R-q)*(x-q); */
//
// return (x / R);
//}
//
///* --------------------------------------------------------------------------
///*/
//
// template <UInt FLAG, UInt Dim, FunctionOrder order>
// inline Weighting<FLAG, Dim, order>::Weighting(LMID geom_id) {
// this->geom = GeometryManager::getManager().getGeometry(geom_id);
//}
/* ------------------------------------------------------------------------ */
// template <UInt FLAG, UInt Dim, FunctionOrder order>
// template <typename T> inline Real ComputeArlequinWeight::computeWeight(T &&X)
// {
// Real result = -1;
// if (geom == NULL) {
// std::stringstream sstr;
// sstr << "weighting of point at coordinate ";
// for (UInt k = 0; k < Dim; ++k)
// sstr << X[k] << " ";
// LM_FATAL(sstr.str() << " impossible since passed geometry is null");
// }
// Geometry &g = *geom;
// switch (g.getType()) {
// case Geometry::CUBE:
// result = weight(X, dynamic_cast<Cube &>(*geom));
// break;
// case Geometry::BALL:
// result = weight(X, dynamic_cast<Ball &>(*geom));
// break;
// default:
// LM_TOIMPLEMENT;
// }
// if (result > 1.0)
// result = 1.0;
// if (result < 0.0)
// result = 0.0;
// if (FLAG == ATOMEFLAG)
// return (1.0 - result);
// return result;
// }
///* --------------------------------------------------------------------------
///*/
//
///*
///--------------------------------------------------------------------------
///*/
//
// template <UInt FLAG, UInt Dim, FunctionOrder order>
// inline Real Weighting<FLAG, Dim, order>::computeWeight(Real x, Real R) {
// return Weight1D::weight<order>(x, R);
//}
//
//
/* --------------------------------------------------------------------------
*/
template <UInt Dim> inline Real ComputeArlequinWeight::weight(Vector<Dim> X) {
auto &&geom = GeometryManager::getGeometry(this->geom);
switch (geom->getType()) {
case Geometry::BALL:
return this->weight(X, *std::dynamic_pointer_cast<Ball>(geom));
return 0.;
case Geometry::CUBE:
return this->weight(X, *std::dynamic_pointer_cast<Cube>(geom));
default:
LM_TOIMPLEMENT;
}
}
/* ------------------------------------------------------------------------ */
template <UInt Dim>
inline Real ComputeArlequinWeight::weight(Vector<Dim> X, Ball &g) {
Real result = 0;
DUMP("employed gemometry " << g, DBG_DETAIL);
auto x = X - g.getCenter().cast<Real>().block<Dim, 1>(0, 0);
Real d = sqrt(x.dot(x));
DUMP("after correction x = " << x[0] << " y = " << x[1] << " z = " << x[2]
<< " d = " << d,
DBG_DETAIL);
result = computeWeight(d - g.rMin(), g.rMax() - g.rMin());
return result;
}
/* ------------------------------------------------------------------------ */
template <UInt Dim>
inline Real ComputeArlequinWeight::weight(Vector<Dim> X, Cube &c) {
Real result = 0;
if (c.getXmax()[1] < 0) {
Real R = c.getXmax()[1] - c.getXmin()[1];
result = computeWeight(c.getXmax()[1] - X[1], R);
} else if (c.getXmin()[1] > 0) {
Real R = c.getXmax()[1] - c.getXmin()[1];
result = computeWeight(X[1] - c.getXmin()[1], R);
} else if (c.getXmax()[2] < 0) {
Real R = c.getXmax()[2] - c.getXmin()[2];
result = computeWeight(c.getXmax()[2] - X[2], R);
} else if (c.getXmin()[2] > 0) {
Real R = c.getXmax()[2] - c.getXmin()[2];
result = computeWeight(X[2] - c.getXmin()[2], R);
} else if (c.getXmax()[0] < 0) {
Real R = c.getXmax()[0] - c.getXmin()[0];
result = computeWeight(c.getXmax()[0] - X[0], R);
} else if (c.getXmin()[0] > 0) {
Real R = c.getXmax()[0] - c.getXmin()[0];
result = computeWeight(X[0] - c.getXmin()[0], R);
} else
LM_FATAL("bridged zone can not be of size 0");
return result;
}
/* ------------------------------------------------------------------------ */
ComputeArlequinWeight::ComputeArlequinWeight(const LMID &id) : LMObject(id) {
this->createInput("input");
this->createOutput("weight");
this->createArrayOutput("weight");
}
/* --------------------------------------------------------------------------
*/
ComputeArlequinWeight::~ComputeArlequinWeight() {}
/* ------------------------------------------------------------------------ */
template <typename Cont>
enable_if_mesh<Cont> ComputeArlequinWeight::build(Cont &meshList) {
constexpr UInt Dim = Cont::Dim;
auto &weightFE = this->getOutputAsArray("weight");
UInt nb = meshList.size();
weightFE.assign(nb, 0);
if (!nb) {
DUMP("There are no nodes in the bridging", DBG_WARNING);
return;
}
#ifndef LM_OPTIMIZED
auto &elemList = meshList.getContainerElems();
DUMP("We found " << nb << " nodes concerned in " << elemList.size()
<< " elements",
DBG_INFO);
#endif // LM_OPTIMIZED
for (auto &&[n, w] : zip(meshList, weightFE)) {
Vector<Dim> pos0 = n.position0();
w = this->weight(pos0);
if (w > 1.)
w = 1.;
if (w < ZERO_LIMIT)
w = quality;
DUMP("w = " << w, DBG_ALL);
}
}
/* ------------------------------------------------------------------------ */
template <typename Cont>
enable_if_point<Cont> ComputeArlequinWeight::build(Cont &pointList) {
constexpr UInt Dim = Cont::Dim;
UInt nb = pointList.size();
auto &weightMD = this->getOutputAsArray("weight");
weightMD.assign(nb, 0);
if (!nb) {
DUMP("We found no atoms in the bridging zone", DBG_WARNING);
return;
}
#ifndef LM_OPTIMIZED
DUMP("We found " << nb << " concerned atoms", DBG_INFO);
//
#endif // LM_OPTIMIZED
for (auto &&[at, w] : zip(pointList, weightMD)) {
Vector<Dim> pos0 = at.position0();
w = 1. - this->weight(pos0);
if (w < ZERO_LIMIT)
w = quality;
DUMP("w = " << w, DBG_ALL);
}
}
/* --------------------------------------------------------------------------
*/
/* LMDESC ARLEQUINWEIGHT
This computes a typical weight function, usable for arlequin
coupling. Beware there are strong geometry assumptions (like syymetry).
*/
/* LMEXAMPLE
COMPUTE weight ARLEQUIN_WEIGHT INPUT md GEOM some_geometry
*/
inline void ComputeArlequinWeight::declareParams() {
ComputeInterface::declareParams();
/* LMKEYWORD QUALITY
Because of the strong sense brought in the formulation
of the Lagrange constraints zero weights are prohibited.
Therefore the so-called quality factor defines the
replacement value for these zeros.
More details on the impact of the factor can be found in
*Ghost force reduction and spectral analysis of the 1D bridging
method*.
**Guillaume Anciaux, Olivier Coulaud, Jean Roman, Gilles Zerah.**
`(Link) <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", geom);
}
/* --------------------------------------------------------------------------
*/
DECLARE_COMPUTE_MAKE_CALL(ComputeArlequinWeight)
__END_LIBMULTISCALE__
c4science · Help