Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/continuum/1d/domain_meca1d.cc
1d1b2315bf45master

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 "domain_meca1d.hh"
#include "ball.hh"
#include "geometry_manager.hh"
#include "lm_common.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::gradient() -> ArrayView {
return ArrayView(mymesh->f.data(), mymesh->f.rows(), 1);
}
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::primal() -> ArrayView {
return ArrayView(mymesh->u.data(), mymesh->u.rows(), 1);
}
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::primalTimeDerivative() -> ArrayView {
return ArrayView(mymesh->v.data(), mymesh->v.rows(), 1);
}
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::acceleration() -> ArrayView {
return ArrayView(mymesh->a.data(), mymesh->a.rows(), 1);
}
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::mass() -> ArrayView {
return ArrayView(mymesh->m.data(), mymesh->m.rows(), 1);
}
/* -------------------------------------------------------------------------- */
Real DomainMeca1D::potentialEnergy() { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
void DomainMeca1D::init() {
if (!this->getCommGroup().amIinGroup())
return;
Real rho_0 = atomic_mass / r0;
auto g = GeometryManager::getManager().getGeometry(geometries[0]);
DUMP("atomic mass " << atomic_mass, DBG_MESSAGE);
DUMP("mass density " << rho_0, DBG_MESSAGE);
DUMP(" Meca1D R0 " << r0 << " RCUT " << rcut << " Rcut/r0 "
<< (Real)(rcut / r0),
DBG_MESSAGE);
mymesh = std::make_shared<MeshMeca1D>();
if (this->mesh_file == "")
mymesh->buildMesh(element_size, *g, rho_0);
else
element_size = mymesh->readMesh(this->mesh_file, rho_0);
mesh_container.getContainerElems().setMesh(mymesh);
mesh_container.getContainerNodes().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));
DUMP("E " << E, DBG_MESSAGE);
DUMP("wave_velocity " << sqrt((E / rho_0) * code_unit_system.e_m2dd_tt),
DBG_MESSAGE);
mymesh->setE(E);
}
/* -------------------------------------------------------------------------- */
auto DomainMeca1D::getField(FieldType field_type) -> ArrayView {
auto create = [](auto &&v, UInt dim) -> ArrayView {
return ArrayView(v.data(), v.rows(), dim);
};
switch (field_type) {
case _position0:
return create(mymesh->p0, 1);
case _velocity:
return create(mymesh->v, 1);
case _displacement:
return create(mymesh->u, 1);
default:
LM_FATAL("Not accessible field: " << field_type);
}
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::enforceCompatibility() {}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::updateGradient() {
STARTTIMER("ElastStepForces");
mymesh->f = 0.;
for (UInt i = 0; i < mymesh->nbElements(); ++i)
this->updateNodeForce(i);
STOPTIMER("ElastStepForces");
// 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->changeRelease();
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::updateAcceleration() {
mymesh->a = code_unit_system.ft_m2v * mymesh->f / mymesh->m;
}
/* -------------------------------------------------------------------------- */
UInt DomainMeca1D::updateNodeForce(UInt i) {
auto &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 (std::isnan(F) || std::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);
mymesh->f(elt.n2) += -P;
mymesh->f(elt.n1) += P;
#ifndef LM_OPTIMIZED
if (std::isnan(P) || std::isinf(P))
LM_FATAL("computing bad stress tensor");
#endif
return 0;
}
/* -------------------------------------------------------------------------- */
// void DomainMeca1D::updateHeatRate(int i) {
// auto &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) { return false; }
/* -------------------------------------------------------------------------- */
Cube DomainMeca1D::getDomainBoundingBox() const {
Cube c;
if (not mymesh)
return c;
c.getXmax()[0] = mymesh->p0(mymesh->nbNodes() - 1);
c.getXmin()[0] = mymesh->p0(0);
return c;
}
/* -------------------------------------------------------------------------- */
Cube DomainMeca1D::getSubDomainBoundingBox() const {
return this->getDomainBoundingBox();
}
/* -------------------------------------------------------------------------- */
void DomainMeca1D::setDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
void DomainMeca1D::setSubDomainBoundingBox(const Cube &) { LM_TOIMPLEMENT; }
/* -------------------------------------------------------------------------- */
/* 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
.. code-block::
Section MultiScale RealUnits
...
MODEL MECA1D fe
...
endSection
Section MECA1D:fe RealUnits
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);
/* L22MKEYWORD CONDUCTIVITY
For thermal transport purpose, this provides the conductivity of the
material.
*/
// this->parseKeyword("CONDUCTIVITY", conductivity, 0.);
/* L22MKEYWORD 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);
/* L22MKEYWORD THERMAL
Activate thermal transport feature
*/
// this->parseTag("THERMAL", thermal_flag, false);
/* LMKEYWORD MASS
Specify the atomic mass. Used to compute the mass density of
the finite element formulation
*/
this->parseKeyword("MASS", atomic_mass);
/* LMKEYWORD MESH_FILE
Specify the file where to load nodal coordinates
*/
this->parseKeyword("MESH_FILE", mesh_file, "");
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help