Page MenuHomec4science
Files F64200099

kobayashi.cc
No OneTemporary
Actions

File Metadata

Created
Sat, May 25, 06:47

kobayashi.cc
View Options

/**
* @file kobayashi.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief Kobayashi's bridging method
*
* @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 "kobayashi.hh"
#include "arlequin_template.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
Kobayashi::Kobayashi(const std::string &name)
: LMObject(name), ArlequinTemplate(name) {
// naming the bridging zones
std::stringstream fname;
fname << name << "-bridging";
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC> void Kobayashi::init() {
ArlequinTemplate::init<DomainA, DomainC>();
auto &bridging_zone = *this->bridging_zone;
auto &MDboundary_zone = *this->MDboundary_zone;
if (lm_world_size > 1)
LM_FATAL("this coupler work only in sequential");
std::stringstream fname;
// initialize the bridging_zone object
bridging_zone.init();
// allocate the vectors necessary to continue
UInt taille = bridging_zone.cast<DomainA, DomainC>().pointList.size();
this->allocate(taille);
// build weights
this->computeContinuumWeights<DomainA, DomainC>();
this->computeAtomWeights<DomainA, DomainC>();
UInt _taille = bridging_zone.cast<DomainA, DomainC>().meshList.size();
DUMP("Allocating vectors", DBG_INFO);
rhs = new Array(taille);
CG_vel = new Array(_taille);
CG_disp = new Array(_taille);
DUMP("Computing initial least square constraint", DBG_INFO);
buildCGMatrix<DomainA, DomainC>();
// treat now MD boundary zone
MDboundary_zone.init();
DUMP("Init ... done", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::coupling(CouplingStage stage) {
auto &bridging_zone = *this->bridging_zone;
auto &MDboundary_zone = *this->MDboundary_zone;
if (stage != COUPLING_STEP4)
return;
MDboundary_zone.projectAtomicDOFsOnMesh(_velocity);
MDboundary_zone.projectAtomicDOFsOnMesh(_displacement);
buildLeastSquareRHS<DomainA, DomainC>();
correctContinuum<DomainA, DomainC>();
correctAtoms<DomainA, DomainC>();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC> void Kobayashi::buildCGMatrix() {
auto &bridging_zone = *this->bridging_zone;
UInt nb_nodes = bridging_zone.cast<DomainA, DomainC>().meshList.size();
UInt nb_atoms = bridging_zone.cast<DomainA, DomainC>().pointList.size();
math::Matrix ls_mat(nb_nodes, nb_nodes);
bridging_zone.buildLeastSquareMatrix(ls_mat);
ls_mat.factoLU();
ls_mat.inverse();
math::Matrix shapes(nb_atoms, nb_nodes);
bridging_zone.extractShapeMatrix(shapes);
shapes.transpose();
A = new math::Matrix(nb_nodes, nb_atoms);
A->matProduct(ls_mat, shapes);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::buildLeastSquareRHS() {
auto &bridging_zone = *this->bridging_zone;
auto &MDboundary_zone = *this->MDboundary_zone;
rhs->zero();
UInt at_index = 0;
for (auto &&at : bridging_zone.cast<DomainA, DomainC>().pointList) {
(*rhs)(at_index) = at.velocity()[0];
++at_index;
}
A->rightMultiply(*rhs, *CG_vel);
at_index = 0;
for (auto &&at : bridging_zone.cast<DomainA, DomainC>().pointList) {
(*rhs)(at_index) = at.displacement()[0];
++at_index;
}
A->rightMultiply(*rhs, *CG_disp);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::correctContinuum() {
auto &bridging_zone = *this->bridging_zone;
auto &MDboundary_zone = *this->MDboundary_zone;
UInt j = 0;
for (auto &&n : bridging_zone.cast<DomainA, DomainC>().meshList) {
if (this->weightFE[j] <= .5) {
Real speed_projection = this->quality;
if (this->weightFE[j] <= 1e-2)
speed_projection = 1;
Real &vel = n.velocity()[0];
Real &disp = n.displacement()[0];
Real corrected = vel + speed_projection * ((*CG_vel)(j)-vel);
n.velocity()[0] = corrected;
corrected = disp + speed_projection * ((*CG_disp)(j)-disp);
n.displacement()[0] = corrected;
}
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC> void Kobayashi::correctAtoms() {
auto &bridging_zone = *this->bridging_zone;
UInt natom = bridging_zone.getNumberLocalMatchedPoints();
constexpr UInt Dim = DomainA::ContainerPoints::Dim;
ContainerArray<Real> atomic_vel;
ContainerArray<Real> atomic_disp;
atomic_vel.assign(natom * Dim, 0);
atomic_disp.assign(natom * Dim, 0);
bridging_zone.projectAtomicDOFsOnMesh(_displacement, atomic_disp);
bridging_zone.projectAtomicDOFsOnMesh(_velocity, atomic_vel);
UInt at_index = 0;
for (auto &&at : bridging_zone.cast<DomainA, DomainC>().pointList) {
if (this->weightMD[at_index] <= .5) {
Real speed_projection = this->quality;
if (this->weightMD[at_index] <= 1e-2)
speed_projection = 1;
Real &vel = at.velocity()[0];
Real disp = at.displacement()[0];
vel = vel + speed_projection * (atomic_vel[at_index] - vel);
Vector<Dim> corrected;
corrected[0] = disp + this->quality * (atomic_disp[at_index] - disp);
at.displacement() = corrected;
}
++at_index;
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC KOBAYASHI
This keyword is used to set a coupling method
which will implement the method of Kobayashi and co-workers \\ \\
\textit{A simple dynamical scale-coupling method for concurrent simulation of
hybridized atomistic/coarse-grained-particle system,} \\
\textbf{Ryo Kobayashi, Takahide Nakamura, Shuji Ogata,} \\
\textit{International Journal for Numerical Methods in Engineering Volume 83,
Issue 2, pages 249-268, 9 July 2010.} \\ \\
This coupling component works only in 1D and in sequential.
*/
/* LMHERITANCE arlequin_template */
/* LMEXAMPLE
COUPLING_CODE md fe KOBAYASHI GEOMETRY 3 BOUNDARY 4 GRID_DIVISIONX 10 QUALITY
1e-3 VERBOSE
*/
void Kobayashi::declareParams() { ArlequinTemplate::declareParams(); }
/* -------------------------------------------------------------------------- */
// DECLARE_ATOMIC_CONTINUUM_TEMPLATE(Kobayashi)
__END_LIBMULTISCALE__

Event Timeline

c4science · Help