Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/continuum/1d/domain_meca1d.cc
b05d5f2aba78features/coupler_compon…

Last Change
Raw File

domain_meca1d.cc
View Options

/**
* @file domain_meca1d.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Tue Jan 21 12:38:14 2014
*
* @brief Domain for the 1D finite element engine of LibMultiScale
*
* @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.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "domain_meca1d.hh"
#include "math_tools.hh"
#include "geometry_manager.hh"
#include "ball.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
void DomainMeca1D::init(){
if (!Communicator::getCommunicator().amIinGroup(this->getGroupID()))
return;
Real rho_0 = mass/r0;
Geometry * g = GeometryManager::getManager().getGeometry(geometries[0]);
DUMP("atomic mass " << mass,DBG_MESSAGE);
DUMP("mass density " << rho_0,DBG_MESSAGE);
DUMP("r0 " << r0,DBG_MESSAGE);
mymesh = new MeshMeca1D();
if (this->mesh_file == "")
mymesh->buildMesh(element_size,*g,rho_0);
else
element_size = mymesh->readMesh(this->mesh_file,rho_0);
elems.setMesh(mymesh);
nodes.setMesh(mymesh);
DUMP("mass density " << rho_0 << " element size " << element_size , DBG_MESSAGE);
UInt order = static_cast<UInt>(rcut/r0);
Real coeff = MathTools::ipow(sigma/r0,6);
Real E = 24.0*epsilon*coeff/r0*(-7.0*MathTools::zeta(6,order)
+ 26.0*coeff*MathTools::zeta(12,order));
mymesh->setE(E);
this->mesh_container.setRelease(INITIAL_MODEL_RELEASE);
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::performStep1(){
STARTTIMER("ElastStepPositions");
if (current_step % multiple_timestep_freq == 0){
for (UInt i=0;i< mymesh->nbNodes();++i)
{
mymesh->v[i] += code_unit_system.ft_m2v*.5*mymesh->f[i]/mymesh->m[i]*timeStep;
mymesh->u[i] += timeStep*mymesh->v[i];
DUMP("POSstep : " << i << " " << timeStep << " " << mymesh->m[i] << " "
<< mymesh->f[i] << " " << mymesh->u[i] << " " << mymesh->v[i],DBG_DETAIL);
#ifndef LM_OPTIMIZED
if (isnan(mymesh->u[i]) || isinf(mymesh->u[i]))
LM_FATAL("computing bad displacement field for node " << i);
#endif
// mymesh->f_old[i]=mymesh->f[i];
mymesh->f[i]=0.0;
}
}
STOPTIMER("ElastStepPositions");
if (thermal_flag && current_step % multiple_timestep_freq == 0){
for (UInt i=0;i< mymesh->nbNodes();++i){
mymesh->T[i] += mymesh->Tdt[i]*timeStep;
}
}
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::performStep3(){
STARTTIMER("ElastStepVelocities");
if (current_step % multiple_timestep_freq == 0){
for (UInt i=0;i< mymesh->nbNodes();++i){
mymesh->v[i] += code_unit_system.ft_m2v*.5/mymesh->m[i]*mymesh->f[i]*timeStep;
DUMP("VELstep : " << timeStep << " " << mymesh->m[i] << " "
<< mymesh->f[i] << " " << mymesh->v[i],DBG_DETAIL);
}
}
STOPTIMER("ElastStepVelocities");
Real density = mass/r0;
if (thermal_flag && current_step % multiple_timestep_freq == 0){
for (unsigned int i=0;i< mymesh->nbNodes();++i){
mymesh->Tdt[i] = mymesh->heat_rate[i]/density/capacity;
}
}
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::performStep2(){
STARTTIMER("ElastStepForces");
if (current_step % multiple_timestep_freq == 0){
for (UInt i=0;i< mymesh->nbElements();++i)
UpdateNodeForce(i);
for (UInt i=0;i< mymesh->nbNodes();++i)
DUMP("force pour " <<i << " :" << mymesh->f[i],DBG_DETAIL);
}
STOPTIMER("ElastStepForces");
if (thermal_flag && current_step % multiple_timestep_freq == 0){
for (UInt i = 0; i < mymesh->nbNodes(); ++i) {
mymesh->heat_rate[i] = mymesh->external_heat_rate[i];
}
for (unsigned int i=0;i< mymesh->nbElements();++i){
updateHeatRate(i);
}
}
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
UInt DomainMeca1D::UpdateNodeForce(UInt i)
{
MeshMeca1D::FEMeca1D * elt = &mymesh->elt[i];
UInt ordre;
Real F = (mymesh->u[elt->n2] - mymesh->u[elt->n1])/
fabs(mymesh->p0[elt->n2] - mymesh->p0[elt->n1]);
#ifndef LM_OPTIMIZED
if (isnan(F) || isinf(F)) LM_FATAL("computing bad deformation tensor");
#endif
F+= 1;
if (F < 0) LM_FATAL("Deformation gradient F cannot be negative");
LM_ASSERT(F>0, "Deformation gradient F cannot be negative");
ordre = (int)(rcut/F/r0);
DUMP("F(" << i << ") = " << F,DBG_DETAIL);
Real coeff = MathTools::ipow(sigma/F/r0,6);
Real z6 = MathTools::zeta(6,ordre);
Real z12 = MathTools::zeta(12,ordre);
Real P = 24*epsilon*coeff/F/r0*(z6 - 2.*coeff*z12);
if(elastic_flag){
coeff = MathTools::ipow(sigma/r0,6);
P = 24.0*epsilon*coeff/r0*(z6 - 2.0*coeff*z12) +
24.0*epsilon*coeff/r0*(-7.0*z6 + 26.0*coeff*z12)*(F-1.0);
}
DUMP("resulting P = " << P,DBG_DETAIL);
if (arlequin_flag){
mymesh->f[elt->n2]+= -1.0*P/2*(mymesh->alpha[elt->n2] + mymesh->alpha[elt->n1]);
mymesh->f[elt->n1]+= P/2*(mymesh->alpha[elt->n2] + mymesh->alpha[elt->n1]);
}
else{
mymesh->f[elt->n2]+= -P;
mymesh->f[elt->n1]+= P;
}
#ifndef LM_OPTIMIZED
if (isnan(P) || isinf(P)) LM_FATAL("computing bad stress tensor");
#endif
return 0;
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::updateHeatRate(int i) {
MeshMeca1D::FEMeca1D * elt = &mymesh->elt[i];
Real TemperatureGradient = (mymesh->T[elt->n2] - mymesh->T[elt->n1])/element_size;
Real Ni = 1./element_size;
mymesh->heat_rate[elt->n2] += -1.0 * conductivity * TemperatureGradient* Ni;
mymesh->heat_rate[elt->n1] += conductivity * TemperatureGradient * Ni;
}
/* -------------------------------------------------------------------------- */
bool DomainMeca1D::isPeriodic(UInt dir){
return false;
}
/* -------------------------------------------------------------------------- */
/* LMDESC MECA1D
This section describe the section that is associated
with the 1D continuum mechanics model Meca1D which maipulates
homogeneous 1D mesh with P1 elements. The constitutive laws
available are atomic crystalline based such as the Cauchy-Born
or the linearized Cauchy Born based on a Lennard Jones potential.
Thus Material parameters requested are Lennard Jones ones.
This is a plugin
not using an external library all coded within LibMultiScale.
*/
/* LMEXAMPLE
Section MultiScale RealUnits\\
... \\
MODEL MECA1D fe \\
... \\
endSection\\\\
Section MECA1D:fe AtomsUnits \\
ELEM_SIZE ELSIZE \\
DOMAIN_GEOMETRY 2 \\
MASS 39.95 \\
R0 r0 \\
RCUT rcut \\
EPSILON 1.6567944e-21/6.94769e-21 \\
SIGMA 1.1 \\
TIMESTEP 10 \\
endSection
*/
/* LMHERITANCE domain_continuum */
inline void DomainMeca1D::declareParams(){
DomainContinuum<ContainerElemsMeca1D,ContainerNodesMeca1D,1u>::declareParams();
/* LMKEYWORD ELEM_SIZE
Specify the element size. Indeed the generated mesh is homogeneously
made of P1 elements. The size is thus required.
*/
this->parseKeyword("ELEM_SIZE",element_size,-1.);
/* LMKEYWORD R0
This provides the inter-atomic spacing.
*/
this->parseKeyword("R0",r0);
/* LMKEYWORD RCUT
This provides the cutoff radius to be used for the Cauchy-Born
constitutive law.
*/
this->parseKeyword("RCUT",rcut);
/* LMKEYWORD CONDUCTIVITY
For thermal transport purpose, this provides the conductivity of the
material.
*/
this->parseKeyword("CONDUCTIVITY",conductivity,0.);
/* LMKEYWORD CAPACITY
For thermal transport purpose, this provides the capacity of the
material.
*/
this->parseKeyword("CAPACITY",capacity,0.);
/* LMKEYWORD EPSILON
This is the Lennard Jones energy parameter
*/
this->parseKeyword("EPSILON",epsilon);
/* LMKEYWORD SIGMA
This is the Lennard Jones distance parameter
*/
this->parseKeyword("SIGMA",sigma);
/* LMKEYWORD ARLEQUIN
Obsolete. Used to activate the weighting of the
DOF when the Arlequin energy weighting is used in its
full form.
*/
this->parseTag("ARLEQUIN",arlequin_flag,false);
/* LMKEYWORD ELASTIC
Activate the linearized lennard jones
and cauchy-born constitutive law.
*/
this->parseTag("ELASTIC",elastic_flag,false);
/* LMKEYWORD THERMAL
Activate thermal transport feature
*/
this->parseTag("THERMAL",thermal_flag,false);
/* LMKEYWORD MULTIPLE_TIMESTEP
Specify the frequency at which the
integration in time are truly performed.
*/
this->parseKeyword("MULTIPLE_TIMESTEP",multiple_timestep_freq,1u);
/* LMKEYWORD MASS
Specify the atomic mass. Used to compute the mass density of
the finite element formulation
*/
this->parseKeyword("MASS",mass);
/* LMKEYWORD MESH_FILE
Specify the file where to load nodal coordinates
*/
this->parseKeyword("MESH_FILE",mesh_file,"");
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help