Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

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

Last Change
Raw File

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__
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
Kobayashi<DomainA, DomainC, Dim>::Kobayashi(const std::string &name,
DomainInterface &A,
DomainInterface &C)
: LMObject(name), ArlequinTemplate<DomainA, DomainC, Dim>(name, A, C) {
// naming the bridging zones
std::stringstream fname;
fname << name << "-bridging";
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::init() {
if (lm_world_size > 1)
LM_FATAL("this coupler work only in sequential");
std::stringstream fname;
// forward name
// fname.str("");
// fname << Parent::name << "-bridging";
// Parent::bridging_zone.setName(fname.str());
// initialize the bridging_zone object
Parent::bridging_zone.init();
// allocate the vectors necessary to continue
UInt taille = this->bridging_zone.pointList.size();
Parent::allocate(taille);
// build weights
Parent::computeContinuumWeights();
Parent::computeAtomWeights();
UInt _taille = Parent::bridging_zone.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();
/***************************
treat now MD boundary zone
****************************/
// //forward name
// fname.str("");
// fname << Parent::name << "-boundary";
// Parent::MDboundary_zone.setName(fname.str());
Parent::MDboundary_zone.init();
DUMP("Init ... done", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::coupling(CouplingStage stage) {
if (stage != COUPLING_STEP4)
return;
Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_velocity);
Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_displacement);
buildLeastSquareRHS();
correctContinuum();
correctAtoms();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::buildCGMatrix() {
UInt nb_nodes = Parent::bridging_zone.meshList.size();
UInt nb_atoms = Parent::bridging_zone.pointList.size();
math::Matrix ls_mat(nb_nodes, nb_nodes);
Parent::bridging_zone.buildLeastSquareMatrix(ls_mat);
ls_mat.factoLU();
ls_mat.inverse();
math::Matrix shapes(nb_atoms, nb_nodes);
Parent::bridging_zone.extractShapeMatrix(shapes);
shapes.transpose();
A = new math::Matrix(nb_nodes, nb_atoms);
A->matProduct(ls_mat, shapes);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::buildLeastSquareRHS() {
rhs->zero();
UInt at_index = 0;
for (auto &&at : Parent::bridging_zone.pointList) {
(*rhs)(at_index) = at.velocity()[0];
++at_index;
}
A->rightMultiply(*rhs, *CG_vel);
at_index = 0;
for (auto &&at : Parent::bridging_zone.pointList) {
(*rhs)(at_index) = at.displacement()[0];
++at_index;
}
A->rightMultiply(*rhs, *CG_disp);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::correctContinuum() {
UInt j = 0;
for (auto &&n : Parent::bridging_zone.meshList) {
if (Parent::weightFE[j] <= .5) {
Real speed_projection = Parent::quality;
if (Parent::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.setVelocity(0, corrected);
corrected = disp + speed_projection * ((*CG_disp)(j)-disp);
n.setDisplacement(0, corrected);
}
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::correctAtoms() {
UInt natom = Parent::bridging_zone.getLocalPoints();
ContainerArray<Real> atomic_vel;
ContainerArray<Real> atomic_disp;
atomic_vel.assign(natom * Dim, 0);
atomic_disp.assign(natom * Dim, 0);
Parent::bridging_zone.projectAtomicDOFsOnMesh(_displacement, atomic_disp);
Parent::bridging_zone.projectAtomicDOFsOnMesh(_velocity, atomic_vel);
UInt at_index = 0;
for (auto &&at : Parent::bridging_zone.pointList) {
if (Parent::weightMD[at_index] <= .5) {
Real speed_projection = Parent::quality;
if (Parent::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 + Parent::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
*/
template <typename DomainA, typename DomainC, UInt Dim>
void Kobayashi<DomainA, DomainC, Dim>::declareParams() {
Parent::declareParams();
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(Kobayashi)
__END_LIBMULTISCALE__
c4science · Help