Homec4science
Diffusion LibMultiScale (features/components) Browse

libmultiscale/src/continuum/akantu/domain_akantu_dynamic.cc
b05d5f2aba78features/components

Last Change
Raw File

domain_akantu_dynamic.cc
View Options

/**
* @file domain_akantu_dynamic.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Tue Jul 22 14:47:56 2014
*
* @brief This is the model wrapping Akantu in its dynamic solve version
*
* @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_dynamic.hh"
#include "lib_geometry.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <mesh_io.hh>
#include <solid_mechanics_model.hh>
#ifdef AKANTU_HEAT_TRANSFER
#include <heat_transfer_model.hh>
#endif
#include <material.hh>
#include <material_elastic.hh>
/* -------------------------------------------------------------------------- */
#include <iomanip>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainAkantuDynamic<Dim>::~DomainAkantuDynamic(){
#ifdef AKANTU_HEAT_TRANSFER
delete (this->heat_transfer_model); this->heat_transfer_model = NULL;
delete (this->temperature ); this->temperature = NULL;
delete (this->temperature_rate ); this->temperature_rate = NULL;
delete (this->heat_flux ); this->heat_flux = NULL;
#endif
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainAkantuDynamic<Dim>::DomainAkantuDynamic(DomainID ID,CommGroup GID):
DomainAkantu<Dim>(ID,GID),
heat_transfer_model(NULL),
heat_transfer_flag(false),
temperature(NULL),
temperature_rate(NULL),
heat_flux(NULL),
heat_boundary_geom(invalidGeom)
{}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::init(){
if (this->material_filename != "")
this->akantu_parser.parse(this->material_filename);
this->model->setParser(this->akantu_parser);
this->model->initFull(akantu::SolidMechanicsModelOptions(akantu::_explicit_lumped_mass,true));
this->initMaterial();
akantu::Array<bool> & bound = this->model->getBlockedDOFs();
UInt full_counter = 0;
for (UInt i = 0 ; i < bound.getSize() ; ++i) {
if (bound(i)) {
full_counter ++;
}
}
this->elems.setSolidMechanicsModel(*this->model);
this->nodes.setSolidMechanicsModel(*this->model);
Real time_step_suggested = this->model->getStableTimeStep() / sqrt(code_unit_system.e_m2dd_tt);
DUMP("suggested timestep is for solid mechanics model " << time_step_suggested, DBG_MESSAGE);
if (this->timeStep > time_step_suggested)
LM_FATAL("provided timestep is larger than the critical timestep: please adapt it");
this->model->setTimeStep(this->timeStep);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::initHeatTransfer(){
#ifndef AKANTU_HEAT_TRANSFER
LM_FATAL("Akantu was not compiled using the heat transfer model");
#else
heat_transfer_model = new akantu::HeatTransferModel(*this->mesh,Dim,this->getID()+"heat");
//initialize PBC
bool pbc_activated = false;
for (UInt i = 0; i < Dim; ++i) {
pbc_activated |= this->pbc[i];
}
if (pbc_activated) {
heat_transfer_model->setPBC(this->pbc[0],this->pbc[1],this->pbc[2]);
std::map<UInt,UInt> & akantu_pairs = heat_transfer_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;
}
}
//initialize everything
if (heat_transfer_material_filename != "")
heat_transfer_model->initFull();
else
LM_FATAL("Heat transfer material file not provided.");
this->elems.setHeatTransferModel(*heat_transfer_model);
this->nodes.setHeatTransferModel(*heat_transfer_model);
this->elems.setHeatTransferFlag(heat_transfer_flag);
this->nodes.setHeatTransferFlag(heat_transfer_flag);
//assemble the lumped capacity
heat_transfer_model->assembleCapacityLumped();
//get stable time step
akantu::Real suggested_time_step_heat_transfer = heat_transfer_model->getStableTimeStep()*0.8;
DUMP("suggested timestep for heat transfer model is " << suggested_time_step_heat_transfer,DBG_MESSAGE);
if (this->timeStep > suggested_time_step_heat_transfer)
LM_FATAL("provided timestep is larger than the critical timestep for the heat transfer model: please adapt it");
heat_transfer_model->setTimeStep(this->timeStep);
/* -------------------- */
// initialize the vectors
/* -------------------- */
temperature = new VecAkantu(heat_transfer_model->getTemperature());
heat_flux = new VecAkantu(heat_transfer_model->getExternalHeatRate());
temperature_rate = new VecAkantu(heat_transfer_model->getTemperatureRate());
heat_transfer_model->updateResidual();
akantu::Array<bool> & boundary = heat_transfer_model->getBlockedDOFs();
UInt nb_nodes = boundary.getSize();
bool * bound = boundary.storage();
if (heat_boundary_geom != invalidGeom){
Geometry * geom_dirichlet = GeometryManager::getManager().getGeometry(heat_boundary_geom);
for (UInt n = 0; n < nb_nodes; ++n, ++bound){
Real pos[Dim];
this->nodes.get(n).getPositions0(pos);
if (geom_dirichlet->contains<Dim>(pos)){
*bound = true;
}
}
}
#endif
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::performStep1(){
STARTTIMER("AkantuStepPosition");
DomainAkantu<Dim>::performStep1();
UInt nb_nodes = this->model->getFEEngine().getMesh().getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
for (UInt i = 0; i < Dim; ++i) {
(*this->velocity)[Dim*n+i] +=
code_unit_system.ft_m2v * .5 * (*this->residual)[Dim*n+i] / (*this->mass)[Dim*n+i] * this->timeStep;
LM_ASSERT(!isnan((*this->velocity)[Dim*n+i]) && !isinf((*this->velocity)[Dim*n+i]),
"problem with velocity on node " << n << " Dim " << i
<< " " << (*this->residual)[Dim*n+i] << " " << (*this->mass)[Dim*n+i]);
(*this->displacement)[Dim*n+i] += this->timeStep * (*this->velocity)[Dim*n+i];
}
}
STOPTIMER("AkantuStepPosition");
#ifdef AKANTU_HEAT_TRANSFER
if(heat_transfer_flag)
heat_transfer_model->explicitPred();
#endif
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::performStep2(){
STARTTIMER("AkantuStepForce");
this->model->initializeUpdateResidualData();
if (this->prestressed_periodicity_flag) {
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();
UInt node_counter = 0;
for (; pbc_pair != end; ++pbc_pair, ++node_counter) {
for (UInt direction = 0 ; direction < Dim ; ++direction) {
(*this->displacement)[Dim*pbc_pair->first + direction] +=
this->slave_displacements[Dim * node_counter+direction];
}
}
}
this->model->updateResidual(false);
// model->updateAcceleration();
STOPTIMER("AkantuStepForce");
#ifdef AKANTU_HEAT_TRANSFER
if(heat_transfer_flag)
heat_transfer_model->updateResidual();
#endif
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::performStep3(){
STARTTIMER("AkantuStepVelocity");
DomainAkantu<Dim>::performStep3();
UInt nb_nodes = this->model->getFEEngine().getMesh().getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
for (UInt i = 0; i < Dim; ++i) {
(*this->velocity)[Dim*n+i] +=
code_unit_system.ft_m2v * .5 * (*this->residual)[Dim*n+i] / (*this->mass)[Dim*n+i] * this->timeStep;
LM_ASSERT(!isnan((*this->velocity)[Dim*n+i]) && !isinf((*this->velocity)[Dim*n+i]),
"problem with velocity on node " << n << " Dim " << i
<< " " << (*this->velocity)[Dim*n+i]
<< " " << (*this->residual)[Dim*n+i] << " " << (*this->mass)[Dim*n+i]
<< " " << this->timeStep);
}
}
STOPTIMER("AkantuStepVelocity");
#ifdef AKANTU_HEAT_TRANSFER
if(heat_transfer_flag)
heat_transfer_model->explicitCorr();
#endif
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
/* LMDESC AKANTU_DYNAMIC
This domain implements the plugin with Akantu Finite Element library.
*/
/* LMHERITANCE domain_akantu */
template <UInt Dim>
void DomainAkantuDynamic<Dim>::declareParams(){
//DomainAkantu<Dim>::declareParams();
/* LMKEYWORD HEAT_TRANSFER
To switch on the heat transfer model object
*/
this->parseTag("HEAT_TRANSFER",this->heat_transfer_flag,false);
/* LMKEYWORD HEAT_TRANSFER_MATERIAL_FILENAME
Specify the material file to be loaded for the Akantu HeatTransferModel object
*/
this->parseKeyword("HEAT_TRANSFER_MATERIAL_FILENAME",heat_transfer_material_filename,"");
/* LMKEYWORD HEAT_BOUNDARY
Specify the geometry for which a Dirichlet boundary condition should be used for
the HeatTransfer model
*/
this->parseKeyword("HEAT_BOUNDARY",heat_boundary_geom,invalidGeom);
}
template class DomainAkantuDynamic<1>;
template class DomainAkantuDynamic<2>;
template class DomainAkantuDynamic<3>;
__END_LIBMULTISCALE__
c4science · Help