diff --git a/src/fem/shape_linked.cc b/src/fem/shape_linked.cc index 68f221492..96697e610 100644 --- a/src/fem/shape_linked.cc +++ b/src/fem/shape_linked.cc @@ -1,47 +1,78 @@ /** * @file shape_linked.cc * * @author Nicolas Richart * @author Fabian Barras * * @date Fri Jul 15 19:41:58 2011 * * @brief ShapeLinked implementation * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "mesh.hh" #include "shape_linked.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template <> ShapeLinked<_ek_structural>::ShapeLinked(Mesh & mesh, const ID & id, const MemoryID & memory_id) : ShapeFunctions(mesh, id, memory_id) { } + +/* -------------------------------------------------------------------------- */ +template <> +ShapeLinked<_ek_structural>::~ShapeLinked() { + for (ghost_type_t::iterator gt = ghost_type_t::begin(); + gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + + // delete all the shapes id + ByElementTypeMultiReal::type_iterator s_type_it = + shapes.firstType(_all_dimensions, ghost_type, _ek_structural); + ByElementTypeMultiReal::type_iterator s_type_end = + shapes.lastType (_all_dimensions, ghost_type, _ek_structural); + + for(; s_type_it != s_type_end; ++s_type_it) { + delete [] shapes(*s_type_it, ghost_type); + } + + + // delete all the shapes derivatives id + ByElementTypeMultiReal::type_iterator sd_type_it = + shapes_derivatives.firstType(_all_dimensions, ghost_type, _ek_structural); + ByElementTypeMultiReal::type_iterator sd_type_end = + shapes_derivatives.lastType (_all_dimensions, ghost_type, _ek_structural); + + for(; sd_type_it != sd_type_end; ++sd_type_it) { + delete [] shapes_derivatives(*sd_type_it, ghost_type); + } + } +} + __END_AKANTU__ diff --git a/src/fem/shape_linked.hh b/src/fem/shape_linked.hh index a6aa5a5c4..6bb7a3e14 100644 --- a/src/fem/shape_linked.hh +++ b/src/fem/shape_linked.hh @@ -1,157 +1,157 @@ /** * @file shape_linked.hh * * @author Nicolas Richart * @author Fabian Barras * * @date Fri Jul 15 19:41:58 2011 * * @brief shape class for element with different set of shapes functions * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #include "shape_functions.hh" #ifndef __AKANTU_SHAPE_LINKED_HH__ #define __AKANTU_SHAPE_LINKED_HH__ __BEGIN_AKANTU__ template class ShapeLinked : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef ByElementType **> ByElementTypeMultiReal; ShapeLinked(Mesh & mesh, const ID & id = "shape_linked", const MemoryID & memory_id = 0); - virtual ~ShapeLinked(){}; + virtual ~ShapeLinked(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array & nodes, const Matrix & control_points, const ElementType & type, const GhostType & ghost_type); /// pre compute all shapes on the element control points from natural coordinates template void precomputeShapesOnControlPoints(const Array & nodes, const GhostType & ghost_type); /// pre compute all shapes on the element control points from natural coordinates template void precomputeShapeDerivativesOnControlPoints(const Array & nodes, const GhostType & ghost_type); /// interpolate nodal values on the control points template void interpolateOnControlPoints(const Array &u, Array &uq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter, bool accumulate = false, UInt id_shape = 0, UInt num_degre_of_freedom_to_interpolate = 0, UInt num_degre_of_freedom_interpolated = 0) const; /// compute the gradient of u on the control points template void gradientOnControlPoints(const Array &u, Array &nablauq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter, bool accumulate = false, UInt id_shape = 0, UInt num_degre_of_freedom_to_interpolate = 0, UInt num_degre_of_freedom_interpolated = 0) const; /// multiply a field by shape functions template void fieldTimesShapes(__attribute__((unused)) const Array & field, __attribute__((unused)) Array & fiedl_times_shapes, __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } private: template void extractNodalToElementField(const Array & nodal_f, Array & elemental_f, UInt num_degre_of_freedom_to_extract, const GhostType & ghost_type, const Array & filter_elements) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get a the shapes vector inline const Array & getShapes(const ElementType & type, const GhostType & ghost_type, UInt id = 0) const; /// get a the shapes derivatives vector inline const Array & getShapesDerivatives(const ElementType & type, const GhostType & ghost_type, UInt id = 0) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// shape functions for all elements ByElementTypeMultiReal shapes; /// shape derivatives for all elements ByElementTypeMultiReal shapes_derivatives; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "shape_linked_inline_impl.cc" #endif /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const ShapeLinked & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_SHAPE_LINKED_HH__ */ diff --git a/src/model/model.cc b/src/model/model.cc index 52c880c09..2292da359 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,214 +1,216 @@ /** * @file model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date Mon Oct 03 12:24:07 2011 * * @brief implementation of model common parts * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #include "model.hh" #include "element_group.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Model::Model(Mesh& m, UInt dim, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(m), spatial_dimension(dim == _all_dimensions ? m.getSpatialDimension() : dim), synch_registry(NULL),is_pbc_slave_node(0,1,"is_pbc_slave_node") { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +Model::~Model() { + AKANTU_DEBUG_IN(); + + FEMMap::iterator it; + for (it = fems.begin(); it != fems.end(); ++it) { + if(it->second) delete it->second; + } + + for (it = fems_boundary.begin(); it != fems_boundary.end(); ++it) { + if(it->second) delete it->second; + } + + delete synch_registry; + + delete dof_synchronizer; + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ void Model::initFull(std::string input_file, const ModelOptions & options) { AKANTU_DEBUG_IN(); if(input_file != "") { parser.parse(input_file); } initModel(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::createSynchronizerRegistry(DataAccessor * data_accessor){ synch_registry = new SynchronizerRegistry(*data_accessor); } /* -------------------------------------------------------------------------- */ void Model::setPBC(UInt x, UInt y, UInt z){ Mesh & mesh = getFEM().getMesh(); mesh.computeBoundingBox(); if (x) MeshUtils::computePBCMap(mesh, 0, pbc_pair); if (y) MeshUtils::computePBCMap(mesh, 1, pbc_pair); if (z) MeshUtils::computePBCMap(mesh, 2, pbc_pair); } /* -------------------------------------------------------------------------- */ void Model::setPBC(SurfacePairList & surface_pairs, ElementType surface_e_type){ Mesh & mesh = getFEM().getMesh(); SurfacePairList::iterator s_it; for(s_it = surface_pairs.begin(); s_it != surface_pairs.end(); ++s_it) { MeshUtils::computePBCMap(mesh, *s_it, surface_e_type, pbc_pair); } } /* -------------------------------------------------------------------------- */ void Model::initPBC() { Mesh & mesh = getFEM().getMesh(); std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); is_pbc_slave_node.resize(mesh.getNbNodes()); #ifndef AKANTU_NDEBUG Real * coords = mesh.getNodes().values; UInt dim = mesh.getSpatialDimension(); #endif while(it != end){ UInt i1 = (*it).first; is_pbc_slave_node(i1) = true; #ifndef AKANTU_NDEBUG UInt i2 = (*it).second; AKANTU_DEBUG_INFO("pairing " << i1 << " (" << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] << ") with " << i2 << " (" << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] << ")"); #endif ++it; } } /* -------------------------------------------------------------------------- */ DistributedSynchronizer & Model::createParallelSynch(MeshPartition * partition, __attribute__((unused)) DataAccessor * data_accessor){ AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); DistributedSynchronizer * synch = NULL; if(prank == 0) synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), partition); else synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), NULL); AKANTU_DEBUG_OUT(); return *synch; } -/* -------------------------------------------------------------------------- */ -Model::~Model() { - AKANTU_DEBUG_IN(); - - FEMMap::iterator it; - for (it = fems.begin(); it != fems.end(); ++it) { - if(it->second) delete it->second; - } - - for (it = fems_boundary.begin(); it != fems_boundary.end(); ++it) { - if(it->second) delete it->second; - } - - delete synch_registry; - - AKANTU_DEBUG_OUT(); -} - /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name, const std::string & dumper_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup() { GroupManager::element_group_iterator bit = mesh.element_group_begin(); GroupManager::element_group_iterator bend = mesh.element_group_end(); for(; bit != bend; ++bit) { bit->second->dump(); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory) { GroupManager::element_group_iterator bit = mesh.element_group_begin(); GroupManager::element_group_iterator bend = mesh.element_group_end(); for (; bit != bend; ++bit) { bit->second->setDirectory(directory); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setGroupBaseName(const std::string & basename, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setBaseName(basename); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Model::getGroupDumper(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); return group.getDumper(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 8b6736853..42559a94b 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1935 +1,1933 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" # include "dumper_iohelper_tmpl.hh" # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.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), Dumpable(), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_index_by_material("element index by material", id), material_selector(new DefaultMaterialSelector(element_index_by_material)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = 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->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; - if(solver) delete solver; - if(mass_matrix) delete mass_matrix; - if(velocity_damping_matrix) delete velocity_damping_matrix; + delete solver; + delete mass_matrix; + delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; - if(jacobian_matrix) delete jacobian_matrix; - - if(dof_synchronizer) delete dof_synchronizer; + delete jacobian_matrix; - if(synch_parallel) delete synch_parallel; + delete synch_parallel; - if(is_default_material_selector) delete material_selector; + delete material_selector; AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) 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(std::string input_file, const ModelOptions & options) { Model::initFull(input_file, options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the materials if(input_file != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary() { FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_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(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, 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 << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(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); element_index_by_material.alloc(nb_element, 2, *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 getFEM().initShapeFunctions(_not_ghost); getFEM().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::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); 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->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// 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->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* 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::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #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::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::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 * Ma = new Array(*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->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_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 * Cv = new Array(*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, *boundary, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_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, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; 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, *boundary, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); 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) // 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"; jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); 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 } /* -------------------------------------------------------------------------- */ /** * 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(*boundary); 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::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); // updateResidualInternal(); solve(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); increment->clear(); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->storage(); Real * displacement_val = displacement->storage(); bool * boundary_val = boundary->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } else { *increment_val = 0.0; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic(Array & boundary_normal, Array & EulerAngles) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); Array * residual_rotated = new Array (nb_nodes, nb_degree_of_freedom, "residual_rotated"); //stiffness_matrix->saveMatrix("stiffness_original.out"); jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated); //jacobian_matrix->saveMatrix("stiffness_rotated_dir.out"); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual_rotated); delete residual_rotated; if (!increment) setIncrementFlagOn(); solver->solve(*increment); Matrix T(nb_degree_of_freedom, nb_degree_of_freedom); Matrix small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Matrix T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool constrain_ij = false; for (UInt j = 0; j < nb_degree_of_freedom; j++) { if (boundary_normal(n, j)) { constrain_ij = true; break; } } if (constrain_ij) { if (nb_degree_of_freedom == 2) { Real Theta = EulerAngles(n, 0); T(0, 0) = cos(Theta); T(0, 1) = -sin(Theta); T(1, 1) = cos(Theta); T(1, 0) = sin(Theta); } else if (nb_degree_of_freedom == 3) { Real Theta_x = EulerAngles(n, 0); Real Theta_y = EulerAngles(n, 1); Real Theta_z = EulerAngles(n, 2); T(0, 0) = cos(Theta_y) * cos(Theta_z); T(0, 1) = -cos(Theta_y) * sin(Theta_z); T(0, 2) = sin(Theta_y); T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y); T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z); T(1, 2) = -cos(Theta_y) * sin(Theta_x); T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y); T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z); T(2, 2) = cos(Theta_x) * cos(Theta_y); } small_rhs.clear(); T_small_rhs.clear(); for (UInt j = 0; j < nb_degree_of_freedom; j++) if(!(boundary_normal(n,j)) ) small_rhs(j,j)=increment_val[j]; T_small_rhs.mul(T,small_rhs); for (UInt j = 0; j < nb_degree_of_freedom; j++){ if(!(boundary_normal(n,j))){ for (UInt k = 0; k < nb_degree_of_freedom; k++) displacement_val[k]+=T_small_rhs(k,j); } } displacement_val += nb_degree_of_freedom; boundary_val += nb_degree_of_freedom; increment_val += nb_degree_of_freedom; } else { for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } } } } 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){ 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 * boundary_val = boundary->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(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().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"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; 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(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().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(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->values; Real * disp_val = displacement->values; bool * boun_val = boundary->values; 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->values; Real * disp_val = displacement->values; Real * prev_disp_val = previous_displacement->values; 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(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::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::iterator< Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Array X(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::iterator< Matrix > X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = el; Real el_size = getFEM().getElementInradius(*X_el, *it); Real el_dt = mat_val[(*eibm)(0)]->getStableTimeStep(el_size, elem); min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ 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->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(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 = getFEM().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::iterator< Vector > vit = vel_on_quad.begin(spatial_dimension); Array::iterator< Vector > vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[element_index_by_material(type)(index, 0)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEM().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = boundary->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); 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::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::iterator mat_it; Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(0)]->getEnergy(energy_id, type, mat(1)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & 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(boundary ) boundary ->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); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_lumped_mass) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(mat_id_vect) mat_id = mat_id_vect[el]; else mat_id = (*material_selector)(elem); Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); Vector id(2); id[0] = mat_id; id[1] = mat_index; if(!element_index_by_material.exists(elem.type, elem.ghost_type)) element_index_by_material.alloc(0, 2, elem.type, elem.ghost_type); element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::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) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & 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(boundary ) mesh.removeNodesFromArray(*boundary , 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(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::reassignMaterial() { AKANTU_DEBUG_IN(); std::vector< Array > element_to_add (materials.size()); std::vector< Array > element_to_remove(materials.size()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type = *it; element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & el_index_by_mat = element_index_by_material(type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt old_material = el_index_by_mat(el, 0); UInt new_material = (*material_selector)(element); if(old_material != new_material) { element_to_add [new_material].push_back(element); element_to_remove[old_material].push_back(element); } } } } std::vector::iterator mat_it; UInt mat_index = 0; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it, ++mat_index) { (*mat_it)->removeElements(element_to_remove[mat_index]); (*mat_it)->addElements (element_to_add[mat_index]); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { 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(field_id); } if (is_internal) { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::HomogenizedField (*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } else { try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} } } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ group.addDumpFieldExternalToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField(mesh, dim, _not_ghost, _ek_regular, elemental_filter)); } else if(field_id == "element_index_by_material") { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, dim, _not_ghost, _ek_regular, elemental_filter)); } else { typedef DumperIOHelper::HomogenizedField Field; group.addDumpFieldExternalToDumper(dumper_name, field_id, new Field(*this, field_id, dim, _not_ghost, _ek_regular, elemental_filter)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter); \ f->setPadding(3); \ group.addDumpFieldExternalToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, dim, _not_ghost, _ek_regular, elemental_filter); field->setPadding(3); group.addDumpFieldExternalToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else if(field_id == "strain") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #endif } /* -------------------------------------------------------------------------- */ 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; getFEM().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); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index ac64555d9..ce09336f6 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,666 +1,668 @@ /** * @file structural_mechanics_model.cc * @author Fabian Barras * @date Thu May 5 15:52:38 2011 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" #include "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), stiffness_matrix(NULL), + jacobian_matrix(NULL), solver(NULL), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEMObject("StructuralMechanicsFEM", mesh, spatial_dimension); this->displacement_rotation = NULL; this->force_momentum = NULL; - this->residual = NULL; - this->boundary = NULL; - this->increment = NULL; + this->residual = NULL; + this->boundary = NULL; + this->increment = NULL; if(spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_DEBUG_TO_IMPLEMENT(); } this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() { AKANTU_DEBUG_IN(); - if(solver) delete solver; - if(stiffness_matrix) delete stiffness_matrix; + delete solver; + delete stiffness_matrix; + delete jacobian_matrix; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFull(std::string material) { initModel(); initArrays(); initSolver(); displacement_rotation->clear(); force_momentum ->clear(); residual ->clear(); boundary ->clear(); increment ->clear(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { computeRotationMatrix(*it); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; std::stringstream sstr_incr; sstr_incr << id << ":increment"; displacement_rotation = &(alloc(sstr_disp.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); force_momentum = &(alloc(sstr_forc.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, nb_degree_of_freedom, false)); increment = &(alloc(sstr_incr.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it); element_material.alloc(nb_element, 1, *it, _not_ghost); set_ID.alloc(nb_element, 1, *it, _not_ghost); UInt size = getTangentStiffnessVoigtSize(*it); stress.alloc(nb_element * nb_quadrature_points, size , *it, _not_ghost); } dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(), nb_degree_of_freedom); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { getFEM().initShapeFunctions(_not_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { AKANTU_DEBUG_IN(); const Mesh & mesh = getFEM().getMesh(); #if !defined(AKANTU_USE_MUMPS) // 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_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * nb_degree_of_freedom, _symmetric, nb_degree_of_freedom, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(options); #endif //AKANTU_HAS_SOLVER AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt StructuralMechanicsModel::getTangentStiffnessVoigtSize(const ElementType & type) { UInt size; #define GET_TANGENT_STIFFNESS_VOIGT_SIZE(type) \ size = getTangentStiffnessVoigtSize(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_TANGENT_STIFFNESS_VOIGT_SIZE); #undef GET_TANGENT_STIFFNESS_VOIGT_SIZE return size; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define ASSEMBLE_STIFFNESS_MATRIX(type) \ assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>(Array & rotations){ ElementType type = _bernoulli_beam_2; Mesh & mesh = getFEM().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::iterator< Vector > nodes_it = mesh.getNodes().begin(spatial_dimension); Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) { Matrix & R = *R_it; Vector & connec = *connec_it; Vector x2; x2 = nodes_it[connec(1)]; // X2 Vector x1; x1 = nodes_it[connec(0)]; // X1 Real le = x1.distance(x2); Real c = (x2(0) - x1(0)) / le; Real s = (x2(1) - x1(1)) / le; /// Definition of the rotation matrix R(0,0) = c; R(0,1) = s; R(0,2) = 0.; R(1,0) = -s; R(1,1) = c; R(1,2) = 0.; R(2,0) = 0.; R(2,1) = 0.; R(2,2) = 1.; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_3>(Array & rotations){ ElementType type = _bernoulli_beam_3; Mesh & mesh = getFEM().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::iterator< Vector > n_it = mesh.getNormals(type).begin(spatial_dimension); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::iterator< Vector > nodes_it = mesh.getNodes().begin(spatial_dimension); Matrix Pe (spatial_dimension, spatial_dimension); Matrix Pg (spatial_dimension, spatial_dimension); Matrix inv_Pg(spatial_dimension, spatial_dimension); Vector x_n(spatial_dimension); // x vect n Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e=0 ; e < nb_element; ++e, ++n_it, ++connec_it, ++R_it) { Vector & n = *n_it; Matrix & R = *R_it; Vector & connec = *connec_it; Vector x; x = nodes_it[connec(1)]; // X2 Vector y; y = nodes_it[connec(0)]; // X1 Real l = x.distance(y); x -= y; // X2 - X1 x_n.crossProduct(x, n); Pe.eye(); Pe(0, 0) *= l; Pe(1, 1) *= -l; Pg(0,0) = x(0); Pg(0,1) = x_n(0); Pg(0,2) = n(0); Pg(1,0) = x(1); Pg(1,1) = x_n(1); Pg(1,2) = n(1); Pg(2,0) = x(2); Pg(2,1) = x_n(2); Pg(2,2) = n(2); inv_Pg.inverse(Pg); Pe *= inv_Pg; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { R(i, j) = Pe(i, j); R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j); } } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEM().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if(!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom*nb_nodes_per_element * nb_degree_of_freedom*nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Arrayrotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) \ computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); Array::iterator< Matrix > T_it = rotation_matrix(type).begin(nb_degree_of_freedom*nb_nodes_per_element, nb_degree_of_freedom*nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { Matrix & T = *T_it; Matrix & R = *R_it; T.clear(); for (UInt k = 0; k < nb_nodes_per_element; ++k){ for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k*nb_degree_of_freedom + i, k*nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define COMPUTE_STRESS_ON_QUAD(type) \ computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); residual->copy(*force_momentum); Array ku(*displacement_rotation, true); ku *= *stiffness_matrix; *residual -= ku; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::solve() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving an implicit step."); UInt nb_nodes = displacement_rotation->getSize(); /// todo residual = force - Kxr * d_bloq jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement_rotation->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); Mesh & mesh= getFEM().getMesh(); UInt nb_nodes = displacement_rotation->getSize(); UInt nb_degree_of_freedom = displacement_rotation->getNbComponent(); Real norm = 0; Real * increment_val = increment->values; bool * boundary_val = boundary->values; 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(!(*boundary_val) && is_local_node) { norm += *increment_val * *increment_val; } boundary_val++; increment_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); error = sqrt(norm); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>(Array & tangent_moduli) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_2); UInt tangent_size = 2; Array::iterator< Matrix > D_it = tangent_moduli.begin(tangent_size, tangent_size); Array & el_mat = element_material(_bernoulli_beam_2, _not_ghost); for (UInt e = 0; e < nb_element; ++e) { UInt mat = el_mat(e); Real E = materials[mat].E; Real A = materials[mat].A; Real I = materials[mat].I; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * I; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_3>(Array & tangent_moduli) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_3); UInt tangent_size = 4; Array::iterator< Matrix > D_it = tangent_moduli.begin(tangent_size, tangent_size); for (UInt e = 0; e < nb_element; ++e) { UInt mat = element_material(_bernoulli_beam_3, _not_ghost)(e); Real E = materials[mat].E; Real A = materials[mat].A; Real Iz = materials[mat].Iz; Real Iy = materials[mat].Iy; Real GJ = materials[mat].GJ; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * Iz; D(2,2) = E * Iy; D(3,3) = GJ; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_2>(Array & b, bool local) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_2); MyFEMType & fem = getFEMClass(); Array::const_iterator< Vector > shape_Np = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_2>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::iterator< Matrix > B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,3) = Np(1); B(1,1) = Mpp(0); B(1,2) = Lpp(0); B(1,4) = Mpp(1); B(1,5) = Lpp(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } // ++R_it; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_3>(Array & b, bool local) { MyFEMType & fem = getFEMClass(); UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_3); Array::const_iterator< Vector > shape_Np = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_3>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::iterator< Matrix > B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,6) = Np(1); B(1,1) = Mpp(0); B(1,5) = Lpp(0); B(1,7) = Mpp(1); B(1,11) = Lpp(1); B(2,2) = Mpp(0); B(2,4) = -Lpp(0); B(2,8) = Mpp(1); B(2,10) = -Lpp(1); B(3,3) = Np(0); B(3,9) = Np(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ internalAddDumpFieldToDumper(dumper_name, id, \ new DumperIOHelper::NodalField(*field, n, stride)) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, n, 0); } else if(field_id == "rotation") { ADD_FIELD(dumper_name, "rotation", displacement_rotation, Real, nb_degree_of_freedom - n, n); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } else if(field_id == "momentum") { ADD_FIELD(dumper_name, "momentum", force_momentum, Real, nb_degree_of_freedom - n, n); } else if(field_id == "residual") { ADD_FIELD(dumper_name, "residual", residual, Real, nb_degree_of_freedom, 0); } else if(field_id == "boundary") { ADD_FIELD(dumper_name, "boundary", boundary, bool, nb_degree_of_freedom, 0); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_material, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, n, stride); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, id, f) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, n, 0); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldTensorToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { } __END_AKANTU__ diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.verified b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.verified new file mode 100644 index 000000000..a50bc1c64 --- /dev/null +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.verified @@ -0,0 +1,4 @@ + d1 = 0.00116667 + d2 = -0.000770833 + M1 = 6405.44 + M2 = 12590.9