Page MenuHomec4science
Files F72921916

kobayashi.cc
No OneTemporary
Actions

File Metadata

Created
Wed, Jul 17, 15:29

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 "lm_common.hh"
#include "lib_bridging.hh"
#include "kobayashi.hh"
#include "arlequin_template.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
Kobayashi<DomainA,DomainC,Dim>::Kobayashi(const std::string & name,
DomainInterface & A,
DomainInterface & C):
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 = Parent::bridging_zone.pointList.nbElem();
Parent::allocate(taille);
//build weights
Parent::computeContinuumWeights();
Parent::computeAtomWeights();
UInt _taille = Parent::bridging_zone.meshList.nbElem();
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.nbElem();
UInt nb_atoms = Parent::bridging_zone.pointList.nbElem();
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();
typename Parent::IteratorAtomsSubset itrec
= Parent::bridging_zone.pointList.getIterator();
typename Parent::RefAtom at;
UInt at_index=0;
for(at_index=0, at = itrec.getFirst(); !itrec.end() ; ++at_index , at = itrec.getNext()){
(*rhs)(at_index) = at.velocity(0);
}
A->rightMultiply(*rhs,*CG_vel);
for(at_index=0, at = itrec.getFirst(); !itrec.end() ; ++at_index , at = itrec.getNext()) {
(*rhs)(at_index) = at.displacement(0);
}
A->rightMultiply(*rhs,*CG_disp);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void Kobayashi<DomainA,DomainC,Dim>::correctContinuum(){
typename Parent::IteratorNodesSubset itrec = Parent::bridging_zone.meshList.getIterator();
typename Parent::RefNode n = itrec.getFirst();
UInt j = 0;
for (n = itrec.getFirst() ; !itrec.end() ; n = itrec.getNext(), ++j)
if (Parent::weightFE[j] <= .5)
{
Real speed_projection = Parent::quality;
if (Parent::weightFE[j] <= 1e-2) speed_projection = 1;
Real corrected = n.velocity(0) + speed_projection * ((*CG_vel)(j)-n.velocity(0));
n.setVelocity(0,corrected);
corrected = n.displacement(0) + speed_projection * ((*CG_disp)(j)-n.displacement(0));
n.setDisplacement(0,corrected);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void Kobayashi<DomainA,DomainC,Dim>::correctAtoms(){
typename Parent::IteratorAtomsSubset itrec
= Parent::bridging_zone.pointList.getIterator();
typename Parent::RefAtom at;
at = itrec.getFirst();
UInt at_index=0;
UInt natom = Parent::bridging_zone.getLocalPoints();
std::vector<Real> atomic_vel;
std::vector<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);
for (at_index = 0 ; !itrec.end(); ++at_index, at = itrec.getNext())
if (Parent::weightMD[at_index] <= .5)
{
Real speed_projection = Parent::quality;
if (Parent::weightMD[at_index] <= 1e-2) speed_projection = 1;
Real corrected = at.velocity(0) + speed_projection * (atomic_vel[at_index] - at.velocity(0));
at.velocity(0) = corrected;
Real disp = at.position(0) - at.position0(0);
corrected = disp + Parent::quality * (atomic_disp[at_index] - disp);
at.position(0) = at.position0(0) + corrected;
}
}
/* -------------------------------------------------------------------------- */
/* 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__

Event Timeline

c4science · Help