diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index cee0bb30f..b134ffc57 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1510 +1,1513 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "aka_math.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "group_manager_inline_impl.cc" #include "static_communicator.hh" #include "synchronizer_registry.hh" #include "sparse_matrix.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_field.hh" #include "dumper_homogenizing_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper.hh" #include "dumper_material_padders.hh" #include "dumper_paraview.hh" #endif #ifdef AKANTU_DAMAGE_NON_LOCAL #include "non_local_manager.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition(), f_m2a(1.0), displacement(NULL), previous_displacement(NULL), displacement_increment(NULL), mass(NULL), velocity(NULL), acceleration(NULL), external_force(NULL), internal_force(NULL), blocked_dofs(NULL), current_position(NULL), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), increment_flag(false), are_materials_instantiated(false), non_local_manager(NULL) { //, pbc_synch(NULL) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); this->mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_material_id); this->registerSynchronizer(synchronizer, _gst_smm_mass); this->registerSynchronizer(synchronizer, _gst_smm_stress); this->registerSynchronizer(synchronizer, _gst_smm_boundary); this->registerSynchronizer(synchronizer, _gst_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); if (is_default_material_selector) { delete material_selector; material_selector = NULL; } for (auto & internal : this->registered_internals) { delete internal.second; } #ifdef AKANTU_DAMAGE_NON_LOCAL delete non_local_manager; non_local_manager = NULL; #endif // delete pbc_synch; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); this->method = smm_options.analysis_method; // initialize the vectors this->initArrays(); if (!this->hasDefaultSolver()) this->initNewSolver(this->method); // initialize pbc if (this->pbc_pair.size() != 0) this->initPBC(); #ifdef AKANTU_DAMAGE_NON_LOCAL /// create the non-local manager object for non-local damage computations std::stringstream nl_manager_name; nl_manager_name << "NLManager" << this->id; this->non_local_manager = new NonLocalManager(*this, nl_manager_name.str(), this->memory_id); #endif // initialize the materials if (this->parser->getLastParsedFile() != "") { this->instantiateMaterials(); } if (!smm_options.no_init_materials) { this->initMaterials(); } - if (increment_flag) - this->initBC(*this, *displacement, *displacement_increment, - *external_force); - else - this->initBC(*this, *displacement, *external_force); + // if (increment_flag) + this->initBC(*this, *displacement, *displacement_increment, *external_force); + // else + // this->initBC(*this, *displacement, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } } return options; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; this->method = method; switch (this->method) { case _explicit_lumped_mass: { solver_name = "explicit_lumped"; tss_type = _tsst_dynamic_lumped; break; } case _explicit_consistent_mass: { solver_name = "explicit"; tss_type = _tsst_dynamic; break; } case _static: { solver_name = "static"; tss_type = _tsst_static; break; } case _implicit_dynamic: { solver_name = "implicit"; tss_type = _tsst_dynamic; break; } } if (!this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); this->setIntegrationScheme(solver_name, "displacement", options.integration_scheme_type["displacement"], options.solution_type["displacement"]); this->setDefaultSolver(solver_name); } } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModel::allocNodalField(Array *& array, __attribute__((unused)) UInt nb_component, const ID & name) { if (array == NULL) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":" << name; array = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, T())); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, __attribute__((unused)) NonLinearSolverType non_linear_solver_type) { DOFManager & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); + this->allocNodalField(this->previous_displacement, spatial_dimension, + "previous_displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); + dof_manager.registerDOFsPrevious("displacement", + *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_static) { if (!dof_manager.hasMatrix("K")) { dof_manager.getNewMatrix("K", _symmetric); } if (!dof_manager.hasMatrix("J")) { dof_manager.getNewMatrix("J", "K"); } } } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initParallel(MeshPartition & partition, // DataAccessor * data_accessor) // { // AKANTU_DEBUG_IN(); // if (data_accessor == NULL) // data_accessor = this; // synch_parallel = &createParallelSynch(partition, *data_accessor); // synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); // synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::initArraysPreviousDisplacment() { - AKANTU_DEBUG_IN(); +// void SolidMechanicsModel::initArraysPreviousDisplacment() { +// AKANTU_DEBUG_IN(); - this->setIncrementFlagOn(); - if (not this->previous_displacement) { - this->allocNodalField(this->previous_displacement, spatial_dimension, - "previous_displacement"); +// this->setIncrementFlagOn(); +// if (not this->previous_displacement) { +// this->allocNodalField(this->previous_displacement, spatial_dimension, +// "previous_displacement"); - this->getDOFManager().registerDOFsPrevious("displacement", - *this->previous_displacement); - } +// this->getDOFManager().registerDOFsPrevious("displacement", +// *this->previous_displacement); +// } - AKANTU_DEBUG_OUT(); -} +// AKANTU_DEBUG_OUT(); +// } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type, ghost_type); this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initPBC() { // Model::initPBC(); // registerPBCSynchronizer(); // // as long as there are ones on the diagonal of the matrix, we can put // // boudandary true for slaves // std::map::iterator it = pbc_pair.begin(); // std::map::iterator end = pbc_pair.end(); // UInt dim = mesh.getSpatialDimension(); // while (it != end) { // for (UInt i = 0; i < dim; ++i) // (*blocked_dofs)((*it).first, i) = true; // ++it; // } // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::registerPBCSynchronizer() { // pbc_synch = new PBCSynchronizer(pbc_pair); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_uv); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res); // synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump); // // changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, -1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleJacobian() { this->assembleStiffnessMatrix(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() {} /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() {} /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as F_{int} = \int_{\Omega} N * \sigma d\Omega@f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->clear(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(_gst_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(_gst_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasStiffnessMatrixChanged(); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").clear(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(spatial_dimension); auto cpos_end = this->current_position->end(spatial_dimension); auto disp_it = this->displacement->begin(spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initializeUpdateResidualData() { // AKANTU_DEBUG_IN(); // UInt nb_nodes = mesh.getNbNodes(); // internal_force->resize(nb_nodes); // /// copy the forces in residual for boundary conditions // this->getDOFManager().assembleToResidual("displacement", // *this->external_force); // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->waitEndSynchronize(_gst_smm_uv); // this->updateCurrentPosition(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::computeStresses() { // if (isExplicit()) { // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->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.computeAllStresses(_not_ghost); // } // #ifdef AKANTU_DAMAGE_NON_LOCAL // /* Computation of the non local part */ // this->non_local_manager->computeAllNonLocalStresses(); // #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); // } // } // } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // /** // * Initialize the solver and create the sparse matrices needed. // * // */ // void SolidMechanicsModel::initSolver(__attribute__((unused)) // SolverOptions & options) { // UInt nb_global_nodes = mesh.getNbGlobalNodes(); // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // _symmetric)); // // jacobian_matrix->buildProfile(mesh, *dof_synchronizer, // spatial_dimension); // if (!isExplicit()) { // delete stiffness_matrix; // std::stringstream sstr_sti; // sstr_sti << id << ":stiffness_matrix"; // stiffness_matrix = &(this->getDOFManager().getNewMatrix("stiffness", // _symmetric)); // } // if (solver) solver->initialize(options); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::initJacobianMatrix() { // // @todo make it more flexible: this is an ugly patch to treat the case of // non // // fix profile of the K matrix // delete jacobian_matrix; // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // "stiffness")); // std::stringstream sstr_solv; sstr_solv << id << ":solver"; // delete solver; // solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); // if(solver) // solver->initialize(_solver_no_options); // } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ // void SolidMechanicsModel::initImplicit(bool dynamic) { // AKANTU_DEBUG_IN(); // method = dynamic ? _implicit_dynamic : _static; // if (!increment) // setIncrementFlagOn(); // initSolver(); // // if(method == _implicit_dynamic) { // // if(integrator) delete integrator; // // integrator = new TrapezoidalRule2(); // // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { // return this->getDOFManager().getNewMatrix("velocity_damping", "jacobian"); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitPred() { // AKANTU_DEBUG_IN(); // if(previous_displacement) previous_displacement->copy(*displacement); // if(method == _implicit_dynamic) // integrator->integrationSchemePred(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs); // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitCorr() { // AKANTU_DEBUG_IN(); // if(method == _implicit_dynamic) { // integrator->integrationSchemeCorrDispl(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs, // *increment); // } else { // UInt nb_nodes = displacement->getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // Real * incr_val = increment->storage(); // Real * disp_val = displacement->storage(); // bool * boun_val = blocked_dofs->storage(); // for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, // ++boun_val){ // *incr_val *= (1. - *boun_val); // *disp_val += *incr_val; // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::updateIncrement() { - AKANTU_DEBUG_IN(); +// void SolidMechanicsModel::updateIncrement() { +// AKANTU_DEBUG_IN(); - auto incr_it = this->displacement_increment->begin(spatial_dimension); - auto incr_end = this->displacement_increment->end(spatial_dimension); - auto disp_it = this->displacement->begin(spatial_dimension); - auto prev_disp_it = this->previous_displacement->begin(spatial_dimension); +// auto incr_it = this->displacement_increment->begin(spatial_dimension); +// auto incr_end = this->displacement_increment->end(spatial_dimension); +// auto disp_it = this->displacement->begin(spatial_dimension); +// auto prev_disp_it = this->previous_displacement->begin(spatial_dimension); - for (; incr_it != incr_end; ++incr_it) { - *incr_it = *disp_it; - *incr_it -= *prev_disp_it; - } +// for (; incr_it != incr_end; ++incr_it) { +// *incr_it = *disp_it; +// *incr_it -= *prev_disp_it; +// } - AKANTU_DEBUG_OUT(); -} +// AKANTU_DEBUG_OUT(); +// } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::updatePreviousDisplacement() { - AKANTU_DEBUG_IN(); +// /* -------------------------------------------------------------------------- +// */ 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?"); +// AKANTU_DEBUG_ASSERT( +// previous_displacement, +// "The previous displacement has to be initialized." +// << " Are you working with Finite or Ineslactic deformations?"); - previous_displacement->copy(*displacement); +// previous_displacement->copy(*displacement); - AKANTU_DEBUG_OUT(); -} +// AKANTU_DEBUG_OUT(); +// } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (material->isInternal(field_name, _ek_regular)) for (UInt g = _not_ghost; g <= _ghost; ++g) { GhostType ghost_type = (GhostType)g; material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); this->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); this->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::setIncrementFlagOn() { - AKANTU_DEBUG_IN(); +// void SolidMechanicsModel::setIncrementFlagOn() { +// AKANTU_DEBUG_IN(); - if (!displacement_increment) { - this->allocNodalField(displacement_increment, spatial_dimension, - "displacement_increment"); +// if (!displacement_increment) { +// this->allocNodalField(displacement_increment, spatial_dimension, +// "displacement_increment"); - this->getDOFManager().registerDOFsIncrement("displacement", - *this->displacement_increment); - } +// this->getDOFManager().registerDOFsIncrement("displacement", +// *this->displacement_increment); +// } - increment_flag = true; +// increment_flag = true; - AKANTU_DEBUG_OUT(); -} +// 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, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(spatial_dimension); auto m_end = this->mass->end(spatial_dimension); auto v_it = this->velocity->begin(spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const Vector & v = *v_it; const Vector & m = *m_it; Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) mv2 += v(i) * v(i) * m(i); } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array Mv(nb_nodes, spatial_dimension); this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv); Array::const_vector_iterator mv_it = Mv.begin(spatial_dimension); Array::const_vector_iterator mv_end = Mv.end(spatial_dimension); Array::const_vector_iterator v_it = this->velocity->begin(spatial_dimension); for (; mv_it != mv_end; ++mv_it, ++v_it) { ekin += v_it->dot(*mv_it); } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } StaticCommunicator::getStaticCommunicator().allReduce(ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto incr_or_velo_it = this->velocity->begin(spatial_dimension); if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(spatial_dimension); } auto ext_force_it = external_force->begin(spatial_dimension); auto int_force_it = internal_force->begin(spatial_dimension); auto boun_it = blocked_dofs->begin(spatial_dimension); Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < this->spatial_dimension; ++i) { if (boun(i)) work -= int_force(i) * incr_or_velo(i); else work += ext_force(i) * incr_or_velo(i); } } } StaticCommunicator::getStaticCommunicator().allReduce(work, _so_sum); if (this->method != _static) work *= this->getTimeStep(); 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.; for (auto & material : materials) energy += material->getEnergy(energy_id); /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(energy, _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); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = this->mesh.getNbElement(type, ghost_type); if (!material_index.exists(type, ghost_type)) { this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } else { this->material_index(type, ghost_type).resize(nb_element); this->material_local_numbering(type, ghost_type).resize(nb_element); } } } ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (auto & elem : element_list) { if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); for (auto & material : materials) material->onElementsAdded(element_list, event); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( __attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ 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, 0.); if (mass) mass->resize(nb_nodes, 0.); if (velocity) velocity->resize(nb_nodes, 0.); if (acceleration) acceleration->resize(nb_nodes, 0.); if (external_force) external_force->resize(nb_nodes, 0.); if (internal_force) internal_force->resize(nb_nodes, 0.); if (blocked_dofs) blocked_dofs->resize(nb_nodes, 0.); if (previous_displacement) previous_displacement->resize(nb_nodes, 0.); if (displacement_increment) displacement_increment->resize(nb_nodes, 0.); for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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 (internal_force) mesh.removeNodesFromArray(*internal_force, new_numbering); if (external_force) mesh.removeNodesFromArray(*external_force, new_numbering); if (blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) mesh.removeNodesFromArray(*displacement_increment, new_numbering); if (previous_displacement) mesh.removeNodesFromArray(*previous_displacement, new_numbering); // if (method != _explicit_lumped_mass) { // this->initSolver(); // } } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind) { /// check if at least one material contains field_id as an internal for (auto & material : materials) { bool is_internal = material->isInternal(field_name, element_kind); if (is_internal) return true; } return false; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind) { if (!(this->isInternal(field_name, element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (auto & material : materials) { if (material->isInternal(field_name, element_kind)) return material->getInternalDataPerElem(field_name, element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type) { std::pair key(field_name, kind); if (this->registered_internals.count(key) == 0) { this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id, this->memory_id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, kind)) { if (internal_flat->exists(type, ghost_type)) { auto & internal = (*internal_flat)(type, ghost_type); // internal.clear(); internal.resize(0); } } for (auto & material : materials) { if (material->isInternal(field_name, kind)) material->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals( const ElementKind & kind) { ElementKind _kind; ID _id; for (auto & internal : this->registered_internals) { std::tie(_id, _kind) = internal.first; if (kind == _kind) this->flattenInternal(_id, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump() { this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * SolidMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = NULL; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, spatial_dimension, kind); else if (field_name == "material_index") field = mesh.createElementalField( material_index, group_name, spatial_dimension, kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy, kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy, kind); field = mesh.createElementalField( internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem); if (field_name == "strain") { dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Von Mises stress") { dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Green strain") { dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal Green strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } // treat the paddings if (padding_flag) { if (field_name == "stress") { if (spatial_dimension == 2) { dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } else if (field_name == "strain" || field_name == "Green strain") { if (spatial_dimension == 2) { dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } } // homogenize the field dumper::ComputeFunctorInterface * foo = dumper::HomogenizerProxy::createHomogenizer(*field); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["displacement"] = this->displacement; real_nodal_fields["mass"] = this->mass; real_nodal_fields["velocity"] = this->velocity; real_nodal_fields["acceleration"] = this->acceleration; real_nodal_fields["external_force"] = this->external_force; real_nodal_fields["internal_force"] = this->internal_force; real_nodal_fields["increment"] = this->displacement_increment; if (field_name == "force") { AKANTU_EXCEPTION("The 'force' field has been renamed in 'external_force'"); } else if (field_name == "residual") { AKANTU_EXCEPTION( "The 'residual' field has been replaced by 'internal_force'"); } dumper::Field * field = NULL; if (padding_flag) field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createElementalField( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, __attribute__((unused)) const ElementKind & kind) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass->printself(stream, indent + 2); velocity->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; for (auto & material : materials) { material->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 1f54828f1..a6f5ffb9e 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,803 +1,803 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "boundary_condition.hh" #include "data_accessor.hh" #include "dumpable.hh" #include "integrator_gauss.hh" #include "material_selector.hh" #include "mesh.hh" #include "model.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class DumperIOHelper; class NonLocalManager; } /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ struct SolidMechanicsModelOptions : public ModelOptions { SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool no_init_materials = false) : analysis_method(analysis_method), no_init_materials(no_init_materials) { } AnalysisMethod analysis_method; bool no_init_materials; }; extern const SolidMechanicsModelOptions default_solid_mechanics_model_options; class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public MeshEventHandler, public BoundaryCondition, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; typedef FEEngineTemplate MyFEEngineType; protected: typedef EventHandlerManager EventManager; public: SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull( const ModelOptions & options = default_solid_mechanics_model_options); /// initialize the fem object needed for boundary conditions void initFEEngineBoundary(); /// register the tags associated with the parallel synchronizer // virtual void initParallel(MeshPartition * partition, // DataAccessor * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// allocate all vectors - void initArraysPreviousDisplacment(); + //void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer // void initPBC(); /// initialize a new solver and sets it as the default one to use void initNewSolver(const AnalysisMethod & method); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// allocate an array if needed template void allocNodalField(Array *& array, UInt nb_component, const ID & name); /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC // void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer // void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces void assembleInternalForces(); private: /// callback for the solver, this adds f_{ext} - f_{int} to the residual virtual void assembleResidual(); /// callback for the solver, this assembles the stiffness matrix virtual void assembleJacobian(); /// callback for the solver, this is called at beginning of solve virtual void predictor(); /// callback for the solver, this is called at end of solve virtual void corrector(); /// Callback for the model to instantiate the matricees when needed virtual void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type); protected: /* ------------------------------------------------------------------------ */ virtual TimeStepSolverType getDefaultSolverType() const; /* ------------------------------------------------------------------------ */ virtual ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ // public: // /// initialize the stuff for the explicit scheme // void initExplicit(AnalysisMethod analysis_method = // _explicit_lumped_mass); // bool isExplicit() { // return method == _explicit_lumped_mass || // method == _explicit_consistent_mass; // } // /// initialize the array needed by updateResidual (residual, // current_position) // void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); // /// assemble the residual for the explicit scheme // virtual void updateResidual(bool need_initialize = true); // /** // * \brief compute the acceleration from the residual // * this function is the explicit equivalent to solveDynamic in implicit // * In the case of lumped mass just divide the residual by the mass // * In the case of not lumped mass call // solveDynamic<_acceleration_corrector> // */ // void updateAcceleration(); /// Update the increment of displacement - void updateIncrement(); + // void updateIncrement(); - /// Copy the actuel displacement into previous displacement - void updatePreviousDisplacement(); + // /// Copy the actuel displacement into previous displacement + // void updatePreviousDisplacement(); // /// Save stress and strain through EventManager // void saveStressAndStrainBeforeDamage(); // /// Update energies through EventManager // void updateEnergiesAfterDamage(); // /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix // void solveLumped(Array & x, const Array & A, // const Array & b, const Array & blocked_dofs, // Real alpha); // /// explicit integration predictor // void explicitPred(); // /// explicit integration corrector // void explicitCorr(); // public: // void solveStep(); // /* // ------------------------------------------------------------------------ // */ // /* Implicit */ // /* // ------------------------------------------------------------------------ // */ // public: // /// initialize the solver and the jacobian_matrix (called by // initImplicit) // void initSolver(); // /// initialize the stuff for the implicit solver // void initImplicit(bool dynamic = false); // /// solve Ma = f to get the initial acceleration // void initialAcceleration(); // /// assemble the stiffness matrix // void assembleStiffnessMatrix(); // public: // /** // * solve a step (predictor + convergence loop + corrector) using the // * the given convergence method (see akantu::SolveConvergenceMethod) // * and the given convergence criteria (see // * akantu::SolveConvergenceCriteria) // **/ // template // bool solveStep(Real tolerance, UInt max_iteration = 100); // template // bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100, // bool do_not_factorize = false); // public: // /** // * solve Ku = f using the the given convergence method (see // * akantu::SolveConvergenceMethod) and the given convergence // * criteria (see akantu::SolveConvergenceCriteria) // **/ // template // bool solveStatic(Real tolerance, UInt max_iteration, // bool do_not_factorize = false); // /// create and return the velocity damping matrix // SparseMatrix & initVelocityDampingMatrix(); // /// implicit time integration predictor // void implicitPred(); // /// implicit time integration corrector // void implicitCorr(); // /// compute the Cauchy stress on user demand. // void computeCauchyStresses(); // // /// compute A and solve @f[ A\delta u = f_ext - f_int @f] // // template // // void solve(Array &increment, Real block_val = 1., // // bool need_factorize = true, bool has_profile_changed = // false); // protected: // /// finish the computation of residual to solve in increment // void updateResidualInternal(); // /// compute the support reaction and store it in force // void updateSupportReaction(); // private: // /// re-initialize the J matrix (to use if the profile of K changed) // void initJacobianMatrix(); /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion); public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation void computeStresses(); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// registers all the custom materials of a given type present in the input /// file template void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type template Material & registerNewEmptyMaterial(const std::string & name); // /// Use a UIntData in the mesh to specify the material to use per element // void setMaterialIDsFromIntData(const std::string & data_name); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database template Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials void assignMaterialToElements(const ElementTypeMapArray * filter = NULL); /// reinitialize dof_synchronizer and solver (either in implicit or /// explicit) when cohesive elements are inserted void reinitializeSolver(); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbData(const Array & elements, const SynchronizationTag & tag) const; inline virtual void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag); inline virtual UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const; inline virtual void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag); protected: inline void splitElementByMaterial(const Array & elements, Array * elements_per_mat) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); virtual void onElementsChanged( __attribute__((unused)) const Array & old_elements_list, __attribute__((unused)) const Array & new_elements_list, __attribute__((unused)) const ElementTypeMapArray & new_numbering, __attribute__((unused)) const ChangedElementsEvent & event){}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif virtual dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step // AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step virtual void setTimeStep(Real time_step, const ID & solver_id = ""); /// void setTimeStep(Real time_step); /// return the of iterations done in the last solveStep // AKANTU_GET_MACRO(NumberIter, n_iter, UInt); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (external forces) AKANTU_GET_MACRO(Force, *external_force, Array &); /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the SolidMechnicsModel::incrementAcceleration vector // AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, // Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ - void setIncrementFlagOn(); + //void setIncrementFlagOn(); // /// get the stiffness matrix // AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); // /// get the global jacobian matrix of the system // AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, const SparseMatrix // &); // /// get the mass matrix // AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); // /// get the velocity damping matrix // AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, // SparseMatrix &); /// get the FEEngine object to integrate or interpolate on the boundary inline FEEngine & getFEEngineBoundary(const ID & name = ""); // /// get integrator // AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder // &); // /// get synchronizer // AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const ElementSynchronizer &); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /// get the non-local manager AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// number of iterations // UInt n_iter; /// time step // Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement; /// increment of displacement Array * displacement_increment; /// lumped mass array Array * mass; /// velocities array Array * velocity; /// accelerations array Array * acceleration; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force; /// internal forces array Array * internal_force; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs; /// array of current position used during update residual Array * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; // /// integration scheme of second order used // IntegrationScheme2ndOrder *integrator; /// flag defining if the increment must be computed or not bool increment_flag; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// tells if the material are instantiated bool are_materials_instantiated; typedef std::map, ElementTypeMapArray *> flatten_internal_map; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// pointer to non-local manager: For non-local continuum damage computations NonLocalManager * non_local_manager; /// pointer to the pbc synchronizer // PBCSynchronizer * pbc_synch; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { typedef FromHigherDim FromStress; typedef FromSameDim FromTraction; } } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_material.cc b/src/model/solid_mechanics/solid_mechanics_model_material.cc index e24712635..62cd461f8 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,333 +1,333 @@ /** * @file solid_mechanics_model_material.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Nov 26 2010 * @date last modification: Mon Nov 16 2015 * * @brief instatiation of materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "material_list.hh" #include "solid_mechanics_model.hh" #ifdef AKANTU_DAMAGE_NON_LOCAL #include "non_local_manager.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL(dim, elem) \ registerNewMaterial> (section) #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL_EACH(r, data, i, elem) \ BOOST_PP_EXPR_IF(BOOST_PP_NOT_EQUAL(0, i), else) \ if (opt_param == BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem))) { \ registerNewMaterial> \ (section); \ } #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL(dim, elem) \ BOOST_PP_SEQ_FOR_EACH_I(AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL_EACH, \ (2, (dim, BOOST_PP_ARRAY_ELEM(1, elem))), \ BOOST_PP_ARRAY_ELEM(2, elem)) \ else { \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL(dim, elem); \ } #define AKANTU_INTANTIATE_MATERIAL_BY_DIM(dim, elem) \ BOOST_PP_IF(BOOST_PP_EQUAL(3, BOOST_PP_ARRAY_SIZE(elem)), \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL, \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL) \ (dim, elem) #define AKANTU_INTANTIATE_MATERIAL(elem) \ switch (spatial_dimension) { \ case 1: { \ AKANTU_INTANTIATE_MATERIAL_BY_DIM(1, elem); \ break; \ } \ case 2: { \ AKANTU_INTANTIATE_MATERIAL_BY_DIM(2, elem); \ break; \ } \ case 3: { \ AKANTU_INTANTIATE_MATERIAL_BY_DIM(3, elem); \ break; \ } \ } #define AKANTU_INTANTIATE_MATERIAL_IF(elem) \ if (mat_type == BOOST_PP_STRINGIZE(BOOST_PP_ARRAY_ELEM(0, elem))) { \ AKANTU_INTANTIATE_MATERIAL(elem); \ } #define AKANTU_INTANTIATE_OTHER_MATERIAL(r, data, elem) \ else AKANTU_INTANTIATE_MATERIAL_IF(elem) #define AKANTU_INSTANTIATE_MATERIALS() \ do { \ AKANTU_INTANTIATE_MATERIAL_IF(BOOST_PP_SEQ_HEAD(AKANTU_MATERIAL_LIST)) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_INTANTIATE_OTHER_MATERIAL, _, \ BOOST_PP_SEQ_TAIL(AKANTU_MATERIAL_LIST)) \ else { \ if (getStaticParser().isPermissive()) \ AKANTU_DEBUG_INFO("Malformed material file " \ << ": unknown material type '" << mat_type << "'"); \ else \ AKANTU_DEBUG_WARNING("Malformed material file " \ << ": unknown material type " << mat_type \ << ". This is perhaps a user" \ << " defined material ?"); \ } \ } while (0) /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::instantiateMaterials() { std::pair sub_sect = this->parser->getSubSections(_st_material); Parser::const_section_iterator it = sub_sect.first; for (; it != sub_sect.second; ++it) { const ParserSection & section = *it; std::string mat_type = section.getName(); std::string opt_param = section.getOption(); AKANTU_INSTANTIATE_MATERIALS(); } are_materials_instantiated = true; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assignMaterialToElements( const ElementTypeMapArray * filter) { Element element; element.ghost_type = _not_ghost; auto element_types = mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); if (filter != NULL) { element_types = filter->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); } // Fill the element material array from the material selector for (auto type : element_types) { UInt nb_element = mesh.getNbElement(type, _not_ghost); const Array * filter_array = NULL; if (filter != NULL) { filter_array = &((*filter)(type, _not_ghost)); nb_element = filter_array->getSize(); } element.type = type; element.kind = mesh.getKind(element.type); Array & mat_indexes = material_index(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { if (filter != NULL) element.element = (*filter_array)(el); else element.element = el; UInt mat_index = (*material_selector)(element); AKANTU_DEBUG_ASSERT( mat_index < materials.size(), "The material selector returned an index that does not exists"); mat_indexes(element.element) = mat_index; } } // synchronize the element material arrays this->synchronize(_gst_material_id); /// fill the element filters of the materials using the element_material /// arrays for (auto ghost_type : ghost_types) { element_types = mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined); if (filter != NULL) { element_types = filter->elementTypes(spatial_dimension, ghost_type, _ek_not_defined); } for (auto type : element_types) { UInt nb_element = mesh.getNbElement(type, ghost_type); const Array * filter_array = NULL; if (filter != NULL) { filter_array = &((*filter)(type, ghost_type)); nb_element = filter_array->getSize(); } Array & mat_indexes = material_index(type, ghost_type); Array & mat_local_num = material_local_numbering(type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { UInt element; if (filter != NULL) element = (*filter_array)(el); else element = el; UInt mat_index = mat_indexes(element); UInt index = materials[mat_index]->addElement(type, element, ghost_type); mat_local_num(element) = index; } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); if (!are_materials_instantiated) instantiateMaterials(); this->assignMaterialToElements(); for (auto & material : materials) { /// init internals properties material->initMaterial(); } this->synchronize(_gst_smm_init_mat); // initialize mass switch (method) { case _explicit_lumped_mass: assembleMassLumped(); break; case _explicit_consistent_mass: case _implicit_dynamic: assembleMass(); break; case _static: break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the previous displacement array if at least on material needs it - for (auto & material : materials) { - if (material->isFiniteDeformation() || material->isInelasticDeformation()) { - initArraysPreviousDisplacment(); - break; - } - } + // for (auto & material : materials) { + // if (material->isFiniteDeformation() || material->isInelasticDeformation()) { + // initArraysPreviousDisplacment(); + // break; + // } + // } #ifdef AKANTU_DAMAGE_NON_LOCAL /// initialize the non-local manager for non-local computations this->non_local_manager->init(); #endif } /* -------------------------------------------------------------------------- */ Int SolidMechanicsModel::getInternalIndexFromID(const ID & id) const { AKANTU_DEBUG_IN(); auto it = materials.begin(); auto end = materials.end(); for (; it != end; ++it) if ((*it)->getID() == id) { AKANTU_DEBUG_OUT(); return (it - materials.begin()); } AKANTU_DEBUG_OUT(); return -1; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::reassignMaterial() { AKANTU_DEBUG_IN(); std::vector> element_to_add(materials.size()); std::vector> element_to_remove(materials.size()); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { element.type = type; element.kind = Mesh::getKind(type); UInt nb_element = mesh.getNbElement(type, ghost_type); Array & mat_indexes = material_index(type, ghost_type); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt old_material = mat_indexes(el); 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); } } } } UInt mat_index = 0; for (auto 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::applyEigenGradU( const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type) { AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() == spatial_dimension * spatial_dimension, "The prescribed grad_u is not of the good size"); for (auto & material : materials) { if (material->getName() == material_name) material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type); } } /* -------------------------------------------------------------------------- */ __END_AKANTU__