Homec4science
Diffusion LibMultiScale (features/new_parralelism) Browse

libmultiscale/src/continuum/akantu/domain_akantu.cc
b05d5f2aba78features/new_parralelism

Last Change
Raw File

domain_akantu.cc
View Options

/**
* @file domain_akantu.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Mon Jul 21 10:32:33 2014
*
* @brief This is the model wrapping Akantu
*
* @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/>.
*
*/
#define TIMER
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "domain_akantu.hh"
#include "lib_geometry.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <solid_mechanics_model.hh>
#include <material.hh>
#include <material_elastic.hh>
/* -------------------------------------------------------------------------- */
#include <iomanip>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <UInt Dim>
DomainAkantu<Dim>::~DomainAkantu(){
delete (model ); model = NULL;
delete (position0 ); position0 = NULL;
delete (velocity ); velocity = NULL;
delete (displacement ); displacement = NULL;
delete (acceleration ); acceleration = NULL;
delete (applied_force ); applied_force = NULL;
delete (residual ); residual = NULL;
delete (mass ); mass = NULL;
delete (mesh ); mesh = NULL;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainAkantu<Dim>::DomainAkantu(DomainID ID,CommGroup GID):
DomainContinuum<ContainerElemsAkantu<Dim>,
ContainerNodesAkantu<Dim>,Dim>(ID,GID),
mesh(NULL),
model(NULL),
position0(NULL),
velocity(NULL),
displacement(NULL),
acceleration(NULL),
applied_force(NULL),
residual(NULL),
mass(NULL),
surface_tag(""), boundary_start(0),
boundary_is_traction(false),
boundary_has_been_applied(false)
{
for (UInt i = 0; i < 3; ++i) {
pbc[i] = 0;
scale[i] = 0.;
shift[i] = 0.;
for (UInt j = 0 ; j < Dim ; ++j) {
this->boundary_stress[Dim*i + j] = std::numeric_limits<Real>::max();
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::scaleMesh(){
for (UInt i = 0; i < Dim; ++i) {
if (!scale[i]) continue;
Real * nodes = mesh->getNodes().storage();
UInt nb_nodes = mesh->getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
nodes[n*Dim+i] *= scale[i];
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::shiftMesh(){
for (UInt i = 0; i < Dim; ++i) {
if (!shift[i]) continue;
Real * nodes = mesh->getNodes().storage();
UInt nb_nodes = mesh->getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
nodes[n*Dim+i] += shift[i];
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::initMesh(){
mesh = new akantu::Mesh(Dim, this->getID(), 0);
UInt prank = Communicator::getCommunicator().groupRank(lm_my_proc_id, this->getGroupID());
if(prank == 0) {
if (mesh_filename != ""){
mesh->read(mesh_filename);
#ifndef LM_OPTIMIZED
akantu::Mesh::type_iterator it = mesh->firstType(Dim);
akantu::Mesh::type_iterator end = mesh->lastType(Dim);
UInt ntypes = 0;
for(; it != end; ++it, ++ntypes) {}
LM_ASSERT(ntypes <= 1,"At current time Akantu plugin cannot support to load"
<< " more than one element type for volume elements");
#endif //LM_OPTIMIZED
scaleMesh();
shiftMesh();
// UInt nb_fragments = mesh->createClusters(Dim);
// DUMP("the number of fragments is " << nb_fragments,DBG_MESSAGE);
}
else{
LM_ASSERT(Dim == 1,
"automatic mesh generation only possible in 1D"
<< "a mesh file should have been provided using the"
<< " MESH_FILENAME keyword");
const akantu::ElementType type = akantu::_segment_2;
mesh->addConnectivityType(type);
Ball * b = dynamic_cast<Ball*>(GeometryManager::getManager().getGeometry(this->geometries[0]));
if (b == NULL) LM_FATAL("geometry given is not a valid ball");
//number of element on each side of the zero
UInt nx = (UInt)(nearbyint((b->rMax() - b->rMin())/elem_size));
Real center = b->getCenter(0);
Real rmin = b->rMin();
Real rmax = b->rMax();
Real range = rmax - rmin;
Real e_size = range/nx;
DUMP("element size " << e_size,DBG_MESSAGE);
UInt nb_nodes;
nb_nodes = 2*nx + 1;
if (rmin > 0) ++nb_nodes;
UInt nb_elems;
nb_elems = 2*nx;
akantu::Array<Real> & nodes = const_cast<akantu::Array<Real> &>(mesh->getNodes());
nodes.resize(nb_nodes);
// create nodes from -1 to 0
for (UInt i = 0; i < nx + 1; ++i){
nodes(i) = static_cast<Real>(i)/static_cast<Real>(nx) - 1.0;
nodes(i) = nodes(i)*range - rmin + center;
DUMP("creating node at position " << nodes(i),DBG_DETAIL)
}
// create nodes from 0 to 1
if (rmin == 0){
for (UInt i = 0; i < nx + 1; ++i){
nodes(i+nx) = static_cast<Real>(i)/static_cast<Real>(nx);
nodes(i+nx) = nodes(i+nx)*range + rmin + center;
DUMP("creating node at position " << nodes(i+nx),DBG_DETAIL)
}
}
else
for (UInt i = 0; i < nx + 1; ++i){
nodes(i+nx+1) = static_cast<Real>(i)/static_cast<Real>(nx);
nodes(i+nx+1) = nodes(i+nx+1)*range + rmin + center;
DUMP("creating node at position " << nodes(i+nx+1),DBG_DETAIL)
}
akantu::Array<UInt> & connectivity =
const_cast<akantu::Array<UInt> &>(mesh->getConnectivity(type));
connectivity.resize(nb_elems);
for (UInt i = 0; i < nx; ++i){
connectivity(i,0) = i;
connectivity(i,1) = i+1;
}
if (rmin == 0)
for (UInt i = 0; i < nx; ++i){
connectivity(i+nx,0) = i+nx;
connectivity(i+nx,1) = i+nx+1;
}
else
for (UInt i = 0; i < nx; ++i){
connectivity(i+nx,0) = i+nx+1;
connectivity(i+nx,1) = i+nx+2;
}
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::initMaterial(){
if (this->material_filename == ""){
LM_ASSERT(Dim == 1,
"automatic material property definition only possible in 1D"
<< "a material file should have been provided using the"
<< " MATERIAL_FILENAME keyword");
DUMP("Automatic parameter generation for LJ and lattice provided parameters",DBG_WARNING);
akantu::Material & mat = model->registerNewEmptyMaterial<akantu::MaterialElastic<1u> >("elastic1DMaterial");
UInt ordre = static_cast<UInt>(rcut/r0);
Real coeff = MathTools::ipow(sigma/r0,6);
Real E = 24.0*epsilon*coeff/r0*(-7.0*MathTools::zeta(6,ordre)
+ 26.0*coeff*MathTools::zeta(12,ordre));
mat.setParam("E", E);
mat.setParam("nu", 0.0);
mat.setParam("rho", Real(atomic_mass/r0));
DUMP("mass density " << (Real)(atomic_mass/r0),DBG_MESSAGE);
DUMP("r0 " << r0,DBG_MESSAGE);
DUMP("atomic_mass " << atomic_mass,DBG_MESSAGE);
}
this->model->initMaterials();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::init(){
if (!Communicator::getCommunicator().amIinGroup(this->getGroupID())) return;
akantu::debug::setDebugLevel(akantu::dbl0);
// akantu::debug::setDebugLevel(akantu::dblTrace);
initMesh();
UInt prank = Communicator::getCommunicator().groupRank(lm_my_proc_id, this->getGroupID());
UInt psize = Communicator::getCommunicator().getNBprocsOnGroup(this->getGroupID());
akantu::MeshPartition * partition = NULL;
if(prank == 0) {
partition = new akantu::MeshPartitionScotch(*mesh, Dim);
partition->partitionate(psize);
}
this->model = new akantu::SolidMechanicsModel(*mesh, Dim,this->getID());
this->model->initParallel(partition);
delete partition;
this->model->getMesh().computeBoundingBox();
this->model->setF_M2A(code_unit_system.ft_m2v);
bool pbc_activated = false;
for (UInt i = 0; i < Dim; ++i) {
pbc_activated |= pbc[i];
}
if (pbc_activated) {
this->model->setPBC(pbc[0],pbc[1],pbc[2]);
std::map<UInt,UInt> & akantu_pairs = this->model->getPBCPairs();
std::map<UInt,UInt>::iterator it = akantu_pairs.begin();
std::map<UInt,UInt>::iterator end = akantu_pairs.end();
while(it != end) {
this->pbc_pairs.push_back((*it));
++it;
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::performStep1(){
if (this->prestressed_periodicity_flag && !this->slave_displacements_computed) {
this->buildSlaveDisplacements();
this->slave_displacements_computed = true;
}
if ((this->haveToApplyNeumann()) && !this->boundary_has_been_applied &&
(current_step >= this->boundary_start)) {
this->enableNeumanBndyConditions();
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::performStep3(){
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
UInt DomainAkantu<Dim>::getCurrentStep(){
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::buildSlaveDisplacements() {
typedef std::vector<std::pair<UInt, UInt> >::iterator pair_it;
pair_it pbc_pair = this->pbc_pairs.begin();
pair_it end = this->pbc_pairs.end();
// upper bounds are masters, lower bounds slaves
for (; pbc_pair != end; ++pbc_pair) {
for (UInt direction = 0 ; direction < Dim ; ++direction) {
this->slave_displacements.push_back
((*this->displacement)[Dim * pbc_pair->first + direction]-
(*this->displacement)[Dim * pbc_pair->second + direction]);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::checkBoundaryInputSanity() {
try{
this->mesh->template createGroupsFromMeshData<std::string>("physical_names");
} catch (...) {
this->mesh->template createGroupsFromMeshData<UInt>("tag_0");
}
bool is_stress;
Real realmax = std::numeric_limits<Real>::max();
this->boundary_is_traction = ((this->boundary_stress[0] != realmax) &&
(this->boundary_stress[Dim] == realmax));
is_stress = (this->boundary_stress[Dim] != realmax);
if (this->boundary_is_traction && is_stress) {
LM_FATAL("you cannot specify both a traction (TRACTION_BNDY_VECTOR) and "
<< "a stress (TRACTION_BNDY_TENSOR) for traction "
<< "boundary conditions");
} else if (!this->boundary_is_traction && !is_stress) {
LM_FATAL("you specified a TRACTION_BNDY_TAG, but neither a stress "
<< "(TRACTION_BNDY_TENSOR) nor a traction "
<< "(TRACTION_BNDY_VECTOR)");
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::enableNeumanBndyConditions() {
try {
if (boundary_is_traction) {
akantu::Vector<Real> surface_traction(Dim);
for (UInt i = 0 ; i < Dim ; ++i) {
surface_traction(i) = this->boundary_stress[i];
}
model->applyBC(akantu::BC::Neumann::FromSameDim(surface_traction),
surface_tag);
} else {
akantu::Matrix<Real> surface_stress(Dim, Dim);
for (UInt i = 0 ; i < Dim ; ++i) {
for (UInt j = 0 ; j < Dim ; ++j) {
surface_stress(i,j) = this->boundary_stress[Dim*i + j];
}
}
model->applyBC(akantu::BC::Neumann::FromHigherDim(surface_stress),
surface_tag);
}
} catch (akantu::debug::Exception & ) {}
this->boundary_has_been_applied = true;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::setTimeStep(Real ts){
model->setTimeStep(ts);
this->timeStep = ts;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::setCurrentStep(UInt step){
current_step = step;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
Real DomainAkantu<Dim>::getEpot(){
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::getSubDomainDimensions(Real * xmin,Real * xmax){
akantu::Vector<Real> _xmin;
akantu::Vector<Real> _xmax;
_xmin = mesh->getLocalLowerBounds();
_xmax = mesh->getLocalUpperBounds();
for (UInt k = 0; k < Dim; ++k) {
xmin[k] = _xmin[k];
xmax[k] = _xmax[k];
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantu<Dim>::getDomainDimensions(Real * xmin,Real * xmax){
akantu::Vector<Real> _xmin;
akantu::Vector<Real> _xmax;
_xmin = mesh->getLowerBounds();
_xmax = mesh->getUpperBounds();
for (UInt k = 0; k < Dim; ++k) {
xmin[k] = _xmin[k];
xmax[k] = _xmax[k];
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC AKANTU_BASE
This domain implements the plugin with Akantu Finite Element library.
*/
/* LMHERITANCE domain_continuum */
template <UInt Dim>
void DomainAkantu<Dim>::declareParams(){
DomainContinuum<ContainerElemsAkantu<Dim>,ContainerNodesAkantu<Dim>,Dim>::declareParams();
/* LMKEYWORD MESH_FILENAME
Specify the mesh file to be loaded using the MeshIO utility of Akantu
*/
this->parseKeyword("MESH_FILENAME",mesh_filename,"");
/* LMKEYWORD MATERIAL_FILENAME
Specify the material file to be loaded for the Akantu SolidMechanicsModel object
*/
this->parseKeyword("MATERIAL_FILENAME",material_filename,"");
/* LMKEYWORD PBC
Specify the directions in which Periodic Boundary Conditions should be activated
*/
this->parseVectorKeyword("PBC",Dim,pbc,VEC_DEFAULTS(0u,0u,0u));
/* LMKEYWORD SCALE
Specify factors for each directions to be used in rescaling the mesh nodal positions
*/
this->parseVectorKeyword("SCALE",Dim,scale,VEC_DEFAULTS(1.,1.,1.));
/* LMKEYWORD SHIFT
Specify factors for each directions to be used in shifting the mesh nodal positions
*/
this->parseVectorKeyword("SHIFT",Dim,shift,VEC_DEFAULTS(0.,0.,0.));
/* LMKEYWORD R0
In the case of 1D mesh allows to use a linearized LJ material definition //
fitting monoatomic chain:\\
This provides the interatomic distance.
*/
this->parseKeyword("R0",r0,0.);
/* LMKEYWORD RCUT
In the case of 1D mesh allows to use a linearized LJ material definition
fitting monoatomic chain:\\
This provides the cutoff radius for the potential evaluation.
*/
this->parseKeyword("RCUT",rcut,0.);
/* LMKEYWORD EPSILON
In the case of 1D mesh allows to use a linearized LJ material definition
fitting monoatomic chain:\\
This provides the $\epsilon$ parameter in the LJ potential definition.
*/
this->parseKeyword("EPSILON",epsilon,0.);
/* LMKEYWORD SIGMA
In the case of 1D mesh allows to use a linearized LJ material definition
fitting monoatomic chain:\\
This provides the $\sigma$ parameter in the LJ potential definition.
*/
this->parseKeyword("SIGMA",sigma,0.);
/* LMKEYWORD MASS
In the case of 1D mesh allows to use a linearized LJ material definition
fitting monoatomic chain:\ \
This provides the atomic mass.
*/
this->parseKeyword("MASS",atomic_mass,0.);
/* LMKEYWORD ELEM_SIZE
In the case of 1D mesh generates automatically a homogeneous
1D first order mesh.
*/
this->parseKeyword("ELEM_SIZE",elem_size,0.);
/* LMKEYWORD TRACTION_BNDY_TENSOR
Specify a full stress tensor to be applied to a boundary */
this->parseVectorKeyword("TRACTION_BNDY_TENSOR", Dim*Dim, this->boundary_stress,VEC_DEFAULTS(0.,0.,0.,0.,0.,0.,0.,0.,0.));
// /* L2MKEYWORD TRACTION_BNDY_VECTOR
// Specify a traction vector to be applied to a boundary */
// this->parseVectorKeyword("TRACTION_BNDY_VECTOR", Dim, this->boundary_stress,
// {0.,0.,0.});
/* LMKEYWORD TRACTION_BNDY_TAG
Specify the tag of the boundary on which the traction should be applied.
This refers to "tag\_1" in the msh format */
this->parseKeyword("TRACTION_BNDY_TAG", this->surface_tag,"");
/* LMKEYWORD TRACTION_BNDY_START
Specify the time step at which to start applying the traction boundary
conditions */
this->parseKeyword("TRACTION_BNDY_START", this->boundary_start,0u);
}
/* -------------------------------------------------------------------------- */
template class DomainAkantu<3>;
template class DomainAkantu<2>;
template class DomainAkantu<1>;
__END_LIBMULTISCALE__
c4science · Help