diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 9a694f432..726e51965 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1494 +1,1477 @@ /** * @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 "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "element_group.hh" #include "static_communicator.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_paraview.hh" #include "dumper_homogenizing_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_elemental_field.hh" #include "dumper_material_padders.hh" #include "dumper_element_partition.hh" #include "dumper_iohelper.hh" #endif #ifdef AKANTU_DAMAGE_NON_LOCAL #include "non_local_manager.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ __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(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), - integrator(NULL), - increment_flag(false), solver(NULL), + increment_flag(false), synch_parallel(NULL), are_materials_instantiated(false), non_local_manager(NULL), pbc_synch(NULL) { AKANTU_DEBUG_IN(); this->createSynchronizerRegistry(this); 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 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 synch_parallel; if (is_default_material_selector) { delete material_selector; material_selector = NULL; } flatten_internal_map::iterator fl_it = this->registered_internals.begin(); flatten_internal_map::iterator fl_end = this->registered_internals.end(); for (; fl_it != fl_end; ++fl_it) { delete fl_it->second; } #ifdef AKANTU_DAMAGE_NON_LOCAL delete non_local_manager; non_local_manager = NULL; #endif delete pbc_synch; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* 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(); // set the initial condition to 0 external_force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pbc if (this->pbc_pair.size() != 0) this->initPBC(); // initialize the time integration schemes // switch (this->method) { // case _explicit_lumped_mass: // this->initExplicit(); // break; // case _explicit_consistent_mass: // this->initSolver(); // this->initExplicit(); // break; // case _implicit_dynamic: // this->initImplicit(true); // break; // case _static: // this->initImplicit(false); // this->initArraysPreviousDisplacment(); // break; // default: // AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); // break; // } #ifdef AKANTU_DAMAGE_NON_LOCAL /// create the non-local manager object for non-local damage computations std::stringstream nl_manager_name; nl_manager_name << "NLManager" << this->id; this->non_local_manager = new NonLocalManager(*this, nl_manager_name.str(), this->memory_id); #endif // initialize the materials if (this->parser->getLastParsedFile() != "") { this->instantiateMaterials(); } if (!smm_options.no_init_materials) { this->initMaterials(); } if (increment_flag) this->initBC(*this, *displacement, *increment, *external_force); else this->initBC(*this, *displacement, *external_force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition, data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { // AKANTU_DEBUG_IN(); // // in case of switch from implicit to explicit // if (!this->isExplicit()) // method = analysis_method; // // if (integrator) delete integrator; // // integrator = new CentralDifference(); // UInt nb_nodes = acceleration->getSize(); // UInt nb_degree_of_freedom = acceleration->getNbComponent(); // std::stringstream sstr; // sstr << id << ":increment_acceleration"; // increment_acceleration = // &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); this->setIncrementFlagOn(); if (not this->previous_displacement) { UInt nb_nodes = this->mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << this->id << ":previous_displacement"; this->previous_displacement = &(this->alloc( sstr_disp_t.str(), nb_nodes, this->spatial_dimension, 0.)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_ifor; sstr_forc << id << ":internal_force"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; /* ------------------------------------------------------------------------ */ // for static this->displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); this->internal_force = &(alloc(sstr_ifor.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); this->external_force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); this->blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); this->getDOFManager().registerDOFs("displacements", *this->displacement, _dst_nodal); this->getDOFManager().registerBlockedDOFs("displacements", *this->blocked_dofs); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; this->current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); /* ------------------------------------------------------------------------ */ // for dynamic std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; this->velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); this->acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); this->getDOFManager().registerDOFsDerivative("displacements", 1, *this->velocity); this->getDOFManager().registerDOFsDerivative("displacements", 2, *this->acceleration); 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); this->material_index.alloc(nb_element, 1, *it, gt); this->material_local_numbering.alloc(nb_element, 1, *it, gt); } } 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(); this->assembleInternalForces(); this->getDOFManager().assembleToResidual("displacements", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacements", *this->internal_force, -1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * 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(); 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.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleInternalForces(_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.assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleJacobian() { this->assembleStiffnessMatrix(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // 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::updateCurrentPosition() { AKANTU_DEBUG_IN(); this->current_position->copy(this->mesh.getNodes()); Array::vector_iterator cpos_it = this->current_position->begin(spatial_dimension); Array::vector_iterator cpos_end = this->current_position->end(spatial_dimension); Array::const_vector_iterator 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("displacements", // *this->external_force); // // start synchronization // synch_registry->asynchronousSynchronize(_gst_smm_uv); // synch_registry->waitEndSynchronize(_gst_smm_uv); // this->updateCurrentPosition(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // 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.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(); // AKANTU_DEBUG_ASSERT( // previous_displacement, // "The previous displacement has to be initialized." // << " Are you working with Finite or Ineslactic deformations?"); // UInt nb_nodes = displacement->getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // Real * incr_val = increment->storage(); // Real * disp_val = displacement->storage(); // Real * prev_disp_val = previous_displacement->storage(); // for (UInt j = 0; j < nb_degree_of_freedom; // ++j, ++disp_val, ++incr_val, ++prev_disp_val) // *incr_val = (*disp_val - *prev_disp_val); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::updatePreviousDisplacement() { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_ASSERT( // previous_displacement, // "The previous displacement has to be initialized." // << " Are you working with Finite or Ineslactic deformations?"); // previous_displacement->copy(*displacement); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry, "Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry, "Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if (!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(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, _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(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_regular); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type, _ek_regular); for (; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin(); Array::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin(); Array X(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::matrix_iterator X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, *it); Real el_c = mat_val[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(ekin, _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(); Real * incr_or_velo = NULL; if (this->method == _static) { incr_or_velo = this->increment->storage(); } else { incr_or_velo = this->velocity->storage(); } Real * forc = external_force->storage(); Real * resi = internal_force->storage(); bool * boun = blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = this->mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < this->spatial_dimension; ++i) { if (count_node) { if (*boun) work -= *resi * *incr_or_velo; else work += *forc * *incr_or_velo; } ++incr_or_velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(work, _so_sum); if (this->method != _static) work *= this->time_step; AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; std::vector::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, _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; 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 (UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType)g; Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); for (; it != end; ++it) { UInt nb_element = this->mesh.getNbElement(*it, gt); if (!material_index.exists(*it, gt)) { this->material_index.alloc(nb_element, 1, *it, gt); this->material_local_numbering.alloc(nb_element, 1, *it, gt); } else { this->material_index(*it, gt).resize(nb_element); this->material_local_numbering(*it, gt).resize(nb_element); } } } Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } // if (method == _explicit_lumped_mass) // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( __attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().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::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 (external_force) external_force->resize(nb_nodes); if (internal_force) internal_force->resize(nb_nodes); if (blocked_dofs) blocked_dofs->resize(nb_nodes); if (previous_displacement) previous_displacement->resize(nb_nodes); if (increment_acceleration) increment_acceleration->resize(nb_nodes); if (increment) increment->resize(nb_nodes); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->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 (increment) // mesh.removeNodesFromArray(*increment, new_numbering); // if (method != _explicit_lumped_mass) { // this->initSolver(); // } } -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::reinitializeSolver() { - - delete dof_synchronizer; - dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); - dof_synchronizer->initLocalDOFEquationNumbers(); - dof_synchronizer->initGlobalDOFEquationNumbers(); - - if (method != _explicit_lumped_mass) { - this->initSolver(); - } - - if (method == _explicit_lumped_mass) - this->assembleMassLumped(); -} - /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind) { bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_name, element_kind); } return is_internal; } /* -------------------------------------------------------------------------- */ 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 (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, element_kind)) return materials[m]->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]; typedef Mesh::type_iterator iterator; iterator tit = this->mesh.firstType(spatial_dimension, ghost_type, kind); iterator end = this->mesh.lastType(spatial_dimension, ghost_type, kind); for (; tit != end; ++tit) { ElementType type = *tit; if (internal_flat->exists(type, ghost_type)) { (*internal_flat)(type, ghost_type).clear(); (*internal_flat)(type, ghost_type).resize(0); } } for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, kind)) materials[m]->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals( const ElementKind & kind) { typedef std::map, ElementTypeMapArray *>::iterator iterator; iterator it = this->registered_internals.begin(); iterator end = this->registered_internals.end(); while (it != end) { if (kind == it->first.second) this->flattenInternal(it->first.first, kind); ++it; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump() { this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * SolidMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = NULL; if (field_name == "partitions") field = mesh.createElementalField( 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["force"] = this->external_force; real_nodal_fields["residual"] = this->internal_force; real_nodal_fields["increment"] = this->increment; 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, 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::computeCauchyStresses() { // AKANTU_DEBUG_IN(); // // call compute stiffness matrix on each local elements // std::vector::iterator mat_it; // for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { // Material & mat = **mat_it; // if (mat.isFiniteDeformation()) // mat.computeAllCauchyStresses(_not_ghost); // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { // EventManager::sendEvent( // SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::updateEnergiesAfterDamage() { // EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass->printself(stream, indent + 2); velocity->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); 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; std::vector::const_iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 1ffc6b9eb..4fc0f6a92 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,761 +1,761 @@ -n/** +/** * @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 "model.hh" #include "data_accessor.hh" #include "mesh.hh" #include "dumpable.hh" #include "boundary_condition.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "material_selector.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 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(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* 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(); /* ------------------------------------------------------------------------ */ /* 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(); // /// 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(); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ 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 getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline virtual void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline virtual void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, 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 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, *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(); // /// 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 DistributedSynchronizer &); 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 *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; /// internal synchronizer for parallel computations DistributedSynchronizer * synch_parallel; /// 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 "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #include "solid_mechanics_model_inline_impl.cc" /// 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_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index 63f668b9e..06943f214 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,437 +1,437 @@ -n/** +/** * @file solid_mechanics_model_inline_impl.cc * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Wed Aug 04 2010 * @date last modification: Wed Nov 18 2015 * * @brief Implementation of the inline functions 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 . * */ /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(const std::string & name) { AKANTU_DEBUG_IN(); std::map::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named "<< name); AKANTU_DEBUG_OUT(); return *materials[it->second]; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getMaterialIndex(const std::string & name) const { AKANTU_DEBUG_IN(); std::map::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named "<< name); AKANTU_DEBUG_OUT(); return it->second; } /* -------------------------------------------------------------------------- */ inline const Material & SolidMechanicsModel::getMaterial(const std::string & name) const { AKANTU_DEBUG_IN(); std::map::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named "<< name); AKANTU_DEBUG_OUT(); return *materials[it->second]; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::setMaterialSelector(MaterialSelector & selector) { if(is_default_material_selector) delete material_selector; material_selector = &selector; is_default_material_selector = false; } /* -------------------------------------------------------------------------- */ inline FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::splitElementByMaterial(const Array & elements, Array * elements_per_mat) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; const Array * mat_indexes = NULL; const Array * mat_loc_num = NULL; Array::const_iterator it = elements.begin(); Array::const_iterator end = elements.end(); for (; it != end; ++it) { Element el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; mat_indexes = &(this->material_index(el.type, el.ghost_type)); mat_loc_num = &(this->material_local_numbering(el.type, el.ghost_type)); } UInt old_id = el.element; el.element = (*mat_loc_num)(old_id); elements_per_mat[(*mat_indexes)(old_id)].push_back(el); } } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; Array::const_iterator it = elements.begin(); Array::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch(tag) { case _gst_material_id: { size += elements.getSize() * sizeof(UInt); break; } case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * spatial_dimension; // mass vector break; } case _gst_smm_for_gradu: { size += nb_nodes_per_element * spatial_dimension * sizeof(Real); // displacement break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case _gst_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * spatial_dimension * sizeof(Real) * 5; break; } default: { } } if(tag != _gst_material_id) { Array * elements_per_mat = new Array[materials.size()]; this->splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { size += materials[i]->getNbDataForElements(elements_per_mat[i], tag); } delete [] elements_per_mat; } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); switch(tag) { case _gst_material_id: { packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case _gst_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if(tag != _gst_material_id) { Array * elements_per_mat = new Array[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->packElementData(buffer, elements_per_mat[i], tag); } delete [] elements_per_mat; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { AKANTU_DEBUG_IN(); switch(tag) { case _gst_material_id: { unpackElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case _gst_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if(tag != _gst_material_id) { Array * elements_per_mat = new Array[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->unpackElementData(buffer, elements_per_mat[i], tag); } delete [] elements_per_mat; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * spatial_dimension; break; } case _gst_for_dump: { size += sizeof(Real) * spatial_dimension * 5; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * spatial_dimension; break; } case _gst_for_dump: { size += sizeof(Real) * spatial_dimension * 5; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { Array::const_vector_iterator it_disp = displacement->begin(spatial_dimension); Array::const_vector_iterator it_velo = velocity->begin(spatial_dimension); Vector disp(it_disp[index]); buffer << disp; Vector velo(it_velo[index]); buffer << velo; break; } case _gst_smm_res: { Array::const_vector_iterator it_res = internal_force->begin(spatial_dimension); Vector resi(it_res[index]); buffer << resi; break; } case _gst_smm_mass: { AKANTU_DEBUG_INFO("pack mass of node " << index << " which is " << (*mass)(index,0)); Array::const_vector_iterator it_mass = mass->begin(spatial_dimension); Vector mass(it_mass[index]); buffer << mass; break; } case _gst_for_dump: { Array::const_vector_iterator it_disp = displacement->begin(spatial_dimension); Array::const_vector_iterator it_velo = velocity->begin(spatial_dimension); Array::const_vector_iterator it_acce = acceleration->begin(spatial_dimension); Array::const_vector_iterator it_resi = internal_force->begin(spatial_dimension); Array::const_vector_iterator it_forc = external_force->begin(spatial_dimension); Vector disp(it_disp[index]); buffer << disp; Vector velo(it_velo[index]); buffer << velo; Vector acce(it_acce[index]); buffer << acce; Vector resi(it_resi[index]); buffer << resi; Vector forc(it_forc[index]); buffer << forc; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { Array::vector_iterator it_disp = displacement->begin(spatial_dimension); Array::vector_iterator it_velo = velocity->begin(spatial_dimension); Vector disp(it_disp[index]); buffer >> disp; Vector velo(it_velo[index]); buffer >> velo; break; } case _gst_smm_res: { Array::vector_iterator it_res = internal_force->begin(spatial_dimension); Vector res(it_res[index]); buffer >> res; break; } case _gst_smm_mass: { AKANTU_DEBUG_INFO("mass of node " << index << " was " << (*mass)(index,0)); Array::vector_iterator it_mass = mass->begin(spatial_dimension); Vector mass_v(it_mass[index]); buffer >> mass_v; AKANTU_DEBUG_INFO("mass of node " << index << " is now " << (*mass)(index,0)); break; } case _gst_for_dump: { Array::vector_iterator it_disp = displacement->begin(spatial_dimension); Array::vector_iterator it_velo = velocity->begin(spatial_dimension); Array::vector_iterator it_acce = acceleration->begin(spatial_dimension); Array::vector_iterator it_resi = internal_force->begin(spatial_dimension); Array::vector_iterator it_forc = external_force->begin(spatial_dimension); Vector disp(it_disp[index]); buffer >> disp; Vector velo(it_velo[index]); buffer >> velo; Vector acce(it_acce[index]); buffer >> acce; Vector resi(it_resi[index]); buffer >> resi; Vector forc(it_forc[index]); buffer >> forc; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ diff --git a/src/solver/solver_mumps.hh b/src/solver/solver_mumps.hh index ec3080f6e..fee890497 100644 --- a/src/solver/solver_mumps.hh +++ b/src/solver/solver_mumps.hh @@ -1,152 +1,152 @@ -n/** +/** * @file sparse_solver_mumps.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 19 2016 * * @brief Solver class implementation for the mumps solver * * @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 "sparse_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ namespace akantu { class DOFManagerDefault; class SparseMatrixAIJ; } __BEGIN_AKANTU__ class SparseSolverMumps : public SparseSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id = "sparse_solver_mumps", const MemoryID & memory_id = 0); virtual ~SparseSolverMumps(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile and do the analysis part virtual void initialize(); // void initializeSlave(SolverOptions & options = _solver_no_options); /// analysis (symbolic facto + permutations) virtual void analysis(); /// factorize the matrix virtual void factorize(); /// solve the system virtual void solve(); private: /// print the error if any happened in mumps void printError(); /// set internal values; void initMumpsData(); /// set the level of verbosity of mumps based on the debug level of akantu void setOutputLevel(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ private: /// access the control variable inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } /// access the results info inline Int & info(UInt i) { return mumps_data.info[i - 1]; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// DOFManager used by the Mumps implementation of the SparseSolver DOFManagerDefault & dof_manager; /// AIJ Matrix, usualy the jacobian matrix SparseMatrixAIJ & matrix; /// Right hand side per processors const Array & rhs; /// solution vector per processors Array & solution; /// Full right hand side on the master processors and solution after solve Array master_rhs_solution; /// mumps data DMUMPS_STRUC_C mumps_data; /// Rank of the current process UInt prank; /// matrix release at last solve UInt last_profile_release; /// matrix release at last solve UInt last_value_release; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverParallelMethod parallel_method; bool rhs_is_local; enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; __END_AKANTU__ #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/src/synchronizer/dof_synchronizer.hh b/src/synchronizer/dof_synchronizer.hh index a2ff6bcd8..90acb22ef 100644 --- a/src/synchronizer/dof_synchronizer.hh +++ b/src/synchronizer/dof_synchronizer.hh @@ -1,246 +1,245 @@ /** * @file dof_synchronizer.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Tue Dec 08 2015 * * @brief Synchronize Array of DOFs * * @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_common.hh" #include "aka_array.hh" #include "static_communicator.hh" #include "synchronizer.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_SYNCHRONIZER_HH__ #define __AKANTU_DOF_SYNCHRONIZER_HH__ __BEGIN_AKANTU__ class Mesh; template class AddOperation { public: inline T operator()(T & b, T & a) { return a + b; }; }; class DOFSynchronizer : public Synchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - DOFSynchronizer(const Mesh & mesh, UInt nb_degree_of_freedom, - StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator()); + DOFSynchronizer(const Mesh & mesh, UInt nb_degree_of_freedom); virtual ~DOFSynchronizer(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// asynchronous synchronization of ghosts virtual void asynchronousSynchronize(DataAccessor & data_accessor,SynchronizationTag tag); /// wait end of asynchronous synchronization of ghosts virtual void waitEndSynchronize(DataAccessor & data_accessor,SynchronizationTag tag); /// compute buffer size for a given tag and data accessor virtual void computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag); /// init the scheme for scatter and gather operation, need extra memory void initScatterGatherCommunicationScheme(); /// initialize the equation number with local ids void initLocalDOFEquationNumbers(); /// initialize the equation number with local ids void initGlobalDOFEquationNumbers(); /** * Gather the DOF value on the root proccessor * * @param to_gather data to gather * @param root processor on which data are gathered * @param gathered Array containing the gathered data, only valid on root processor */ template /// Gather the DOF value on the root proccessor void gather(const Array & to_gather, UInt root, Array * gathered = NULL) const; /** * Scatter a DOF Array form root to all processors * * @param scattered data to scatter, only valid on root processor * @param root processor scattering data * @param to_scatter result of scattered data */ template /// Scatter a DOF Array form root to all processors void scatter(Array & scattered, UInt root, const Array * to_scatter = NULL) const; template void synchronize(Array & vector) const ; template