Homec4science
Diffusion akantu (features/contact) Browse

rAKA/src/model/solid_mechanics/solid_mechanics_model.cc
08c32f4fe99dfeatures/contact

Last Change
Raw File

solid_mechanics_model.cc
View Options

/**
* @file solid_mechanics_model.cc
*
* @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Jan 19 2016
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Akantu 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.
*
* Akantu 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 Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "integration_scheme_2nd_order.hh"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include "element_group.hh"
#include <cmath>
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_PETSC
#include "solver_petsc.hh"
#include "petsc_matrix.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
#include "dumper_field.hh"
#include "dumper_paraview.hh"
#include "dumper_homogenizing_field.hh"
#include "dumper_internal_material_field.hh"
#include "dumper_elemental_field.hh"
#include "dumper_material_padders.hh"
#include "dumper_element_partition.hh"
#include "dumper_iohelper.hh"
#endif
#ifdef AKANTU_DAMAGE_NON_LOCAL
#include "non_local_manager.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const SolidMechanicsModelOptions
default_solid_mechanics_model_options(_explicit_lumped_mass, false);
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
Model(mesh, dim, id, memory_id),
BoundaryCondition<SolidMechanicsModel>(),
time_step(NAN), f_m2a(1.0),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
stiffness_matrix(NULL),
jacobian_matrix(NULL),
material_index("material index", id, memory_id),
material_local_numbering("material local numbering", id, memory_id),
material_selector(new DefaultMaterialSelector(material_index)),
is_default_material_selector(true),
integrator(NULL),
increment_flag(false), solver(NULL),
synch_parallel(NULL),
are_materials_instantiated(false),
non_local_manager(NULL),
pbc_synch(NULL) {
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
spatial_dimension);
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->blocked_dofs = NULL;
this->increment = NULL;
this->increment_acceleration = NULL;
this->dof_synchronizer = NULL;
this->previous_displacement = NULL;
materials.clear();
mesh.registerEventHandler(*this);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete integrator;
delete solver;
delete mass_matrix;
delete velocity_damping_matrix;
if (stiffness_matrix && stiffness_matrix != jacobian_matrix)
delete stiffness_matrix;
delete jacobian_matrix;
delete synch_parallel;
if (is_default_material_selector) {
delete material_selector;
material_selector = NULL;
}
flatten_internal_map::iterator fl_it = this->registered_internals.begin();
flatten_internal_map::iterator fl_end = this->registered_internals.end();
for (; fl_it != fl_end; ++fl_it) {
delete fl_it->second;
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
delete non_local_manager;
non_local_manager = NULL;
#endif
delete pbc_synch;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void SolidMechanicsModel::initFull(const ModelOptions & options) {
Model::initFull(options);
const SolidMechanicsModelOptions & smm_options =
dynamic_cast<const SolidMechanicsModelOptions &>(options);
this->method = smm_options.analysis_method;
// initialize the vectors
this->initArrays();
// set the initial condition to 0
this->force->clear();
this->velocity->clear();
this->acceleration->clear();
this->displacement->clear();
// initialize pbc
if (this->pbc_pair.size() != 0)
this->initPBC();
// initialize the time integration schemes
switch (this->method) {
case _explicit_lumped_mass:
this->initExplicit();
break;
case _explicit_consistent_mass:
this->initSolver();
this->initExplicit();
break;
case _implicit_dynamic:
this->initImplicit(true);
break;
case _static:
this->initImplicit(false);
this->initArraysPreviousDisplacment();
break;
default:
AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel");
break;
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/// create the non-local manager object for non-local damage computations
std::stringstream nl_manager_name;
nl_manager_name << "NLManager" << this->id;
this->non_local_manager = new NonLocalManager(*this, nl_manager_name.str(), this->memory_id);
#endif
// initialize the materials
if (this->parser->getLastParsedFile() != "") {
this->instantiateMaterials();
}
if (!smm_options.no_init_materials) {
this->initMaterials();
}
if (this->increment_flag)
this->initBC(*this, *this->displacement, *this->increment, *this->force);
else
this->initBC(*this, *this->displacement, *this->force);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initParallel(MeshPartition * partition,
DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL)
data_accessor = this;
synch_parallel = &createParallelSynch(partition, data_accessor);
synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary);
synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initFEEngineBoundary() {
FEEngine & fem_boundary = getFEEngineBoundary();
fem_boundary.initShapeFunctions(_not_ghost);
fem_boundary.initShapeFunctions(_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) {
AKANTU_DEBUG_IN();
// in case of switch from implicit to explicit
if (!this->isExplicit())
method = analysis_method;
if (integrator)
delete integrator;
integrator = new CentralDifference();
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
std::stringstream sstr;
sstr << id << ":increment_acceleration";
increment_acceleration =
&(alloc<Real>(sstr.str(), nb_nodes, nb_degree_of_freedom, Real()));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initArraysPreviousDisplacment() {
AKANTU_DEBUG_IN();
this->setIncrementFlagOn();
if (not this->previous_displacement) {
UInt nb_nodes = this->mesh.getNbNodes();
std::stringstream sstr_disp_t;
sstr_disp_t << this->id << ":previous_displacement";
this->previous_displacement = &(this->alloc<Real>(
sstr_disp_t.str(), nb_nodes, this->spatial_dimension, 0.));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*
*/
void SolidMechanicsModel::initArrays() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp;
sstr_disp << id << ":displacement";
// std::stringstream sstr_mass; sstr_mass << id << ":mass";
std::stringstream sstr_velo;
sstr_velo << id << ":velocity";
std::stringstream sstr_acce;
sstr_acce << id << ":acceleration";
std::stringstream sstr_forc;
sstr_forc << id << ":force";
std::stringstream sstr_resi;
sstr_resi << id << ":residual";
std::stringstream sstr_boun;
sstr_boun << id << ":blocked_dofs";
displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
// mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension,
// 0));
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension,
REAL_INIT_VALUE));
blocked_dofs =
&(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
std::stringstream sstr_curp;
sstr_curp << id << ":current_position";
current_position =
&(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType)g;
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end =
mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for (; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it, gt);
material_index.alloc(nb_element, 1, *it, gt);
material_local_numbering.alloc(nb_element, 1, *it, gt);
}
}
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initPBC() {
Model::initPBC();
registerPBCSynchronizer();
// as long as there are ones on the diagonal of the matrix, we can put
// boudandary true for slaves
std::map<UInt, UInt>::iterator it = pbc_pair.begin();
std::map<UInt, UInt>::iterator end = pbc_pair.end();
UInt dim = mesh.getSpatialDimension();
while (it != end) {
for (UInt i = 0; i < dim; ++i)
(*blocked_dofs)((*it).first, i) = true;
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::registerPBCSynchronizer() {
pbc_synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_uv);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass);
synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res);
synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump);
changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
current_position->resize(nb_nodes);
Real * current_position_val = current_position->storage();
Real * position_val = mesh.getNodes().storage();
Real * displacement_val = displacement->storage();
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val,
nb_nodes * spatial_dimension * sizeof(Real));
for (UInt n = 0; n < nb_nodes * spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeUpdateResidualData() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
residual->resize(nb_nodes);
/// copy the forces in residual for boundary conditions
memcpy(residual->storage(), force->storage(),
nb_nodes * spatial_dimension * sizeof(Real));
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
updateCurrentPosition();
AKANTU_DEBUG_OUT();
}
/*----------------------------------------------------------------------------*/
void SolidMechanicsModel::reInitialize() {}
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* This function compute the second member of the motion equation. That is to
* say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given
* by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$
* F_{int} = \int_{\Omega} N \sigma d\Omega@f$
*
*/
void SolidMechanicsModel::updateResidual(bool need_initialize) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
// f = f_ext - f_int
// f = f_ext
if (need_initialize)
initializeUpdateResidualData();
AKANTU_DEBUG_INFO("Compute local stresses");
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
this->non_local_manager->computeAllNonLocalStresses();
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the stress
AKANTU_DEBUG_INFO("Send data for residual assembly");
synch_registry->asynchronousSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant stresses");
// finalize communications
synch_registry->waitEndSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeStresses() {
if (isExplicit()) {
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
// compute stresses on all local elements for each materials
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* Computation of the non local part */
this->non_local_manager->computeAllNonLocalStresses();
#endif
} else {
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStressesFromTangentModuli(_not_ghost);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateResidualInternal() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Update the residual");
// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;
if (method != _static) {
// f -= Ma
if (mass_matrix) {
// if full mass_matrix
Array<Real> * Ma = new Array<Real>(*acceleration, true, "Ma");
*Ma *= *mass_matrix;
/// \todo check unit conversion for implicit dynamics
// *Ma /= f_m2a
*residual -= *Ma;
delete Ma;
} else if (mass) {
// else lumped mass
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->storage();
Real * accel_val = acceleration->storage();
Real * res_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if (!(*blocked_dofs_val)) {
*res_val -= *accel_val ** mass_val / f_m2a;
}
blocked_dofs_val++;
res_val++;
mass_val++;
accel_val++;
}
} else {
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
}
// f -= Cv
if (velocity_damping_matrix) {
Array<Real> * Cv = new Array<Real>(*velocity);
*Cv *= *velocity_damping_matrix;
*residual -= *Cv;
delete Cv;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateAcceleration() {
AKANTU_DEBUG_IN();
updateResidualInternal();
if (method == _explicit_lumped_mass) {
/* residual = residual_{n+1} - M * acceleration_n therefore
solution = increment acceleration not acceleration */
solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs,
f_m2a);
} else if (method == _explicit_consistent_mass) {
solve<NewmarkBeta::_acceleration_corrector>(*increment_acceleration);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveLumped(Array<Real> & x, const Array<Real> & A,
const Array<Real> & b,
const Array<bool> & blocked_dofs,
Real alpha) {
Real * A_val = A.storage();
Real * b_val = b.storage();
Real * x_val = x.storage();
bool * blocked_dofs_val = blocked_dofs.storage();
UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent();
for (UInt n = 0; n < nb_degrees_of_freedom; ++n) {
if (!(*blocked_dofs_val)) {
*x_val = alpha *(*b_val / *A_val);
}
x_val++;
A_val++;
b_val++;
blocked_dofs_val++;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitPred() {
AKANTU_DEBUG_IN();
if (increment_flag) {
if (previous_displacement)
increment->copy(*previous_displacement);
else
increment->copy(*displacement);
}
AKANTU_DEBUG_ASSERT(integrator, "itegrator should have been allocated: "
<< "have called initExplicit ? "
<< "or initImplicit ?");
integrator->integrationSchemePred(time_step, *displacement, *velocity,
*acceleration, *blocked_dofs);
if (increment_flag) {
Real * inc_val = increment->storage();
Real * dis_val = displacement->storage();
UInt nb_degree_of_freedom =
displacement->getSize() * displacement->getNbComponent();
for (UInt n = 0; n < nb_degree_of_freedom; ++n) {
*inc_val = *dis_val - *inc_val;
inc_val++;
dis_val++;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity,
*acceleration, *blocked_dofs,
*increment_acceleration);
if (previous_displacement)
previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStep() {
AKANTU_DEBUG_IN();
EventManager::sendEvent(
SolidMechanicsModelEvent::BeforeSolveStepEvent(method));
this->explicitPred();
this->updateResidual();
this->updateAcceleration();
this->explicitCorr();
EventManager::sendEvent(
SolidMechanicsModelEvent::AfterSolveStepEvent(method));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* Initialize the solver and create the sparse matrices needed.
*
*/
void SolidMechanicsModel::initSolver(__attribute__((unused))
SolverOptions & options) {
#if !defined(AKANTU_USE_MUMPS) && \
!defined(AKANTU_USE_PETSC) // or other solver in the future \todo add
// AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#else
UInt nb_global_nodes = mesh.getNbGlobalNodes();
delete jacobian_matrix;
std::stringstream sstr;
sstr << id << ":jacobian_matrix";
#ifdef AKANTU_USE_PETSC
jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension,
_symmetric, sstr.str(), memory_id);
#else
jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension,
_symmetric, sstr.str(), memory_id);
#endif // AKANTU_USE PETSC
jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
if (!isExplicit()) {
delete stiffness_matrix;
std::stringstream sstr_sti;
sstr_sti << id << ":stiffness_matrix";
#ifdef AKANTU_USE_PETSC
stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension,
_symmetric, sstr.str(), memory_id);
stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
#else
stiffness_matrix =
new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
#endif // AKANTU_USE_PETSC
}
delete solver;
std::stringstream sstr_solv;
sstr_solv << id << ":solver";
#ifdef AKANTU_USE_PETSC
solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str(), memory_id);
#elif defined(AKANTU_USE_MUMPS)
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str(), memory_id);
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif // AKANTU_USE_MUMPS
if (solver)
solver->initialize(options);
#endif // AKANTU_HAS_SOLVER
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initJacobianMatrix() {
#if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)
// @todo make it more flexible: this is an ugly patch to treat the case of non
// fix profile of the K matrix
delete jacobian_matrix;
std::stringstream sstr_sti;
sstr_sti << id << ":jacobian_matrix";
jacobian_matrix =
new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id);
std::stringstream sstr_solv;
sstr_solv << id << ":solver";
delete solver;
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
if (solver)
solver->initialize(_solver_no_options);
#else
AKANTU_DEBUG_ERROR("You need to activate the solver mumps.");
#endif
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
void SolidMechanicsModel::initImplicit(bool dynamic,
SolverOptions & solver_options) {
AKANTU_DEBUG_IN();
method = dynamic ? _implicit_dynamic : _static;
if (!increment)
setIncrementFlagOn();
initSolver(solver_options);
if (method == _implicit_dynamic) {
if (integrator)
delete integrator;
integrator = new TrapezoidalRule2();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initialAcceleration() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ma = f");
Solver * acc_solver = NULL;
std::stringstream sstr;
sstr << id << ":tmp_mass_matrix";
SparseMatrix * tmp_mass =
new SparseMatrix(*mass_matrix, sstr.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solver;
sstr << id << ":solver_mass_matrix";
acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif // AKANTU_USE_MUMPS
acc_solver->initialize();
tmp_mass->applyBoundary(*blocked_dofs);
acc_solver->setRHS(*residual);
acc_solver->solve(*acceleration);
delete acc_solver;
delete tmp_mass;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
stiffness_matrix->clear();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->assembleStiffnessMatrix(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() {
if (!velocity_damping_matrix)
velocity_damping_matrix = new SparseMatrix(
*jacobian_matrix, id + ":velocity_damping_matrix", memory_id);
return *velocity_damping_matrix;
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error,
StaticCommunicator & comm){
AKANTU_DEBUG_IN();
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
error = 0;
Real norm[2] = {0., 0.};
Real * increment_val = increment->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
Real * displacement_val = displacement->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if (!(*blocked_dofs_val) && is_local_node) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
blocked_dofs_val++;
increment_val++;
displacement_val++;
}
}
comm.allReduce(norm, 2, _so_sum);
norm[0] = sqrt(norm[0]);
norm[1] = sqrt(norm[1]);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]),
"Something goes wrong in the solve phase");
if (norm[1] < Math::getTolerance()) {
error = norm[0];
AKANTU_DEBUG_OUT();
// cout<<"Error 1: "<<error<<endl;
return error < tolerance;
}
AKANTU_DEBUG_OUT();
if (norm[1] > Math::getTolerance())
error = norm[0] / norm[1];
else
error = norm[0]; // In case the total displacement is zero!
// cout<<"Error 2: "<<error<<endl;
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm, StaticCommunicator & comm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
norm = 0;
Real * residual_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if (is_local_node) {
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if (!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val;
}
blocked_dofs_val++;
residual_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
}
}
comm.allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm),
"Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance,
Real & norm,
StaticCommunicator & comm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->storage();
Real * mass_val = this->mass->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if (is_local_node) {
for (UInt d = 0; d < spatial_dimension; ++d) {
if (!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val / (*mass_val * *mass_val);
}
blocked_dofs_val++;
residual_val++;
mass_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
mass_val += spatial_dimension;
}
}
comm.allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm),
"Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) {
AKANTU_DEBUG_IN();
Real error = 0;
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance,
Real & error) {
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) {
AKANTU_DEBUG_IN();
Real error = 0;
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance,
Real & error) {
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitPred() {
AKANTU_DEBUG_IN();
if (method == _implicit_dynamic)
integrator->integrationSchemePred(time_step, *displacement, *velocity,
*acceleration, *blocked_dofs);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
if (method == _implicit_dynamic) {
integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity,
*acceleration, *blocked_dofs,
*increment);
} else {
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
bool * boun_val = blocked_dofs->storage();
for (UInt j = 0; j < nb_degree_of_freedom;
++j, ++disp_val, ++incr_val, ++boun_val) {
*incr_val *= (1. - *boun_val);
*disp_val += *incr_val;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateIncrement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(
previous_displacement,
"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
Real * prev_disp_val = previous_displacement->storage();
for (UInt j = 0; j < nb_degree_of_freedom;
++j, ++disp_val, ++incr_val, ++prev_disp_val)
*incr_val = (*disp_val - *prev_disp_val);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updatePreviousDisplacement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(
previous_displacement,
"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,
"Synchronizer registry was not initialized."
<< " Did you call initParallel?");
synch_registry->synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeResidual() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,
"Synchronizer registry was not initialized."
<< " Did you call initPBC?");
synch_registry->synchronize(_gst_smm_res);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if (!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc;
sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
elem.kind = _ek_regular;
Mesh::type_iterator it =
mesh.firstType(spatial_dimension, ghost_type, _ek_regular);
Mesh::type_iterator end =
mesh.lastType(spatial_dimension, ghost_type, _ek_regular);
for (; it != end; ++it) {
elem.type = *it;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
UInt nb_element = mesh.getNbElement(*it);
Array<UInt>::const_scalar_iterator mat_indexes =
material_index(*it, ghost_type).begin();
Array<UInt>::const_scalar_iterator mat_loc_num =
material_local_numbering(*it, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element * spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position, X, *it,
_not_ghost);
Array<Real>::matrix_iterator X_el =
X.begin(spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element;
++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, *it);
Real el_c = mat_val[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
if (!mass && !mass_matrix)
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
Real * vel_val = velocity->storage();
Real * mass_val = mass->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node)
mv2 += *vel_val * *vel_val * *mass_val;
vel_val++;
mass_val++;
}
ekin += mv2;
}
StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type,
UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
spatial_dimension, type,
_not_ghost, filter_element);
Array<Real>::vector_iterator vit = vel_on_quad.begin(spatial_dimension);
Array<Real>::vector_iterator vend = vel_on_quad.end(spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5 * getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
Real * incr_or_velo = NULL;
if (this->method == _static) {
incr_or_velo = this->increment->storage();
} else
incr_or_velo = this->velocity->storage();
Real * forc = this->force->storage();
Real * resi = this->residual->storage();
bool * boun = this->blocked_dofs->storage();
Real work = 0.;
UInt nb_nodes = this->mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = this->mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < this->spatial_dimension; ++i) {
if (count_node) {
if (*boun)
work -= *resi * *incr_or_velo;
else
work += *forc * *incr_or_velo;
}
++incr_or_velo;
++forc;
++resi;
++boun;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum);
if (this->method != _static)
work *= this->time_step;
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
} else if (energy_id == "external work") {
return getExternalWork();
}
Real energy = 0.;
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getEnergy(energy_id);
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
std::vector<Material *>::iterator mat_it;
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy =
this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
for (UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType)g;
Mesh::type_iterator it =
this->mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end =
this->mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for (; it != end; ++it) {
UInt nb_element = this->mesh.getNbElement(*it, gt);
if (!material_index.exists(*it, gt)) {
this->material_index.alloc(nb_element, 1, *it, gt);
this->material_local_numbering.alloc(nb_element, 1, *it, gt);
} else {
this->material_index(*it, gt).resize(nb_element);
this->material_local_numbering(*it, gt).resize(nb_element);
}
}
}
Array<Element>::const_iterator<Element> it = element_list.begin();
Array<Element>::const_iterator<Element> end = element_list.end();
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(), this->getMemoryID());
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
if (!filter.exists(elem.type, elem.ghost_type))
filter.alloc(0, 1, elem.type, elem.ghost_type);
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
this->assignMaterialToElements(&filter);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsAdded(element_list, event);
}
// if (method == _explicit_lumped_mass)
// this->assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(
__attribute__((unused)) const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
__attribute__((unused))
const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if (displacement)
displacement->resize(nb_nodes);
if (mass)
mass->resize(nb_nodes);
if (velocity)
velocity->resize(nb_nodes);
if (acceleration)
acceleration->resize(nb_nodes);
if (force)
force->resize(nb_nodes);
if (residual)
residual->resize(nb_nodes);
if (blocked_dofs)
blocked_dofs->resize(nb_nodes);
if (previous_displacement)
previous_displacement->resize(nb_nodes);
if (increment_acceleration)
increment_acceleration->resize(nb_nodes);
if (increment)
increment->resize(nb_nodes);
if (current_position)
current_position->resize(nb_nodes);
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onNodesAdded(nodes_list, event);
}
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
if (method != _explicit_lumped_mass) {
this->initSolver();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(__attribute__((unused))
const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
__attribute__((unused))
const RemovedNodesEvent & event) {
if (displacement)
mesh.removeNodesFromArray(*displacement, new_numbering);
if (mass)
mesh.removeNodesFromArray(*mass, new_numbering);
if (velocity)
mesh.removeNodesFromArray(*velocity, new_numbering);
if (acceleration)
mesh.removeNodesFromArray(*acceleration, new_numbering);
if (force)
mesh.removeNodesFromArray(*force, new_numbering);
if (residual)
mesh.removeNodesFromArray(*residual, new_numbering);
if (blocked_dofs)
mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
if (increment_acceleration)
mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if (increment)
mesh.removeNodesFromArray(*increment, new_numbering);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
if (method != _explicit_lumped_mass) {
this->initSolver();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::reinitializeSolver() {
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
if (method != _explicit_lumped_mass) {
this->initSolver();
}
if (method == _explicit_lumped_mass)
this->assembleMassLumped();
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::isInternal(const std::string & field_name,
const ElementKind & element_kind) {
bool is_internal = false;
/// check if at least one material contains field_id as an internal
for (UInt m = 0; m < materials.size() && !is_internal; ++m) {
is_internal |= materials[m]->isInternal<Real>(field_name, element_kind);
}
return is_internal;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap<UInt>
SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
const ElementKind & element_kind) {
if (!(this->isInternal(field_name, element_kind)))
AKANTU_EXCEPTION("unknown internal " << field_name);
for (UInt m = 0; m < materials.size(); ++m) {
if (materials[m]->isInternal<Real>(field_name, element_kind))
return materials[m]->getInternalDataPerElem<Real>(field_name,
element_kind);
}
return ElementTypeMap<UInt>();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray<Real> &
SolidMechanicsModel::flattenInternal(const std::string & field_name,
const ElementKind & kind,
const GhostType ghost_type) {
std::pair<std::string, ElementKind> key(field_name, kind);
if (this->registered_internals.count(key) == 0) {
this->registered_internals[key] =
new ElementTypeMapArray<Real>(field_name, this->id, this->memory_id);
}
ElementTypeMapArray<Real> * internal_flat = this->registered_internals[key];
typedef Mesh::type_iterator iterator;
iterator tit = this->mesh.firstType(spatial_dimension, ghost_type, kind);
iterator end = this->mesh.lastType(spatial_dimension, ghost_type, kind);
for (; tit != end; ++tit) {
ElementType type = *tit;
if (internal_flat->exists(type, ghost_type)) {
(*internal_flat)(type, ghost_type).clear();
(*internal_flat)(type, ghost_type).resize(0);
}
}
for (UInt m = 0; m < materials.size(); ++m) {
if (materials[m]->isInternal<Real>(field_name, kind))
materials[m]->flattenInternal(field_name, *internal_flat, ghost_type,
kind);
}
return *internal_flat;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::flattenAllRegisteredInternals(
const ElementKind & kind) {
typedef std::map<std::pair<std::string, ElementKind>,
ElementTypeMapArray<Real> *>::iterator iterator;
iterator it = this->registered_internals.begin();
iterator end = this->registered_internals.end();
while (it != end) {
if (kind == it->first.second)
this->flattenInternal(it->first.first, kind);
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onDump() {
this->flattenAllRegisteredInternals(_ek_regular);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper::Field * SolidMechanicsModel::createElementalField(
const std::string & field_name, const std::string & group_name,
bool padding_flag, const UInt & spatial_dimension,
const ElementKind & kind) {
dumper::Field * field = NULL;
if (field_name == "partitions")
field = mesh.createElementalField<UInt, dumper::ElementPartitionField>(
mesh.getConnectivities(), group_name, spatial_dimension, kind);
else if (field_name == "material_index")
field = mesh.createElementalField<UInt, Vector, dumper::ElementalField>(
material_index, group_name, spatial_dimension, kind);
else {
// this copy of field_name is used to compute derivated data such as
// strain and von mises stress that are based on grad_u and stress
std::string field_name_copy(field_name);
if (field_name == "strain" || field_name == "Green strain" ||
field_name == "principal strain" ||
field_name == "principal Green strain")
field_name_copy = "grad_u";
else if (field_name == "Von Mises stress")
field_name_copy = "stress";
bool is_internal = this->isInternal(field_name_copy, kind);
if (is_internal) {
ElementTypeMap<UInt> nb_data_per_elem =
this->getInternalDataPerElem(field_name_copy, kind);
ElementTypeMapArray<Real> & internal_flat =
this->flattenInternal(field_name_copy, kind);
field = mesh.createElementalField<Real, dumper::InternalMaterialField>(
internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem);
if (field_name == "strain") {
dumper::ComputeStrain<false> * foo =
new dumper::ComputeStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "Von Mises stress") {
dumper::ComputeVonMisesStress * foo =
new dumper::ComputeVonMisesStress(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "Green strain") {
dumper::ComputeStrain<true> * foo =
new dumper::ComputeStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "principal strain") {
dumper::ComputePrincipalStrain<false> * foo =
new dumper::ComputePrincipalStrain<false>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
} else if (field_name == "principal Green strain") {
dumper::ComputePrincipalStrain<true> * foo =
new dumper::ComputePrincipalStrain<true>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
// treat the paddings
if (padding_flag) {
if (field_name == "stress") {
if (spatial_dimension == 2) {
dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
} else if (field_name == "strain" || field_name == "Green strain") {
if (spatial_dimension == 2) {
dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
}
}
// homogenize the field
dumper::ComputeFunctorInterface * foo =
dumper::HomogenizerProxy::createHomogenizer(*field);
field = dumper::FieldComputeProxy::createFieldCompute(field, *foo);
}
}
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string, Array<Real> *> real_nodal_fields;
real_nodal_fields["displacement"] = displacement;
real_nodal_fields["mass"] = mass;
real_nodal_fields["velocity"] = velocity;
real_nodal_fields["acceleration"] = acceleration;
real_nodal_fields["force"] = force;
real_nodal_fields["residual"] = residual;
real_nodal_fields["increment"] = increment;
dumper::Field * field = NULL;
if (padding_flag)
field = mesh.createNodalField(real_nodal_fields[field_name], group_name, 3);
else
field = mesh.createNodalField(real_nodal_fields[field_name], group_name);
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(
const std::string & field_name, const std::string & group_name,
__attribute__((unused)) bool padding_flag) {
std::map<std::string, Array<bool> *> uint_nodal_fields;
uint_nodal_fields["blocked_dofs"] = blocked_dofs;
dumper::Field * field = NULL;
field = mesh.createNodalField(uint_nodal_fields[field_name], group_name);
return field;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createElementalField(
__attribute__((unused)) const std::string & field_name,
__attribute__((unused)) const std::string & group_name,
__attribute__((unused)) bool padding_flag,
__attribute__((unused)) const UInt & spatial_dimension,
__attribute__((unused)) const ElementKind & kind) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field *
SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
#endif
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, Real time,
UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump() {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeCauchyStresses() {
AKANTU_DEBUG_IN();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
if (mat.isFiniteDeformation())
mat.computeAllCauchyStresses(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::saveStressAndStrainBeforeDamage() {
EventManager::sendEvent(
SolidMechanicsModelEvent::BeginningOfDamageIterationEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateEnergiesAfterDamage() {
EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass->printself(stream, indent + 2);
velocity->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force->printself(stream, indent + 2);
residual->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + materials [" << std::endl;
std::vector<Material *>::const_iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
const Material & mat = *(*mat_it);
mat.printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
c4science · Help