diff --git a/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc b/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc index 254478c3f..3e868e31a 100644 --- a/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc +++ b/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc @@ -1,275 +1,276 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" #include "aka_random_generator.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "non_local_manager.hh" #include "parser.hh" #include "sparse_matrix.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), ASRTools(dynamic_cast(*this)), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), stiffness_changed(true) { AKANTU_DEBUG_IN(); RandomGenerator::seed(1); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); mesh.makePeriodic(_x, "left", "right"); mesh.makePeriodic(_y, "bottom", "top"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 50); solver.set("threshold", 1e-4); solver.set("convergence_type", SolveConvergenceCriteria::_solution); // this->initMaterials(); /// compute the volume of the RVE this->computeModelVolume(); std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id; // << this->memory_id - 1; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); this->addDumpField("delta_T"); this->addDumpField("extra_volume"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::assembleInternalForces() { /// displacement correction auto disp_begin = make_view(*this->displacement, spatial_dimension).begin(); auto u_top_corr = Vector(disp_begin[this->corner_nodes(2)]) - Vector(disp_begin[this->corner_nodes(1)]); auto u_right_corr = Vector(disp_begin[this->corner_nodes(2)]) - Vector(disp_begin[this->corner_nodes(3)]); auto correct_disp = [&](auto && node, auto && correction) { const auto & pair = mesh.getPeriodicMasterSlaves().equal_range(node); for (auto && data : range(pair.first, pair.second)) { auto slave = data.second; auto slave_disp = Vector(disp_begin[slave]); slave_disp = Vector(disp_begin[node]); slave_disp += correction; } }; for (auto node : bottom_nodes) { correct_disp(node, u_top_corr); } for (auto node : left_nodes) { correct_disp(node, u_right_corr); } SolidMechanicsModel::assembleInternalForces(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) { (*prestrain_it) = prestrain; } } /// advance the damage MaterialDamageIterativeInterface & mat_paste = dynamic_cast( this->getMaterial("paste")); MaterialDamageIterativeInterface & mat_aggregate = dynamic_cast( this->getMaterial("aggregate")); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; this->stiffness_changed = false; /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); /// save nodal fields (disp, boun & ext_force) this->storeNodalFields(); do { /// restore nodals and update grad_u accordingly this->restoreNodalFields(); /// restore historical internal fields (sigma_v for visc) this->restoreInternalFields(); this->solveStep(); /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) { nb_damaged_elements = mat_aggregate.updateDamage(); } else if (max_eq_stress_aggregate < max_eq_stress_paste) { nb_damaged_elements = mat_paste.updateDamage(); } else { nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); } - // mark the flag to update stiffness if elements were damaged - if (nb_damaged_elements) { + + /// mark the flag to update stiffness if elements were damaged + if (nb_damaged_elements != 0) { this->stiffness_changed = true; } std::cout << "Proc " << prank << " the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); // if (this->nb_dumps % 10 == 0) { // this->dump(); // } // this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) { instantiateMaterials(); } if (use_RVE_mat_selector) { auto tmp = std::make_shared( *this, "gel", this->nb_gel_pockets, "aggregate"); tmp->setFallback(material_selector); material_selector = tmp; } this->assignMaterialToElements(); // synchronize the element material arrays this->synchronize(SynchronizationTag::_material_id); for (auto & material : materials) { /// init internals properties const auto type = material->getID(); if (type.find("material_FE2") != std::string::npos) continue; material->initMaterial(); } this->synchronize(SynchronizationTag::_smm_init_mat); if (this->non_local_manager) { this->non_local_manager->initialize(); } // SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_orthotropic_inline_impl.hh b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_orthotropic_inline_impl.hh index 705414050..a2dcdc73e 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_orthotropic_inline_impl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_orthotropic_inline_impl.hh @@ -1,347 +1,346 @@ /** * @file material_iterative_stiffness_reduction.cc * @author Aurelia Isabel Cuba Ramos * @date Thu Feb 18 16:03:56 2016 * * @brief Implementation of material iterative stiffness reduction * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "material_damage_iterative_orthotropic.hh" #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialDamageIterativeOrthotropic:: MaterialDamageIterativeOrthotropic(SolidMechanicsModel & model, const ID & id) : parent(model, id), nb_state_changes("nb_state_changes", *this), damage_prev_iteration("damage_prev_iteration", *this), in_tension("in_tension", *this) { this->registerParam("max_state_changes_allowed", max_state_changes_allowed, Real(5), _pat_parsmod, "How many times an element can change between tension " "and compression stiffness"); this->registerParam("contact", contact, true, _pat_parsmod, "Account for contact by recovering initial stiffness in " "direction orthogonal to crack"); this->registerParam("iso_damage", iso_damage, false, _pat_parsmod, "Reduce stiffness in all directions"); this->registerParam("loading_test", loading_test, false, _pat_modifiable, "Indicate to material that it's in the loading test"); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeOrthotropic::initMaterial() { AKANTU_DEBUG_IN(); parent::initMaterial(); this->nb_state_changes.setDefaultValue(0.); this->nb_state_changes.initialize(1); this->damage_prev_iteration.initialize(1); this->in_tension.setDefaultValue(true); this->in_tension.initialize(1); this->E = this->E1; this->nu = this->nu12; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeOrthotropic::computeStress( const ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); PlaneStressToolbox>:: computeStress(el_type, ghost_type); /// if in loading test stop updating stiffness after 1st iteration Int nb_iter = this->model.getDOFManager().getNonLinearSolver("static").get( "nb_iterations"); if (not(nb_iter > 0 and this->loading_test)) { computeC(el_type, ghost_type); } auto C_it = this->C_field(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); auto sigma_th_it = make_view(this->sigma_th(el_type, ghost_type)).begin(); const auto & elem_filter = this->element_filter(el_type); auto & extra_vol = this->extra_volume(el_type); Array volumetric_strain( elem_filter.size() * this->fem.getNbIntegrationPoints(el_type), 1, 0.); auto vol_strain_it = volumetric_strain.begin(); auto dam_it = this->damage(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); /// compute stress according to anisotropic material law this->computeStressOnQuad(grad_u, sigma, *C_it, *sigma_th_it); /// compute volumetric strain in record it into extra_volume array if (*dam_it) { *vol_strain_it = grad_u.trace(); *vol_strain_it *= int(*vol_strain_it > 0); } ++C_it; ++sigma_th_it; ++vol_strain_it; ++dam_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; /// integrate volumetric strain across elements, subtract initial volume this->fem.integrate(volumetric_strain, extra_vol, 1, el_type, ghost_type, elem_filter); this->computeNormalizedEquivalentStress(el_type, ghost_type); this->norm_max_equivalent_stress = 0; this->findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialDamageIterativeOrthotropic::computeC( const ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * dam_prev_iter = this->damage_prev_iteration(el_type, ghost_type).storage(); auto E1_it = this->E1_field(el_type, ghost_type).begin(); auto E2_it = this->E2_field(el_type, ghost_type).begin(); auto E3_it = this->E3_field(el_type, ghost_type).begin(); auto nu12_it = this->nu12_field(el_type, ghost_type).begin(); auto nu13_it = this->nu13_field(el_type, ghost_type).begin(); auto nu23_it = this->nu23_field(el_type, ghost_type).begin(); auto G12_it = this->G12_field(el_type, ghost_type).begin(); auto G13_it = this->G13_field(el_type, ghost_type).begin(); auto G23_it = this->G23_field(el_type, ghost_type).begin(); auto Cprime_it = this->Cprime_field(el_type, ghost_type) .begin(spatial_dimension * spatial_dimension, spatial_dimension * spatial_dimension); auto C_it = this->C_field(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); auto eigC_it = this->eigC_field(el_type, ghost_type).begin(voigt_h::size); auto dir_vecs_it = this->dir_vecs_field(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto flick_it = this->nb_state_changes(el_type, ghost_type).begin(); auto tension_it = this->in_tension(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); /// parameters reduction & update of C and Cprime only if damage changed // if (*dam != *dam_prev_iter) { /// reduce or recover elastic moduli due to damage updateElasticModuli(sigma, grad_u, *dam, *E1_it, *E2_it, *E3_it, *nu12_it, *nu13_it, *nu23_it, *G12_it, *G13_it, *G23_it, *dir_vecs_it, *flick_it, *tension_it); /// construct the stiffness matrix with update parameters this->updateInternalParametersOnQuad( *E1_it, *E2_it, *E3_it, *nu12_it, *nu13_it, *nu23_it, *G12_it, *G13_it, *G23_it, *Cprime_it, *C_it, *eigC_it, *dir_vecs_it); // } /// update damage at previous iteration value *dam_prev_iter = *dam; ++dam; ++dam_prev_iter; ++E1_it; ++E2_it; ++E3_it; ++nu12_it; ++nu13_it; ++nu23_it; ++G12_it; ++G13_it; ++G23_it; ++Cprime_it; ++C_it; ++eigC_it; ++dir_vecs_it; ++flick_it; ++tension_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialDamageIterativeOrthotropic::updateElasticModuli( Matrix & /*sigma*/, Matrix & grad_u, Real & dam, Real & E1, Real & /*E2*/, Real & /*E3*/, Real & nu12, Real & nu13, Real & /*nu23*/, Real & G12, Real & G13, Real & /*G23*/, Matrix & dir_vecs, Real & nb_flicks, bool & in_tension) { if (dam == 0.) { E1 = this->E1; nu12 = this->nu12; nu13 = this->nu13; } else { Matrix strain(spatial_dimension, spatial_dimension); strain = 0.5 * (grad_u + grad_u.transpose()); Vector normal_to_crack(spatial_dimension); /// normals are stored row-wise for (auto i : arange(spatial_dimension)) { normal_to_crack(i) = dir_vecs(0, i); } auto def_on_crack = strain * normal_to_crack; auto def_normal_to_crack = normal_to_crack.dot(def_on_crack); // auto traction_on_crack = sigma * normal_to_crack; // auto stress_normal_to_crack = normal_to_crack.dot(traction_on_crack); /// coefficients of the cubic spline function // if (nb_flicks == this->max_state_changes_allowed) { // _E1 = std::sqrt(this->E1 * this->E1 * (1 - dam)); // _nu12 = std::sqrt(this->nu12 * this->nu12 * (1 - dam)); // _nu13 = std::sqrt(this->nu13 * this->nu13 * (1 - dam)); // // _G12 = std::sqrt(this->G12 * this->G12 * (1 - dam)); // // _G13 = std::sqrt(this->G13 * this->G13 * (1 - dam)); // std::cout << "Max number of state changes" << std::endl; // } else { if (def_normal_to_crack >= 0) { if (not in_tension) { ++nb_flicks; } in_tension = true; E1 = this->E1 * (1 - dam); nu12 = this->nu12 * (1 - dam); nu13 = this->nu13 * (1 - dam); G12 = this->G12 * (1 - dam); G13 = this->G13 * (1 - dam); } else { if (in_tension) { ++nb_flicks; } in_tension = false; // _E1 = std::min(-def_normal_to_crack * 1e13 + this->E1 * (1 - dam), // this->E1); // _nu12 = std::min(this->nu12 * _E1 / this->E1, this->nu12); // _nu13 = std::min(this->nu13 * _E1 / this->E1, this->nu13); E1 = this->E1; nu12 = this->nu12; nu13 = this->nu13; // _G12 = this->G12; // _G13 = this->G13; } } // } } -/* -------------------------------------------------------------------------- - */ +/* -------------------------------------------------------------------------- */ template inline void MaterialDamageIterativeOrthotropic::reduceInternalParameters( Matrix & sigma, Real & dam, Real & E1, Real & E2, Real & /*E3*/, Real & nu12, Real & nu13, Real & nu23, Real & /*G12*/, Real & /*G13*/, Real & /*G23*/, Matrix & dir_vecs, Real & nb_flicks) { /// detect compression in the normal to crack plane direction Vector normal_to_crack(spatial_dimension); /// normals are stored row-wise for (auto i : arange(spatial_dimension)) { normal_to_crack(i) = dir_vecs(0, i); } auto traction_on_crack = sigma * normal_to_crack; auto stress_normal_to_crack = normal_to_crack.dot(traction_on_crack); /// elements who exceed max allowed number of state changes get the average /// elastic properties if (nb_flicks == this->max_state_changes_allowed) { E1 = std::sqrt(this->E1 * this->E1 * (1 - dam)); nu12 = std::sqrt(this->nu12 * this->nu12 * (1 - dam)); nu13 = std::sqrt(this->nu13 * this->nu13 * (1 - dam)); // G12 = std::sqrt(this->G12 * this->G12 * (1 - dam)); // G13 = std::sqrt(this->G13 * this->G13 * (1 - dam)); } else { /// recover stiffness only when compressive stress is considerable if (this->contact && std::abs(stress_normal_to_crack) > this->E / 1e9 && stress_normal_to_crack < 0) { E1 = this->E1; nu12 = this->nu12; nu13 = this->nu13; // G12 = this->G12; // G13 = this->G13; ++nb_flicks; } else { E1 = this->E1 * (1 - dam); nu12 = this->nu12 * (1 - dam); nu13 = this->nu13 * (1 - dam); // _G12 = this->G12 * (1 - dam); // _G13 = this->G13 * (1 - dam); // _G12 = std::sqrt(this->G12 * this->G12 * (1 - dam)); // _G13 = std::sqrt(this->G13 * this->G13 * (1 - dam)); if (this->iso_damage) { E2 = this->E2 * (1 - dam); nu23 = this->nu23 * (1 - dam); } /// update shear moduli (Stable orthotropic materials (Li & Barbic)) // _G12 = std::sqrt(_E1 * _E2) / 2 / (1 + std::sqrt(_nu12 * // this->nu12)); _G13 = std::sqrt(_E1 * _E3) / 2 / (1 + std::sqrt(_nu13 // * this->nu13)); } } } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeOrthotropic< spatial_dimension>::computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto release = this->model.getDisplacementRelease(); - if (release != this->last_displacement_release) { + if (release != this->gradu_release(el_type, ghost_type)) { /// if in loading test stop updating stiffness after 1st iteration Int nb_iter = this->model.getDOFManager().getNonLinearSolver("static").get( "nb_iterations"); if (not(nb_iter > 0 and this->loading_test)) { computeC(el_type, ghost_type); } } OrthotropicParent::computeTangentModuli(el_type, tangent_matrix, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeOrthotropic::beforeSolveStep() { parent::beforeSolveStep(); for (auto & el_type : this->element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { this->nb_state_changes(el_type, _not_ghost).zero(); } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction_inline_impl.hh b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction_inline_impl.hh index 00f3d049a..dac7b184d 100644 --- a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction_inline_impl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction_inline_impl.hh @@ -1,192 +1,191 @@ /** * @file material_iterative_stiffness_reduction.cc * @author Aurelia Isabel Cuba Ramos * @date Thu Feb 18 16:03:56 2016 * * @brief Implementation of material iterative stiffness reduction * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "material_iterative_stiffness_reduction.hh" #include "material_viscoelastic_maxwell.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialIterativeStiffnessReduction:: MaterialIterativeStiffnessReduction(SolidMechanicsModel & model, const ID & id) : IsotropicParent(model, id), eps_u("eps_u", *this), Sc_init("Sc_init", *this), D("D", *this), Gf(0.), crack_band_width(0.), reduction_constant(0.) { AKANTU_DEBUG_IN(); this->registerParam("Gf", Gf, _pat_parsable | _pat_modifiable, "fracture energy"); this->registerParam("crack_band_width", crack_band_width, _pat_parsable | _pat_modifiable, "crack_band_width"); this->registerParam("reduction_constant", reduction_constant, 2., _pat_parsable | _pat_modifiable, "reduction constant"); this->eps_u.initialize(1); this->Sc_init.initialize(1); this->D.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIterativeStiffnessReduction::initMaterial() { AKANTU_DEBUG_IN(); IsotropicParent::initMaterial(); this->Sc_init.copy(this->Sc); for (auto ghost_type : ghost_types) { /// loop over all types in the filter for (auto & el_type : this->element_filter.elementTypes(_ghost_type = ghost_type)) { /// get the stiffness on each quad point auto Sc_it = this->Sc(el_type, ghost_type).begin(); /// get the tangent of the tensile softening on each quad point auto D_it = this->D(el_type, ghost_type).begin(); auto D_end = this->D(el_type, ghost_type).end(); /// get the ultimate strain on each quad auto eps_u_it = this->eps_u(el_type, ghost_type).begin(); // compute the tangent and the ultimate strain for each quad for (; D_it != D_end; ++Sc_it, ++D_it, ++eps_u_it) { *eps_u_it = ((2. * this->Gf) / (*Sc_it * this->crack_band_width)); *D_it = *(Sc_it) / ((*eps_u_it) - ((*Sc_it) / this->E)); } } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- - */ + +/* -------------------------------------------------------------------------- */ template UInt MaterialIterativeStiffnessReduction::updateDamage() { UInt nb_damaged_elements = 0; if (this->norm_max_equivalent_stress >= 1.) { AKANTU_DEBUG_IN(); /// update the damage only on non-ghosts elements! Doesn't make sense to /// update on ghost. GhostType ghost_type = _not_ghost; for (auto && el_type : this->element_filter.elementTypes(spatial_dimension, ghost_type)) { /// loop over all the elements /// get iterators on the needed internal fields auto equivalent_stress_it = this->equivalent_stress(el_type, ghost_type).begin(); auto equivalent_stress_end = this->equivalent_stress(el_type, ghost_type).end(); auto dam_it = this->damage(el_type, ghost_type).begin(); auto reduction_it = this->reduction_step(el_type, ghost_type).begin(); auto eps_u_it = this->eps_u(el_type, ghost_type).begin(); auto Sc_it = this->Sc(el_type, ghost_type).begin(); auto Sc_init_it = this->Sc_init(el_type, ghost_type).begin(); auto D_it = this->D(el_type, ghost_type).begin(); /// EXPERIMENTAL Real E_ef_sum = this->E; Real E_sum = E_ef_sum; if (aka::is_of_type>( *this)) { E_ef_sum = this->getParam("Einf"); E_sum = this->getParam("Einf"); const Vector & Eta = this->getParam("Eta"); const Vector & Ev = this->getParam("Ev"); for (UInt k = 0; k < Eta.size(); ++k) { E_sum += Ev(k); Real E_ef; Real exp_dt_lambda; dynamic_cast &>(*this) .computeEffectiveModulus(k, E_ef, exp_dt_lambda); E_ef_sum += E_ef; } } /// loop over all the quads of the given element type for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it, ++reduction_it, ++eps_u_it, ++Sc_it, ++Sc_init_it, ++D_it) { if (aka::is_of_type>( *this)) { *eps_u_it = 2. * this->Gf / (*Sc_init_it * this->crack_band_width) + *Sc_init_it * (1 / E_ef_sum - 1 / E_sum); *D_it = *(Sc_init_it) / ((*eps_u_it) - ((*Sc_init_it) / E_ef_sum)); } /// check if damage occurs if (*equivalent_stress_it >= (1 - this->dam_tolerance) * this->norm_max_equivalent_stress) { /// check if this element can still be damaged if (*reduction_it == this->max_reductions) { continue; } /// increment the counter of stiffness reduction steps *reduction_it += 1; if (*reduction_it == this->max_reductions) { *dam_it = this->max_damage; } else { /// update the damage on this quad *dam_it = 1. - (1. / std::pow(this->reduction_constant, *reduction_it)); /// update the stiffness on this quad *Sc_it = (*eps_u_it) * (1. - (*dam_it)) * E_ef_sum * (*D_it) / ((1. - (*dam_it)) * E_ef_sum + (*D_it)); } nb_damaged_elements += 1; } } } } // auto rve_model = dynamic_cast(&this->model); // if (rve_model == NULL) { const auto & comm = this->model.getMesh().getCommunicator(); comm.allReduce(nb_damaged_elements, SynchronizerOperation::_sum); //} AKANTU_DEBUG_OUT(); return nb_damaged_elements; } // namespace akantu -/* -------------------------------------------------------------------------- - */ +/* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 41e6768f8..1e9b00221 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1407 +1,1416 @@ /** * @file material.cc * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common part of the material class * * * Copyright (©) 2010-2018 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 "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this), eigen_grad_u(model.getSpatialDimension(), model.getSpatialDimension(), 0.) { AKANTU_DEBUG_IN(); this->registerParam("eigen_grad_u", eigen_grad_u, _pat_parsable, "EigenGradU"); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), is_init(false), fem(fe_engine), finite_deformation(false), model(model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigengradu("eigen_grad_u", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), piola_kirchhoff_2("piola_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() = default; /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(0.), _pat_parsable | _pat_modifiable, "Density"); registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation", finite_deformation, false, _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false, _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); this->model.registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if (finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if (use_previous_stress) { this->piola_kirchhoff_2.initializeHistory(); } this->green_strain.initialize(spatial_dimension * spatial_dimension); } if (use_previous_stress) { this->stress.initializeHistory(); } if (use_previous_gradu) { this->gradu.initializeHistory(); } this->resizeInternals(); auto dim = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(_element_kind = _ek_regular)) { for (auto eigen_gradu : make_view(eigengradu(type), dim, dim)) { eigen_gradu = eigen_grad_u; } } + + for(auto && ghost_type : ghost_types) { + for (const auto & type : + element_filter.elementTypes(_element_kind = _ek_regular, + _ghost_type = ghost_type)) { + gradu_release(type, ghost_type) = -1; + } + } + is_init = true; updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->saveCurrentValues(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::restorePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->restorePreviousValues(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resetInternalsWithHistory() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) if (pair.second->hasHistory()) { pair.second->reset(); pair.second->resetPreviousValues(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the internal forces by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] ghost_type compute the internal forces for _ghost or _not_ghost * element */ void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { auto & internal_force = const_cast &>(model.getInternalForce()); // Mesh & mesh = fem.getMesh(); for (auto && type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ auto * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); fem.computeBtD(stress(type, ghost_type), *sigma_dphi_dx, type, ghost_type, elem_filter); /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by * @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ auto * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem.integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, type, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model.getDOFManager().assembleElementalArrayLocalArray( *int_sigma_dphi_dx, internal_force, type, ghost_type, -1, elem_filter); delete int_sigma_dphi_dx; } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeGradU(ElementType type, GhostType ghost_type) { auto release = model.getDisplacementRelease(); - if (release == last_displacement_release) { + if (release == gradu_release(type, ghost_type)) { return; } const auto & elem_filter = element_filter(type, ghost_type); auto spatial_dimension = model.getSpatialDimension(); auto & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); - ++gradu_release; - last_displacement_release = release; + + gradu_release(type, ghost_type) = release; } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } this->computeGradU(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ this->computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->StoCauchy<1>(type, ghost_type); break; case 2: this->StoCauchy<2>(type, ghost_type); break; case 3: this->StoCauchy<3>(type, ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::StoCauchy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); auto gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); auto piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); auto stress_it = this->stress(el_type, ghost_type).begin(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; auto F_tensor = gradUToF(grad_u); this->StoCauchy(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { AKANTU_DEBUG_OUT(); return; } // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); this->computeGradU(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt tangent_size = getTangentStiffnessVoigtSize(dim); auto * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->zero(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; auto * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); fem.computeBtDB(*tangent_stiffness_matrix, *bt_d_b, 4, type, ghost_type, elem_filter); delete tangent_stiffness_matrix; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; auto * bt_s_b = new Array(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); auto Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); auto piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { auto & Bt_S_B = *Bt_S_B_it; const auto & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.template mul(B, S); Bt_S_B.template mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * K_e = new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); auto * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->zero(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); auto * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; auto * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); auto grad_u_it = gradu_vect.begin(dim, dim); auto D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); auto D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { const auto & grad_u = *grad_u_it; const auto & D = *D_it; auto & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto * shapesd_filtered = new Array( nb_element, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; auto * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); auto grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); auto grad_u_end = this->gradu(type, ghost_type).end(dim, dim); auto stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { auto & grad_u = *grad_u_it; auto & r = *bt_s_it; auto & S = *stress_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; auto S_vect = Material::stressToVoigt(S); Matrix S_voigt(S_vect.storage(), stress_size, 1); r.template mul(B_tensor, S_voigt); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto * r_e = new Array(nb_element, bt_s_size, "r_e"); fem.integrate(*bt_s, *r_e, bt_s_size, type, ghost_type, elem_filter); delete bt_s; model.getDOFManager().assembleElementalArrayLocalArray( *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType /*unused*/) { AKANTU_DEBUG_IN(); AKANTU_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") { return getPotentialEnergy(); } AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") { return getPotentialEnergy(type, index); } AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison{type, 0, ghost_type}; const Array> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + static_cast(is_second_element) * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & /*vect_id*/, ElementType /*type*/, GhostType /*ghost_type*/) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(const ID & /*vect_id*/, ElementType /*type*/, GhostType /*ghost_type*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) const { std::stringstream sstr; std::string ghost_id; if (ghost_type == _ghost) { ghost_id = ":ghost"; } sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) { std::stringstream sstr; std::string ghost_id; if (ghost_type == _ghost) { ghost_id = ":ghost"; } sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) const { std::stringstream sstr; std::string ghost_id; if (ghost_type == _ghost) { ghost_id = ":ghost"; } sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, ElementType type, GhostType ghost_type) { std::stringstream sstr; std::string ghost_id; if (ghost_type == _ghost) { ghost_id = ":ghost"; } sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template const InternalField & Material::getInternal([[gnu::unused]] const ID & int_id) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal([[gnu::unused]] const ID & int_id) { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model.getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for (; el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Array & mat_loc_num = model.getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if (el_begin == el_end) { return; } ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (const auto & type : element_filter.elementTypes( _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model.getMaterialLocalNumbering(type, ghost_type); if (!material_local_new_numbering.exists(type, ghost_type)) { material_local_new_numbering.alloc(elem_filter.size(), 1, type, ghost_type); } Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.size(); Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if (std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->resize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, [[gnu::unused]] const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); auto el_begin = element_list.begin(); auto el_end = element_list.end(); for (auto && gt : ghost_types) { for (auto && type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { if (not element_filter.exists(type, gt) || element_filter(type, gt).empty()) { continue; } auto & elem_filter = element_filter(type, gt); auto & mat_indexes = this->model.getMaterialByElement(type, gt); auto & mat_loc_num = this->model.getMaterialLocalNumbering(type, gt); auto nb_element = this->model.getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if (!material_local_new_numbering.exists(type, gt)) { material_local_new_numbering.alloc(elem_filter.size(), 1, type, gt); } auto & mat_renumbering = material_local_new_numbering(type, gt); const auto & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el{type, 0, gt}; for (UInt i = 0; i < elem_filter.size(); ++i) { el.element = elem_filter(i); if (std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->removeIntegrationPoints(material_local_new_numbering); } } /* -------------------------------------------------------------------------- */ void Material::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::afterSolveStep(bool converged) { if (not converged) { this->restorePreviousState(); return; } for (const auto & type : element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergies(type); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { for (const auto & type : element_filter.elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergiesAfterDamage(type); } } /* -------------------------------------------------------------------------- */ void Material::onDump() { if (this->isFiniteDeformation()) { this->computeAllCauchyStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, [[gnu::unused]] const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (element_filter(type, ghost_type).empty()) { continue; } for (auto & eigengradu : make_view(this->eigengradu(type, ghost_type), spatial_dimension, spatial_dimension)) { eigengradu = prescribed_eigen_grad_u; } } } /* -------------------------------------------------------------------------- */ MaterialFactory & Material::getFactory() { return MaterialFactory::getInstance(); } } // namespace akantu diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 75fe95518..cf2351a8d 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,728 +1,728 @@ /** * @file material.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Feb 21 2018 * * @brief Mother class for all materials * * * Copyright (©) 2010-2018 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_factory.hh" #include "aka_voigthelper.hh" #include "data_accessor.hh" #include "integration_point.hh" #include "parsable.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #include "mesh_events.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_HH_ #define AKANTU_MATERIAL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class Material; } // namespace akantu namespace akantu { using MaterialFactory = Factory; /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = * ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(ElementType el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler, protected SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(const Material & mat) = delete; Material & operator=(const Material & mat) = delete; /// Initialize material with defaults Material(SolidMechanicsModel & model, const ID & id = ""); /// Initialize material with custom mesh & fe_engine Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor ~Material() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Function that materials can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: void computeGradU(ElementType type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeStress(ElementType /* el_type */, GhostType /* ghost_type */ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentModuli(ElementType /*el_type*/, Array & /*tangent_matrix*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(ElementType /*type*/, UInt /*index*/, Vector & /*epot_on_quad_points*/) { AKANTU_TO_IMPLEMENT(); } virtual void updateEnergies(ElementType /*el_type*/) {} virtual void updateEnergiesAfterDamage(ElementType /*el_type*/) {} /// set the material to steady state (to be implemented for materials that /// need it) virtual void setToSteadyState(ElementType /*el_type*/, GhostType /*ghost_type*/ = _not_ghost) {} /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} public: /// extrapolate internal values virtual void extrapolateInternal(const ID & id, const Element & element, const Matrix & points, Matrix & extrapolated); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(const Element & /*element*/) const { AKANTU_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(const Element & /*element*/) const { AKANTU_TO_IMPLEMENT(); } /// get a material celerity to compute the stable time step (default: is the /// push wave speed) virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(InternalField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } template void unregisterInternal(InternalField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material // virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleInternalForces(GhostType ghost_type); /// save the stress in the previous_stress if needed virtual void savePreviousState(); /// restore the stress from previous_stress if needed virtual void restorePreviousState(); /// resets all internals with history and their previous values too virtual void resetInternalsWithHistory(); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); // virtual void // computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// add an element to the local mesh filter inline UInt addElement(ElementType type, UInt element, GhostType ghost_type); inline UInt addElement(const Element & element); /// add many elements at once void addElements(const Array & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ void interpolateStress(ElementTypeMapArray & result, GhostType ghost_type = _not_ghost); /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points and store the * results per facet */ void interpolateStressOnFacets(ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates); /* ------------------------------------------------------------------------ */ /* Common part */ /* ------------------------------------------------------------------------ */ protected: /* ------------------------------------------------------------------------ */ static inline UInt getTangentStiffnessVoigtSize(UInt dim); /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays virtual void resizeInternals(); /* ------------------------------------------------------------------------ */ /* Finite deformation functions */ /* This functions area implementing what is described in the paper of Bathe */ /* et al, in IJNME, Finite Element Formulations For Large Deformation */ /* Dynamic Analysis, Vol 9, 353-386, 1975 */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleInternalForces(GhostType ghost_type); template void computeAllStressesFromTangentModuli(ElementType type, GhostType ghost_type); template void assembleStiffnessMatrix(ElementType type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type); template void assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ public: /// Size of the Stress matrix for the case of finite deformation see: Bathe et /// al, IJNME, Vol 9, 353-386, 1975 static inline UInt getCauchyStressMatrixSize(UInt dim); /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, /// 1975 template static inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & sigma); /// write the stress tensor in the Voigt notation. template static inline decltype(auto) stressToVoigt(const Matrix & stress) { return VoigtHelper::matrixToVoigt(stress); } /// write the strain tensor in the Voigt notation. template static inline decltype(auto) strainToVoigt(const Matrix & strain) { return VoigtHelper::matrixToVoigtWithFactors(strain); } /// write a voigt vector to stress template static inline void voigtToStress(const Vector & voigt, Matrix & stress) { VoigtHelper::voigtToMatrix(voigt, stress); } /// Computation of Cauchy stress tensor in the case of finite deformation from /// the 2nd Piola-Kirchhoff for a given element type template void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost); /// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation /// gradient template inline void StoCauchy(const Matrix & F, const Matrix & S, Matrix & sigma, const Real & C33 = 1.0) const; template static inline void gradUToF(const Matrix & grad_u, Matrix & F); template static inline decltype(auto) gradUToF(const Matrix & grad_u); static inline void rightCauchy(const Matrix & F, Matrix & C); static inline void leftCauchy(const Matrix & F, Matrix & B); template static inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); template static inline decltype(auto) gradUToEpsilon(const Matrix & grad_u); template static inline void gradUToE(const Matrix & grad_u, Matrix & epsilon); template static inline decltype(auto) gradUToE(const Matrix & grad_u); static inline Real stressToVonMises(const Matrix & stress); protected: /// converts global element to local element inline Element convertToLocalElement(const Element & global_element) const; /// converts local element to global element inline Element convertToGlobalElement(const Element & local_element) const; /// converts global quadrature point to local quadrature point inline IntegrationPoint convertToLocalPoint(const IntegrationPoint & global_point) const; /// converts local quadrature point to global quadrature point inline IntegrationPoint convertToGlobalPoint(const IntegrationPoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ void onNodesAdded(const Array & /*unused*/, const NewNodesEvent & /*unused*/) override{}; void onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) override{}; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void beforeSolveStep(); virtual void afterSolveStep(bool converged = true); void onDamageIteration() override; void onDamageUpdate() override; void onDump() override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// return the potential energy for the subset of elements contained by the /// material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained /// by the material virtual Real getEnergy(const std::string & type); /// return the energy (identified by id) for the provided element virtual Real getEnergy(const std::string & energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(FEEngine, fem, FEEngine &); bool isNonLocal() const { return is_non_local; } template const Array & getArray(const ID & id, ElementType type, GhostType ghost_type = _not_ghost) const; template Array & getArray(const ID & id, ElementType type, GhostType ghost_type = _not_ghost); template const InternalField & getInternal(const ID & id) const; template InternalField & getInternal(const ID & id); template inline bool isInternal(const ID & id, ElementKind element_kind) const; template ElementTypeMap getInternalDataPerElem(const ID & id, ElementKind element_kind) const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; /// apply a constant eigengrad_u everywhere in the material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, GhostType /*ghost_type*/ = _not_ghost); bool hasMatrixChanged(const ID & id) { if (id == "K") { return hasStiffnessMatrixChanged() or finite_deformation; } return true; } MatrixType getMatrixType(const ID & id) { if (id == "K") { return getTangentType(); } if (id == "M") { return _symmetric; } return _mt_not_defined; } /// specify if the matrix need to be recomputed for this material virtual bool hasStiffnessMatrixChanged() { return true; } /// specify the type of matrix, if not overloaded the material is not valid /// for static or implicit computations virtual MatrixType getTangentType() { return _mt_not_defined; } /// static method to reteive the material factory static MaterialFactory & getFactory(); /// retrieve maps of internal fields std::map *> getInternalVectorsReal() const { return internal_vectors_real; } std::map *> getInternalVectorsUInt() const { return internal_vectors_uint; } std::map *> getInternalVectorsBool() const { return internal_vectors_bool; } protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; std::map *> internal_vectors_bool; protected: /// Link to the fem object in the model FEEngine & fem; /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel & model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ElementTypeMapArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// eigengrad_u arrays ordered by element types InternalField eigengradu; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types /// (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// tell if the material need the previous strain state bool use_previous_gradu; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; /// vector that contains the names of all the internals that need to /// be transferred when material interfaces move std::vector internals_to_transfer; - Int last_displacement_release{-1}; - Int gradu_release{0}; + /// release of grad u stored per type + ElementTypeMap gradu_release; private: /// eigen_grad_u for the parser Matrix eigen_grad_u; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "material_inline_impl.hh" #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ /// This can be used to automatically write the loop on quadrature points in /// functions such as computeStress. This macro in addition to write the loop /// provides two tensors (matrices) sigma and grad_u #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ auto && grad_u_view = \ make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \ this->spatial_dimension); \ \ auto stress_view = \ make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \ this->spatial_dimension); \ \ if (this->isFiniteDeformation()) { \ stress_view = make_view(this->piola_kirchhoff_2(el_type, ghost_type), \ this->spatial_dimension, this->spatial_dimension); \ } \ \ for (auto && data : zip(grad_u_view, stress_view)) { \ [[gnu::unused]] Matrix & grad_u = std::get<0>(data); \ [[gnu::unused]] Matrix & sigma = std::get<1>(data) #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END } /// This can be used to automatically write the loop on quadrature points in /// functions such as computeTangentModuli. This macro in addition to write the /// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix /// where the elemental tangent moduli should be stored in Voigt Notation #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ auto && grad_u_view = \ make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \ this->spatial_dimension); \ \ auto && stress_view = \ make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \ this->spatial_dimension); \ \ auto tangent_size = \ this->getTangentStiffnessVoigtSize(this->spatial_dimension); \ \ auto && tangent_view = make_view(tangent_mat, tangent_size, tangent_size); \ \ for (auto && data : zip(grad_u_view, stress_view, tangent_view)) { \ [[gnu::unused]] Matrix & grad_u = std::get<0>(data); \ [[gnu::unused]] Matrix & sigma = std::get<1>(data); \ Matrix & tangent = std::get<2>(data); #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END } /* -------------------------------------------------------------------------- */ #define INSTANTIATE_MATERIAL_ONLY(mat_name) \ template class mat_name<1>; /* NOLINT */ \ template class mat_name<2>; /* NOLINT */ \ template class mat_name<3> /* NOLINT */ #define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \ [](UInt dim, const ID &, SolidMechanicsModel & model, \ const ID & id) /* NOLINT */ \ -> std::unique_ptr< \ Material> { /* NOLINT */ \ switch (dim) { \ case 1: \ return std::make_unique>(/* NOLINT */ \ model, id); \ case 2: \ return std::make_unique>(/* NOLINT */ \ model, id); \ case 3: \ return std::make_unique>(/* NOLINT */ \ model, id); \ default: \ AKANTU_EXCEPTION( \ "The dimension " \ << dim \ << "is not a valid dimension for the material " \ << #id); \ } \ } #define INSTANTIATE_MATERIAL(id, mat_name) \ INSTANTIATE_MATERIAL_ONLY(mat_name); \ static bool material_is_alocated_##id [[gnu::unused]] = \ MaterialFactory::getInstance().registerAllocator( \ #id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name)) #endif /* AKANTU_MATERIAL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index e24796600..d8e032a88 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1243 +1,1246 @@ /** * @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: Wed Feb 21 2018 * * @brief Implementation of the SolidMechanicsModel class * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * 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 * @param memory_id the id of the memory * @param model_type this is an internal parameter for inheritance purposes */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif material_selector = std::make_shared(material_index), this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param options * \parblock * contains the different options to initialize the model * \li \c analysis_method specify the type of solver to use * \endparblock */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { material_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the materials if (not this->parser.getLastParsedFile().empty()) { this->instantiateMaterials(); this->initMaterials(); } this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ /** * 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::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForces(); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") { return _mt_not_defined; } if (matrix_id == "K") { auto matrix_type = _unsymmetric; for (auto & material : materials) { matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id)); } } return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) { material->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep(bool converged) { for (auto & material : materials) { material->afterSolveStep(converged); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as \f$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->zero(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) { this->non_local_manager->computeAllNonLocalStresses(); } // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasMatrixChanged("K"); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").zero(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) { return; } this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) { continue; } for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const auto & v = *v_it; const auto & m = *m_it; Real mv2 = 0.; auto is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); auto count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) { mv2 += v(i) * v(i) * m(i); } } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension), make_view(*this->velocity, spatial_dimension))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } } else { AKANTU_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); auto vit = vel_on_quad.begin(Model::spatial_dimension); auto vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) { work -= int_force(i) * incr_or_velo(i); } else { work += ext_force(i) * incr_or_velo(i); } } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) { work *= this->getTimeStep(); } AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) { energy += material->getEnergy(energy_id); } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->material_index.initialize(mesh, _element_kind = _ek_not_defined, _spatial_dimension = spatial_dimension, _with_nb_element = true, _default_value = UInt(-1)); this->material_local_numbering.initialize( mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (const auto & elem : element_list) { if (mesh.getSpatialDimension(elem.type) != spatial_dimension) { continue; } 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 fails in parallel if the event is sent on facet between constructor // and initFull \todo: to debug... this->assignMaterialToElements(&filter); - for (auto & material : materials) { - material->onElementsAdded(element_list, event); - } + splitByMaterial(element_list, [&event](auto && mat, auto && elements) { + mat.onElementsAdded(elements, event); + }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) { mass->resize(nb_nodes, 0.); } if (velocity) { velocity->resize(nb_nodes, 0.); } if (acceleration) { acceleration->resize(nb_nodes, 0.); } if (external_force) { external_force->resize(nb_nodes, 0.); } if (internal_force) { internal_force->resize(nb_nodes, 0.); } if (blocked_dofs) { blocked_dofs->resize(nb_nodes, false); } if (current_position) { current_position->resize(nb_nodes, 0.); } if (previous_displacement) { previous_displacement->resize(nb_nodes, 0.); } if (displacement_increment) { displacement_increment->resize(nb_nodes, 0.); } for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array & /*element_list*/, const Array & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) { mesh.removeNodesFromArray(*mass, new_numbering); } if (velocity) { mesh.removeNodesFromArray(*velocity, new_numbering); } if (acceleration) { mesh.removeNodesFromArray(*acceleration, new_numbering); } if (internal_force) { mesh.removeNodesFromArray(*internal_force, new_numbering); } if (external_force) { mesh.removeNodesFromArray(*external_force, new_numbering); } if (blocked_dofs) { mesh.removeNodesFromArray(*blocked_dofs, new_numbering); } // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) { mesh.removeNodesFromArray(*displacement_increment, new_numbering); } if (previous_displacement) { mesh.removeNodesFromArray(*previous_displacement, new_numbering); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) { velocity->printself(stream, indent + 2); } if (acceleration) { acceleration->printself(stream, indent + 2); } if (mass) { mass->printself(stream, indent + 2); } external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + materials [" << std::endl; for (const auto & material : materials) { material->printself(stream, indent + 2); } stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeNonLocal() { this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( GhostType ghost_type) { for (auto & mat : materials) { MaterialNonLocalInterface * mat_non_local; if ((mat_non_local = dynamic_cast(mat.get())) == nullptr) { continue; } ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (const auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } mat_non_local->initMaterialNonLocal(); mat_non_local->insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) { for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) { material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { if (not aka::is_of_type(*mat)) { continue; } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return getFEEngineClassBoundary(name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::splitElementByMaterial( const Array & elements, std::vector> & elements_per_mat) const { for (const auto & el : elements) { + if(Mesh::getSpatialDimension(el.type) != spatial_dimension) { + continue; + } Element mat_el = el; mat_el.element = this->material_local_numbering(el); elements_per_mat[this->material_index(el)].push_back(mat_el); } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_material_id: { size += elements.size() * sizeof(UInt); break; } case SynchronizationTag::_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case SynchronizationTag::_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case SynchronizationTag::_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case SynchronizationTag::_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { packElementalDataHelper( material_index, buffer, elements, false, getFEEngine()); break; } case SynchronizationTag::_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_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 SynchronizationTag::_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { for (auto && element : elements) { UInt recv_mat_index; buffer >> recv_mat_index; UInt & mat_index = material_index(element); if (mat_index != UInt(-1)) { continue; } // add ghosts element to the correct material mat_index = recv_mat_index; UInt index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; } break; } case SynchronizationTag::_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_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 SynchronizationTag::_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: { } } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case SynchronizationTag::_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case SynchronizationTag::_smm_res: /* FALLTHRU */ case SynchronizationTag::_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 301bc4f27..671601fb7 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,583 +1,584 @@ /** * @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: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * Copyright (©) 2010-2018 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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, 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; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// 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, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database 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 virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* 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(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* 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, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; 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); void dump() override; 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, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// 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 array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical 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(); } /// 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, ElementType type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, 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); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } - AKANTU_GET_MACRO(DisplacementRelease, displacement_release, Int); + AKANTU_GET_MACRO_NOT_CONST(DisplacementRelease, displacement_release, Int &); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array - UInt displacement_release{0}; + Int displacement_release{0}; /// release version of the current_position array - UInt current_position_release{0}; + Int current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; + /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// 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; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc index 0806dd7aa..6fec34b19 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc @@ -1,214 +1,198 @@ /** * @file material_cohesive_bilinear.cc * * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Feb 20 2018 * * @brief Bilinear cohesive constitutive law * * * Copyright (©) 2010-2018 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 "material_cohesive_bilinear.hh" //#include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveBilinear::MaterialCohesiveBilinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model, id) { AKANTU_DEBUG_IN(); this->registerParam("delta_0", delta_0, Real(0.), _pat_parsable | _pat_readable, "Elastic limit displacement"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::initMaterial() { AKANTU_DEBUG_IN(); this->sigma_c_eff.setRandomDistribution(this->sigma_c.getRandomParameter()); MaterialCohesiveLinear::initMaterial(); this->delta_max.setDefaultValue(delta_0); this->insertion_stress.setDefaultValue(0); this->delta_max.reset(); this->insertion_stress.reset(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::onElementsAdded(element_list, event); bool scale_traction = false; // don't scale sigma_c if volume_s hasn't been specified by the user if (!Math::are_float_equal(this->volume_s, 0.)) { scale_traction = true; } - for (auto & el: element_list) { - // filter not ghost cohesive elements - if (el.ghost_type != _not_ghost || - Mesh::getKind(el.type) != _ek_cohesive) { - continue; - } - - ElementType type = el.type; - auto & elem_filter = this->element_filter(type, el.ghost_type); - auto filter_begin = elem_filter.begin(); - auto filter_end = elem_filter.end(); - - UInt index = el.element; - if (std::find(filter_begin, filter_end, index) == filter_end) { - continue; - } - - UInt nb_quad_per_element = this->fem_cohesive.getNbIntegrationPoints(type); - auto sigma_c_begin = make_view(this->sigma_c_eff(type), nb_quad_per_element).begin(); - Vector sigma_c_vec = sigma_c_begin[index]; + for (const auto & el: element_list) { + Vector sigma_c_vec = this->sigma_c_eff.get(el); + Vector delta_c_vec = this->delta_c_eff.get(el); - auto delta_c_begin = make_view(this->delta_c_eff(type), nb_quad_per_element).begin(); - Vector delta_c_vec = delta_c_begin[index]; + auto global_el = el; + global_el.element = this->element_filter(el); if (scale_traction) { - scaleTraction(el, sigma_c_vec); + scaleTraction(global_el, sigma_c_vec); } /** * Recompute sigma_c as * @f$ {\sigma_c}_\textup{new} = * \frac{{\sigma_c}_\textup{old} \delta_c} {\delta_c - \delta_0} @f$ */ + for(auto && data : zip(delta_c_vec, sigma_c_vec)) { + auto & delta_c = std::get<0>(data); + auto & sigma_c = std::get<1>(data); + delta_c = 2 * this->G_c / sigma_c; - for (UInt q = 0; q < nb_quad_per_element; ++q) { - delta_c_vec(q) = 2 * this->G_c / sigma_c_vec(q); - - if (delta_c_vec(q) - delta_0 < Math::getTolerance()) { + if (delta_c - delta_0 < Math::getTolerance()) { AKANTU_ERROR("delta_0 = " << delta_0 << " must be lower than delta_c = " - << delta_c_vec(q) + << delta_c << ", modify your material file"); } - sigma_c_vec(q) *= delta_c_vec(q) / (delta_c_vec(q) - delta_0); + sigma_c *= delta_c / (delta_c - delta_0); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::scaleTraction( const Element & el, Vector & sigma_c_vec) { AKANTU_DEBUG_IN(); Real base_sigma_c = this->sigma_c_eff; const Mesh & mesh_facets = this->model->getMeshFacets(); const FEEngine & fe_engine = this->model->getFEEngine(); auto coh_element_to_facet_begin = mesh_facets.getSubelementToElement(el.type).begin(2); const Vector & coh_element_to_facet = coh_element_to_facet_begin[el.element]; // compute bounding volume Real volume = 0; // loop over facets for (UInt f = 0; f < 2; ++f) { const Element & facet = coh_element_to_facet(f); const Array> & facet_to_element = mesh_facets.getElementToSubelement(facet.type, facet.ghost_type); const std::vector & element_list = facet_to_element(facet.element); auto elem = element_list.begin(); auto elem_end = element_list.end(); // loop over elements connected to each facet for (; elem != elem_end; ++elem) { // skip cohesive elements and dummy elements if (*elem == ElementNull || Mesh::getKind(elem->type) == _ek_cohesive) { continue; } // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } } // scale sigma_c sigma_c_vec -= base_sigma_c; sigma_c_vec *= std::pow(this->volume_s / volume, 1. / this->m_s); sigma_c_vec += base_sigma_c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::computeTraction(normal, el_type, ghost_type); // adjust damage - auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); - auto delta_max_it = this->delta_max(el_type, ghost_type).begin(); - auto damage_it = this->damage(el_type, ghost_type).begin(); - auto damage_end = this->damage(el_type, ghost_type).end(); - - for (; damage_it != damage_end; ++damage_it, ++delta_max_it, ++delta_c_it) { - *damage_it = - std::max((*delta_max_it - delta_0) / (*delta_c_it - delta_0), Real(0.)); - *damage_it = std::min(*damage_it, Real(1.)); + for(auto && data : zip(this->delta_c_eff(el_type, ghost_type), + this->delta_max(el_type, ghost_type), + this->damage(el_type, ghost_type))) { + const auto & delta_c = std::get<0>(data); + const auto & delta_max = std::get<1>(data); + auto & dam = std::get<2>(data); + + dam = std::max((delta_max - delta_0) / (delta_c - delta_0), Real(0.)); + dam = std::min(dam, Real(1.)); } } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_bilinear, MaterialCohesiveBilinear); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc index 371326f6b..cf0c3bd51 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc @@ -1,445 +1,445 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Wed Feb 21 2018 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * * Copyright (©) 2010-2018 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 "dof_synchronizer.hh" #include "material_cohesive_linear.hh" +#include "dof_synchronizer.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam( "contact_after_breaking", contact_after_breaking, false, _pat_parsable | _pat_readable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high " "enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); if (not Math::are_float_equal(delta_c, 0.)) { delta_c_eff.setDefaultValue(delta_c); } else { delta_c_eff.setDefaultValue(2 * G_c / sigma_c); } if (model->getIsExtrinsic()) { scaleInsertionTraction(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::updateInternalParameters() { /// compute scalars beta2_kappa2 = beta * beta / kappa / kappa; beta2_kappa = beta * beta / kappa; if (Math::are_float_equal(beta, 0)) { beta2_inv = 0; } else { beta2_inv = 1. / beta / beta; } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) { return; } const Mesh & mesh_facets = model->getMeshFacets(); const auto & fe_engine = model->getFEEngine(); const auto & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); Real base_sigma_c = sigma_c; for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { const Array> & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.size(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet auto sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; auto elem = element_list.begin(); auto elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) { continue; } // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion( bool check_only) { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const auto & facets_check = inserter.getCheckFacets(type_facet); auto & f_insertion = inserter.getInsertionFacets(type_facet); auto & f_filter = facet_filter(type_facet); auto & sig_c_eff = sigma_c_eff(type_cohesive); auto & del_c = delta_c_eff(type_cohesive); auto & ins_stress = insertion_stress(type_cohesive); auto & trac_old = tractions.previous(type_cohesive); const auto & f_stress = model->getStressOnFacets(type_facet); const auto & sigma_lim = sigma_c(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); #ifndef AKANTU_NDEBUG UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, "The cohesive element and the corresponding facet do " "not have the same numbers of integration points"); #endif UInt nb_facet = f_filter.size(); // if (nb_facet == 0) continue; auto sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const auto & tangents = model->getTangents(type_facet); const auto & normals = model->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); auto normal_begin = normals.begin(spatial_dimension); auto tangent_begin = tangents.begin(tangents.getNbComponent()); auto facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector> new_normal_traction; std::vector new_delta_c; // loop over each facet belonging to this material for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { UInt facet = f_filter(f); // skip facets where check shouldn't be realized if (!facets_check(facet)) { continue; } // compute the effective norm on each quadrature point of the facet for (UInt q = 0; q < nb_quad_facet; ++q) { UInt current_quad = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector normal_traction_vec(normal_traction(q)); // compute normal and effective stress stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, normal_traction_vec); } // verify if the effective stress overcomes the threshold Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) { final_stress = *std::max_element( stress_check.storage(), stress_check.storage() + nb_quad_facet); } if (final_stress > *sigma_lim_it) { f_insertion(facet) = true; if (check_only) { continue; } // store the new cohesive material parameters for each quadrature // point for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = stress_check(q); Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) { normal_traction_vec *= -1.; } new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) { new_delta = 2 * G_c / new_sigma; } else { new_delta = (*sigma_lim_it) / new_sigma * delta_c; } new_delta_c.push_back(new_delta); } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.size(); UInt new_nb_quad_points = new_sigmas.size(); sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); ins_stress.resize(old_nb_quad_points + new_nb_quad_points); trac_old.resize(old_nb_quad_points + new_nb_quad_points); del_c.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); auto normal_opening_norm_it = normal_opening_norm(el_type, ghost_type).begin(); auto tangential_opening_norm_it = tangential_opening_norm(el_type, ghost_type).begin(); auto normal_it = normal.begin(spatial_dimension); auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_max_it = delta_max(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (not this->model->isDefaultSolverExplicit()) { this->delta_max(el_type, ghost_type) .copy(this->delta_max.previous(el_type, ghost_type)); } /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++normal_opening_norm_it, ++tangential_opening_norm_it) { bool penetration{false}; this->computeTractionOnQuad( *traction_it, *opening_it, *normal_it, *delta_max_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, normal_opening, tangential_opening, *normal_opening_norm_it, *tangential_opening_norm_it, *damage_it, penetration, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction( ElementType el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto normal_opening_norm_it = normal_opening_norm(el_type, ghost_type).begin(); auto tangential_opening_norm_it = tangential_opening_norm(el_type, ghost_type).begin(); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++normal_opening_norm_it, ++tangential_opening_norm_it) { // Real normal_opening_norm{0}; // Real tangential_opening_norm{0}; bool penetration{false}; this->computeTangentTractionOnQuad( *tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, *normal_opening_norm_it, *tangential_opening_norm_it, *damage_it, penetration, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear, MaterialCohesiveLinear); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc index 15842f003..476abf9c2 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc @@ -1,585 +1,577 @@ /** * @file material_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Mon Feb 19 2018 * * @brief Specialization of the material class for cohesive elements * * * Copyright (©) 2010-2018 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 "material_cohesive.hh" #include "aka_random_generator.hh" #include "dof_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive( model.getFEEngineClass("CohesiveFEEngine")), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), eigen_opening("eigen_opening", *this), tractions("tractions", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), normal_opening_norm("normal_opening_norm", *this), tangential_opening_norm("tangential_opening_norm", *this), delta_max("delta_max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normals("normals", *this) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->element_filter.initialize(this->model->getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_cohesive); - // this->model->getMesh().initElementTypeMapArray( - // this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) { this->facet_filter.initialize(this->model->getMeshFacets(), _spatial_dimension = spatial_dimension - 1, _element_kind = _ek_regular); } // this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, // spatial_dimension - // 1); this->reversible_energy.initialize(1); this->total_energy.initialize(1); this->tractions.initialize(spatial_dimension); this->tractions.initializeHistory(); this->contact_tractions.initialize(spatial_dimension); this->contact_opening.initialize(spatial_dimension); this->normal_opening_norm.initialize(1); this->normal_opening_norm.initializeHistory(); this->tangential_opening_norm.initialize(1); this->tangential_opening_norm.initializeHistory(); this->opening.initialize(spatial_dimension); this->opening.initializeHistory(); this->normals.initialize(spatial_dimension); this->eigen_opening.initialize(spatial_dimension); this->eigen_opening.setDefaultValue(0.); this->delta_max.initialize(1); this->damage.initialize(1); if (this->model->getIsExtrinsic()) { this->sigma_c.initialize(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() = default; /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) { this->delta_max.initializeHistory(); } if (this->use_previous_opening) { this->opening.initializeHistory(); this->normal_opening_norm.initializeHistory(); this->tangential_opening_norm.initializeHistory(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto & traction = tractions(type, ghost_type); auto & contact_traction = contact_tractions(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ auto * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto contact_traction_it = contact_traction.begin(spatial_dimension, 1); auto shapes_filtered_begin = shapes.begin(1, size_of_shapes); auto traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ auto * partial_int_t_N = new Array( nb_element, spatial_dimension * size_of_shapes, "int_t_N"); fem_cohesive.integrate(*traction_cpy, *partial_int_t_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); delete traction_cpy; auto * int_t_N = new Array( nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); Real * int_t_N_val = int_t_N->storage(); Real * partial_int_t_N_val = partial_int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val); std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val + size_of_shapes * spatial_dimension); for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) { int_t_N_val[n] *= -1.; } int_t_N_val += nb_nodes_per_element * spatial_dimension; partial_int_t_N_val += size_of_shapes * spatial_dimension; } delete partial_int_t_N; /// assemble model->getDOFManager().assembleElementalArrayLocalArray( *int_t_N, internal_force, type, ghost_type, 1, elem_filter); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & shapes = fem_cohesive.getShapes(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0U) { continue; } UInt size_of_shapes = shapes.getNbComponent(); auto * shapes_filtered = new Array(nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { auto * shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension * size_of_shapes) = -1; } /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} /// @f$ auto * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); computeNormal(model->getMesh().getNodes(), normals(type, ghost_type), type, ghost_type); // computeNormal(model->getCurrentPosition(), normals(type, ghost_type), // type, // ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix->zero(); computeTangentTraction(type, *tangent_stiffness_matrix, normals(type, ghost_type), ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; auto * at_nt_d_n_a = new Array(nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator> shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.zero(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt n = 0; n < size_of_shapes; ++n) { N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); } } /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} *{\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; auto * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type, ghost_type, elem_filter); delete at_nt_d_n_a; model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif for (const auto & type : element_filter.elementTypes( spatial_dimension, ghost_type, _ek_cohesive)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getMesh().getNodes(), normals(type, ghost_type), type, ghost_type); // computeNormal(model->getCurrentPosition(), normals(type, ghost_type), // type, // ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normals(type, ghost_type), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); normal.zero(); #define COMPUTE_NORMAL(type) \ fem_cohesive.getShapeFunctions() \ .computeNormalsOnIntegrationPoints( \ position, normal, ghost_type, element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); #define COMPUTE_OPENING(type) \ fem_cohesive.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ displacement, opening, spatial_dimension, ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING opening -= eigen_opening(type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::updateEnergies(ElementType type) { AKANTU_DEBUG_IN(); if (Mesh::getKind(type) != _ek_cohesive) { return; } - Vector b(spatial_dimension); - Vector h(spatial_dimension); - auto erev = reversible_energy(type).begin(); - auto etot = total_energy(type).begin(); - auto traction_it = tractions(type).begin(spatial_dimension); - auto traction_old_it = tractions.previous(type).begin(spatial_dimension); - auto opening_it = opening(type).begin(spatial_dimension); - auto opening_old_it = opening.previous(type).begin(spatial_dimension); - - auto traction_end = tractions(type).end(spatial_dimension); - /// loop on each quadrature point - for (; traction_it != traction_end; ++traction_it, ++traction_old_it, - ++opening_it, ++opening_old_it, ++erev, - ++etot) { - /// trapezoidal integration - b = *opening_it; - b -= *opening_old_it; + for (auto && data : + zip(make_view(tractions(type), spatial_dimension), + make_view(tractions.previous(type), spatial_dimension), + make_view(opening(type), spatial_dimension), + make_view(opening.previous(type), spatial_dimension), + reversible_energy(type), total_energy(type))) { + auto & traction = std::get<0>(data); + auto & opening = std::get<2>(data); + auto & prev_opening = std::get<3>(data); + + auto & erev = std::get<4>(data); + auto & etot = std::get<5>(data); - h = *traction_old_it; - h += *traction_it; + /// trapezoidal integration + auto h = traction + std::get<1>(data); + auto b = opening - prev_opening; + etot += b.dot(h) / 2.; - *etot += .5 * b.dot(h); - *erev = .5 * traction_it->dot(*opening_it); + erev = 1. / 2. * traction.dot(opening); } - /// update old values - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { erev += fem_cohesive.integrate(reversible_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { Array dissipated_energy(total_energy(type, _not_ghost)); dissipated_energy -= reversible_energy(type, _not_ghost); edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { auto & el_filter = element_filter(type, _not_ghost); UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost); UInt nb_quad_points = el_filter.size() * nb_quad_per_el; Array contact_energy(nb_quad_points); auto contact_traction_it = contact_tractions(type, _not_ghost).begin(spatial_dimension); auto contact_opening_it = contact_opening(type, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt q = 0; q < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++q) { contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(const std::string & type) { if (type == "reversible") { return getReversibleEnergy(); } if (type == "dissipated") { return getDissipatedEnergy(); } if (type == "cohesive contact") { return getContactEnergy(); } return 0.; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index f2f3f27cd..56546e8cb 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,724 +1,724 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Fabian Barras * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Solid mechanics model for cohesive elements * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), *model.cohesive_synchronizer) {} /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto *cohesive_nodes_event = dynamic_cast(&nodes_event); if (cohesive_nodes_event == nullptr) { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_flags = mesh_accessor.getNodesFlags(); auto nb_old_nodes = nodes_flags.size(); nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && data : zip(old_nodes, new_nodes)) { UInt old_node; UInt new_node; std::tie(old_node, new_node) = data; nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(elements_event); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } Vector nb_new_stuff = {nb_new_nodes, elements_event.getList().size()}; const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_stuff, SynchronizerOperation::_sum); if (nb_new_stuff(1) > 0) { mesh.sendEvent(elements_event); } if (nb_new_stuff(0) > 0) { mesh.sendEvent(nodes_event); // mesh.sendEvent(global_ids_updater.getChangedNodeEvent()); } return std::make_tuple(nb_new_stuff(0), nb_new_stuff(1)); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->material_selector); this->material_selector = tmp_material_selector; #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) - // this->mesh.getDumper("cohesive elements").setTimeStep(time_step); + this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = aka::as_type(options); this->is_extrinsic = smmc_options.is_extrinsic; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = aka::as_type(mesh_facets.getElementSynchronizer()); // synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, SynchronizationTag::_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (not are_materials_instantiated) { instantiateMaterials(); } /// find the first cohesive material UInt cohesive_index = UInt(-1); for (auto && material : enumerate(materials)) { if (dynamic_cast(std::get<1>(material).get()) != nullptr) { cohesive_index = std::get<0>(material); break; } } if (cohesive_index == UInt(-1)) { AKANTU_EXCEPTION("No cohesive materials in the material input file"); } material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = material_selector->getFallbackValue()); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = (*material_selector)(element); auto & mat = aka::as_type(*materials[mat_index]); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initMaterials(); if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.empty()) { continue; } type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (const auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) { continue; } /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) { continue; } Matrix quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; // for (UInt q = 0; q < nb_quad_per_facet; ++q) { // for (UInt s = 0; s < Model::spatial_dimension; ++s) { // el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + // f * nb_quad_per_facet + q, // s) = quad_f(global_facet * nb_quad_per_facet + q, // s); // } // } //} } } } /// loop over non cohesive materials for (auto && material : materials) { if (aka::is_of_type(material)) { continue; } /// initialize the interpolation function material->initElementalFieldInterpolation(elements_quad_facets); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { auto & mat = aka::as_type(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { if (not aka::is_of_type(material)) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } this->synchronize(SynchronizationTag::_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); for (auto & mat : materials) { if (aka::is_of_type(mat)) { /// check which not ghost cohesive elements are to be created auto * mat_cohesive = aka::as_type(mat.get()); mat_cohesive->checkInsertion(); } } /// communicate data among processors // this->synchronize(SynchronizationTag::_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onElementsAdded(element_list, event); if (is_extrinsic) { resizeFacetStress(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); SolidMechanicsModel::onNodesAdded(new_nodes, event); const CohesiveNewNodesEvent * cohesive_event; if ((cohesive_event = dynamic_cast(&event)) == nullptr) { return; } const auto & old_nodes = cohesive_event->getOldNodesList(); auto copy = [this, &new_nodes, &old_nodes](auto & arr) { UInt new_node; UInt old_node; auto view = make_view(arr, spatial_dimension); auto begin = view.begin(); for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; auto old_ = begin + old_node; auto new_ = begin + new_node; *new_ = *old_; } }; copy(*displacement); copy(*blocked_dofs); if (velocity) { copy(*velocity); } if (acceleration) { copy(*acceleration); } if (current_position) { copy(*current_position); } if (previous_displacement) { copy(*previous_displacement); } // if (external_force) // copy(*external_force); // if (internal_force) // copy(*internal_force); if (displacement_increment) { copy(*displacement_increment); } copy(getDOFManager().getSolution("displacement")); // this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep(bool converged) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ if (converged) { for (auto & mat : materials) { if (mat->isFiniteDeformation()) { mat->computeAllCauchyStresses(_not_ghost); } } } SolidMechanicsModel::afterSolveStep(converged); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << "\n"; SolidMechanicsModel::printself(stream, indent + 2); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); // for (auto && ghost_type : ghost_types) { // for (const const auto & type : // mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { // UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); // UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") // .getNbIntegrationPoints(type, // ghost_type); // UInt new_size = nb_facet * nb_quadrature_points; // facet_stress(type, ghost_type).resize(new_size); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh index 867e8b485..e01716ad5 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_fixture.hh @@ -1,347 +1,346 @@ /** * @file test_cohesive_fixture.hh * * @author Nicolas Richart * * @date creation: Wed Jan 10 2018 * @date last modification: Tue Feb 20 2018 * * @brief Coehsive element test fixture * * * Copyright (©) 2016-2018 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 "communicator.hh" #include "solid_mechanics_model_cohesive.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_TEST_COHESIVE_FIXTURE_HH_ #define AKANTU_TEST_COHESIVE_FIXTURE_HH_ using namespace akantu; template <::akantu::AnalysisMethod t> using analysis_method_t = std::integral_constant<::akantu::AnalysisMethod, t>; class StrainIncrement : public BC::Functor { public: StrainIncrement(const Matrix & strain, BC::Axis dir) : strain_inc(strain), dir(dir) {} void operator()(UInt /*node*/, Vector & flags, Vector & primal, const Vector & coord) const { if (std::abs(coord(dir)) < 1e-8) { return; } flags.set(true); primal += strain_inc * coord; } static const BC::Functor::Type type = BC::Functor::_dirichlet; private: Matrix strain_inc; BC::Axis dir; }; template class TestSMMCFixture : public ::testing::Test { public: static constexpr ElementType cohesive_type = std::tuple_element_t<0, param_>::value; static constexpr ElementType type_1 = std::tuple_element_t<1, param_>::value; static constexpr ElementType type_2 = std::tuple_element_t<2, param_>::value; static constexpr size_t dim = ElementClass::getSpatialDimension(); void SetUp() { mesh = std::make_unique(this->dim); if (Communicator::getWorldCommunicator().whoAmI() == 0) { EXPECT_NO_THROW({ mesh->read(this->mesh_name); }); } mesh->distribute(); } void TearDown() { model.reset(nullptr); mesh.reset(nullptr); } void createModel() { model = std::make_unique(*mesh); model->initFull(_analysis_method = this->analysis_method, _is_extrinsic = this->is_extrinsic); auto time_step = this->model->getStableTimeStep() * 0.01; this->model->setTimeStep(time_step); if (dim == 1) { surface = 1; group_size = 1; return; } auto facet_type = mesh->getFacetType(this->cohesive_type); auto & fe_engine = model->getFEEngineBoundary(); const auto & group = mesh->getElementGroup("insertion"); group_size = group.size(_ghost_type = _not_ghost); const auto & elements = group.getElements(facet_type); Array ones(fe_engine.getNbIntegrationPoints(facet_type) * group_size); ones.set(1.); surface = fe_engine.integrate(ones, facet_type, _not_ghost, elements); - mesh->getCommunicator().allReduce(surface, SynchronizerOperation::_sum); - - + mesh->getCommunicator().allReduce(surface, SynchronizerOperation::_sum); mesh->getCommunicator().allReduce(group_size, SynchronizerOperation::_sum); #define debug_ 0 #if debug_ this->model->addDumpFieldVector("displacement"); this->model->addDumpFieldVector("velocity"); this->model->addDumpFieldVector("internal_force"); this->model->addDumpFieldVector("external_force"); this->model->addDumpField("blocked_dofs"); this->model->addDumpField("stress"); this->model->addDumpField("strain"); this->model->assembleInternalForces(); this->model->setBaseNameToDumper("cohesive elements", "cohesive_elements"); this->model->addDumpFieldVectorToDumper("cohesive elements", "displacement"); this->model->addDumpFieldToDumper("cohesive elements", "damage"); this->model->addDumpFieldToDumper("cohesive elements", "tractions"); this->model->addDumpFieldToDumper("cohesive elements", "opening"); this->model->dump(); this->model->dump("cohesive elements"); #endif } void setInitialCondition(const Matrix & strain) { for (auto && data : zip(make_view(this->mesh->getNodes(), this->dim), make_view(this->model->getDisplacement(), this->dim))) { const auto & pos = std::get<0>(data); auto & disp = std::get<1>(data); disp = strain * pos; } + ++model->getDisplacementRelease(); } bool checkDamaged() { UInt nb_damaged = 0; auto & damage = model->getMaterial("insertion").getArray("damage", cohesive_type); for (auto d : damage) { if (d >= .99) { ++nb_damaged; } } return (nb_damaged == group_size); } void steps(const Matrix & strain) { StrainIncrement functor((1. / 300) * strain, this->dim == 1 ? _x : _y); for (auto _ [[gnu::unused]] : arange(nb_steps)) { this->model->applyBC(functor, "loading"); this->model->applyBC(functor, "fixed"); if (this->is_extrinsic) { this->model->checkCohesiveStress(); } this->model->solveStep(); #if debug_ this->model->dump(); this->model->dump("cohesive elements"); #endif } } void checkInsertion() { auto nb_cohesive_element = this->mesh->getNbElement(cohesive_type); mesh->getCommunicator().allReduce(nb_cohesive_element, SynchronizerOperation::_sum); EXPECT_EQ(nb_cohesive_element, group_size); } void checkDissipated(Real expected_density) { Real edis = this->model->getEnergy("dissipated"); EXPECT_NEAR(this->surface * expected_density, edis, 5e-1); } void testModeI() { this->createModel(); auto & mat_el = this->model->getMaterial("body"); auto speed = mat_el.getPushWaveSpeed(Element()); auto direction = _y; if (dim == 1) { direction = _x; } auto length = mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction); nb_steps = length / speed / model->getTimeStep(); SCOPED_TRACE(std::to_string(this->dim) + "D - " + std::to_string(type_1) + ":" + std::to_string(type_2)); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); Matrix strain; if (dim == 1) { strain = {{1.}}; } else if (dim == 2) { strain = {{-nu, 0.}, {0., 1. - nu}}; strain *= (1. + nu); } else if (dim == 3) { strain = {{-nu, 0., 0.}, {0., 1., 0.}, {0., 0., -nu}}; } strain *= sigma_c / E; this->setInitialCondition((1 - 1e-5) * strain); this->steps(2e-2 * strain); } void testModeII() { this->createModel(); auto & mat_el = this->model->getMaterial("body"); Real speed; try { speed = mat_el.getShearWaveSpeed(Element()); // the slowest speed if exists } catch (...) { speed = mat_el.getPushWaveSpeed(Element()); } auto direction = _y; if (dim == 1) direction = _x; auto length = mesh->getUpperBounds()(direction) - mesh->getLowerBounds()(direction); nb_steps = 2 * length / 2. / speed / model->getTimeStep(); SCOPED_TRACE(std::to_string(this->dim) + "D - " + std::to_string(type_1) + ":" + std::to_string(type_2)); if (this->dim > 1) this->model->applyBC(BC::Dirichlet::FlagOnly(_y), "sides"); if (this->dim > 2) this->model->applyBC(BC::Dirichlet::FlagOnly(_z), "sides"); auto & mat_co = this->model->getMaterial("insertion"); Real sigma_c = mat_co.get("sigma_c"); Real beta = mat_co.get("beta"); // Real G_c = mat_co.get("G_c"); Real E = mat_el.get("E"); Real nu = mat_el.get("nu"); Matrix strain; if (dim == 1) { strain = {{1.}}; } else if (dim == 2) { strain = {{0., 1.}, {0., 0.}}; strain *= (1. + nu); } else if (dim == 3) { strain = {{0., 1., 0.}, {0., 0., 0.}, {0., 0., 0.}}; strain *= (1. + nu); } strain *= 2 * beta * beta * sigma_c / E; // nb_steps *= 5; this->setInitialCondition((1. - 1e-5) * strain); this->steps(0.005 * strain); } protected: std::unique_ptr mesh; std::unique_ptr model; std::string mesh_name{std::to_string(cohesive_type) + std::to_string(type_1) + (type_1 == type_2 ? "" : std::to_string(type_2)) + ".msh"}; bool is_extrinsic; AnalysisMethod analysis_method; Real surface{0}; UInt nb_steps{1000}; UInt group_size{10000}; }; /* -------------------------------------------------------------------------- */ template constexpr ElementType TestSMMCFixture::cohesive_type; template constexpr ElementType TestSMMCFixture::type_1; template constexpr ElementType TestSMMCFixture::type_2; template constexpr size_t TestSMMCFixture::dim; /* -------------------------------------------------------------------------- */ using IsExtrinsicTypes = std::tuple; using AnalysisMethodTypes = std::tuple>; using coh_types = gtest_list_t, std::tuple<_element_type_cohesive_2d_4, _element_type_triangle_3, _element_type_triangle_3>, std::tuple<_element_type_cohesive_2d_4, _element_type_quadrangle_4, _element_type_quadrangle_4>, std::tuple<_element_type_cohesive_2d_4, _element_type_triangle_3, _element_type_quadrangle_4>, std::tuple<_element_type_cohesive_2d_6, _element_type_triangle_6, _element_type_triangle_6>, std::tuple<_element_type_cohesive_2d_6, _element_type_quadrangle_8, _element_type_quadrangle_8>, std::tuple<_element_type_cohesive_2d_6, _element_type_triangle_6, _element_type_quadrangle_8>, std::tuple<_element_type_cohesive_3d_6, _element_type_tetrahedron_4, _element_type_tetrahedron_4>, std::tuple<_element_type_cohesive_3d_12, _element_type_tetrahedron_10, _element_type_tetrahedron_10> /*, std::tuple<_element_type_cohesive_3d_8, _element_type_hexahedron_8, _element_type_hexahedron_8>, std::tuple<_element_type_cohesive_3d_16, _element_type_hexahedron_20, _element_type_hexahedron_20>*/>>; TYPED_TEST_SUITE(TestSMMCFixture, coh_types, ); #endif /* AKANTU_TEST_COHESIVE_FIXTURE_HH_ */ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation.cc index f9cd18b61..322a03eaa 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation.cc @@ -1,169 +1,169 @@ /** * @file test_material_standard_linear_solid_deviatoric_relaxation.cc * * @author David Simon Kammer * @author Nicolas Richart * @author Vladislav Yastrebov * * @date creation: Mon Aug 09 2010 * @date last modification: Mon Jun 12 2017 * * @brief test of the viscoelastic material: relaxation * * * Copyright (©) 2010-2018 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 #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" using namespace akantu; int main(int argc, char * argv[]) { akantu::initialize("material_standard_linear_solid_deviatoric_relaxation.dat", argc, argv); akantu::debug::setDebugLevel(akantu::dblWarning); // sim data Real T = 10.; Real eps = 0.001; const UInt dim = 2; Real sim_time = 25.; Real time_factor = 0.1; Real tolerance = 1e-7; Mesh mesh(dim); mesh.read("test_material_standard_linear_solid_deviatoric_relaxation.msh"); - const ElementType element_type = _quadrangle_4; + const auto element_type = _quadrangle_4; SolidMechanicsModel model(mesh); /* ------------------------------------------------------------------------ */ /* Initialization */ /* ------------------------------------------------------------------------ */ model.initFull(); std::cout << model.getMaterial(0) << std::endl; model.assembleMassLumped(); std::stringstream filename_sstr; filename_sstr << "test_material_standard_linear_solid_deviatoric_relaxation.out"; std::ofstream output_data; output_data.open(filename_sstr.str().c_str()); output_data << "#[1]-time [2]-sigma_analytic [3+]-sigma_measurements" << std::endl; Material & mat = model.getMaterial(0); - const Array & stress = mat.getStress(element_type); - Real Eta = mat.get("Eta"); Real EV = mat.get("Ev"); Real Einf = mat.get("Einf"); Real nu = mat.get("nu"); Real Ginf = Einf / (2 * (1 + nu)); Real G = EV / (2 * (1 + nu)); Real G0 = G + Ginf; Real gamma = G / G0; Real tau = Eta / EV; Real gammainf = Ginf / G0; UInt nb_nodes = mesh.getNbNodes(); - const Array & coordinate = mesh.getNodes(); - Array & displacement = model.getDisplacement(); + const auto & coordinate = mesh.getNodes(); + auto & displacement = model.getDisplacement(); /// Setting time step Real time_step = model.getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model.setTimeStep(time_step); UInt max_steps = sim_time / time_step; UInt out_interval = 1; Real time = 0.; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for (UInt s = 0; s <= max_steps; ++s) { - if (s % 1000 == 0) + if (s % 1000 == 0) { std::cerr << "passing step " << s << "/" << max_steps << std::endl; + } time = s * time_step; // impose displacement Real epsilon = 0.; if (time < T) { epsilon = eps * time / T; } else { epsilon = eps; } + for (UInt n = 0; n < nb_nodes; ++n) { displacement(n, 0) = epsilon * coordinate(n, 1); displacement(n, 1) = epsilon * coordinate(n, 0); } + ++model.getDisplacementRelease(); // compute stress model.assembleInternalForces(); // print output if (s % out_interval == 0) { // analytical solution Real solution = 0.; if (time < T) { solution = 2 * G0 * eps / T * (gammainf * time + gamma * tau * (1 - exp(-time / tau))); } else { solution = 2 * G0 * eps * (gammainf + gamma * tau / T * (exp((T - time) / tau) - exp(-time / tau))); } output_data << s * time_step << " " << solution; // data output - Array::const_matrix_iterator stress_it = stress.begin(dim, dim); - Array::const_matrix_iterator stress_end = stress.end(dim, dim); - for (; stress_it != stress_end; ++stress_it) { - output_data << " " << (*stress_it)(0, 1) << " " << (*stress_it)(1, 0); + for (auto && stress : make_view(mat.getStress(element_type), dim, dim)) { + output_data << " " << stress(0, 1) << " " << stress(1, 0); // test error - Real rel_error_1 = std::abs(((*stress_it)(0, 1) - solution) / solution); - Real rel_error_2 = std::abs(((*stress_it)(1, 0) - solution) / solution); + Real rel_error_1 = std::abs((stress(0, 1) - solution) / solution); + Real rel_error_2 = std::abs((stress(1, 0) - solution) / solution); + if (rel_error_1 > tolerance || rel_error_2 > tolerance) { std::cerr << "Relative error: " << rel_error_1 << " " << rel_error_2 << std::endl; return EXIT_FAILURE; } } output_data << std::endl; } } finalize(); std::cout << "Test successful!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation_tension.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation_tension.cc index 685f130d9..1cd8c4a85 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation_tension.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic/test_material_standard_linear_solid_deviatoric_relaxation_tension.cc @@ -1,179 +1,181 @@ /** * @file test_material_standard_linear_solid_deviatoric_relaxation_tension.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Mon Aug 09 2010 * @date last modification: Mon Jun 12 2017 * * @brief test of the viscoelastic material: relaxation * * * Copyright (©) 2010-2018 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 #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" using namespace akantu; int main(int argc, char * argv[]) { akantu::initialize("material_standard_linear_solid_deviatoric_relaxation.dat", argc, argv); // sim data Real T = 10.; Real eps = 0.001; // const UInt dim = 3; const UInt dim = 2; Real sim_time = 25.; // Real sim_time = 250.; Real time_factor = 0.1; Real tolerance = 1e-5; Mesh mesh(dim); mesh.read("test_material_standard_linear_solid_deviatoric_relaxation.msh"); // mesh_io.read("hexa_structured.msh",mesh); // const ElementType element_type = _hexahedron_8; const ElementType element_type = _quadrangle_4; SolidMechanicsModel model(mesh); /* ------------------------------------------------------------------------ */ /* Initialization */ /* ------------------------------------------------------------------------ */ model.initFull(); std::cout << model.getMaterial(0) << std::endl; model.assembleMassLumped(); model.assembleInternalForces(); model.getMaterial(0).setToSteadyState(); std::stringstream filename_sstr; filename_sstr << "test_material_standard_linear_solid_deviatoric_relaxation_" "tension.out"; std::ofstream output_data; output_data.open(filename_sstr.str().c_str()); output_data << "#[1]-time [2]-sigma_analytic [3+]-sigma_measurements" << std::endl; Material & mat = model.getMaterial(0); const Array & stress = mat.getStress(element_type); Real Eta = mat.get("Eta"); Real EV = mat.get("Ev"); Real Einf = mat.get("Einf"); Real E0 = mat.get("E"); Real kpa = mat.get("kapa"); Real mu = mat.get("mu"); Real gamma = EV / E0; Real gammainf = Einf / E0; Real tau = Eta / EV; std::cout << "relaxation time = " << tau << std::endl; UInt nb_nodes = mesh.getNbNodes(); const Array & coordinate = mesh.getNodes(); Array & displacement = model.getDisplacement(); /// Setting time step Real time_step = model.getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model.setTimeStep(time_step); UInt max_steps = sim_time / time_step; UInt out_interval = 1; Real time = 0.; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for (UInt s = 0; s <= max_steps; ++s) { - - if (s % 1000 == 0) + if (s % 1000 == 0) { std::cerr << "passing step " << s << "/" << max_steps << std::endl; + } time = s * time_step; // impose displacement Real epsilon = 0.; if (time < T) { epsilon = eps * time / T; } else { epsilon = eps; } for (UInt n = 0; n < nb_nodes; ++n) { - for (UInt d = 0; d < dim; ++d) + for (UInt d = 0; d < dim; ++d) { displacement(n, d) = epsilon * coordinate(n, d); + } } + ++model.getDisplacementRelease(); // compute stress model.assembleInternalForces(); // print output if (s % out_interval == 0) { // analytical solution Real epskk = dim * eps; Real solution = 0.; if (time < T) { solution = 2 * mu * (eps - epskk / 3.) / T * (gammainf * time + gamma * tau * (1 - exp(-time / tau))) + gammainf * kpa * epskk * time / T; } else { solution = 2 * mu * (eps - epskk / 3.) * (gammainf + gamma * tau / T * (exp((T - time) / tau) - exp(-time / tau))) + gammainf * kpa * epskk; } output_data << s * time_step << " " << solution; // data output Array::const_matrix_iterator stress_it = stress.begin(dim, dim); Array::const_matrix_iterator stress_end = stress.end(dim, dim); for (; stress_it != stress_end; ++stress_it) { output_data << " " << (*stress_it)(1, 1); // test error Real rel_error_1 = std::abs(((*stress_it)(1, 1) - solution) / solution); if (rel_error_1 > tolerance) { std::cerr << "Relative error: " << rel_error_1 << std::endl; return EXIT_FAILURE; } } output_data << std::endl; } } finalize(); std::cout << "Test successful!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_material_eigenstrain.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_material_eigenstrain.cc index 2399256bb..f56b5ef76 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_material_eigenstrain.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_material_eigenstrain.cc @@ -1,200 +1,200 @@ /** * @file test_solid_mechanics_model_material_eigenstrain.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Sat Apr 16 2011 * @date last modification: Thu Feb 01 2018 * * @brief test the internal field prestrain * * * Copyright (©) 2010-2018 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 "mesh_utils.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; Real alpha[3][4] = {{0.01, 0.02, 0.03, 0.04}, {0.05, 0.06, 0.07, 0.08}, {0.09, 0.10, 0.11, 0.12}}; /* -------------------------------------------------------------------------- */ template static Matrix prescribed_strain() { UInt spatial_dimension = ElementClass::getSpatialDimension(); Matrix strain(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = alpha[i][j + 1]; } } return strain; } template static Matrix prescribed_stress(Matrix prescribed_eigengradu) { UInt spatial_dimension = ElementClass::getSpatialDimension(); Matrix stress(spatial_dimension, spatial_dimension); // plane strain in 2d Matrix strain(spatial_dimension, spatial_dimension); Matrix pstrain; pstrain = prescribed_strain(); Real nu = 0.3; Real E = 2.1e11; Real trace = 0; /// symetric part of the strain tensor - for (UInt i = 0; i < spatial_dimension; ++i) - for (UInt j = 0; j < spatial_dimension; ++j) + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); - + } + } // elastic strain is equal to elastic strain minus the eigenstrain strain -= prescribed_eigengradu; for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i, i); Real lambda = nu * E / ((1 + nu) * (1 - 2 * nu)); Real mu = E / (2 * (1 + nu)); if (spatial_dimension == 1) { stress(0, 0) = E * strain(0, 0); } else { - for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt i = 0; i < spatial_dimension; ++i){ for (UInt j = 0; j < spatial_dimension; ++j) { stress(i, j) = (i == j) * lambda * trace + 2 * mu * strain(i, j); } + } } return stress; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_elastic_plane_strain.dat", argc, argv); UInt dim = 3; const ElementType element_type = _tetrahedron_4; Matrix prescribed_eigengradu(dim, dim); prescribed_eigengradu.set(0.1); /// load mesh Mesh mesh(dim); mesh.read("cube_3d_tet_4.msh"); /// declaration of model SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _static); // model.getNewSolver("static", TimeStepSolverType::_static, // NonLinearSolverType::_newton_raphson_modified); auto & solver = model.getNonLinearSolver("static"); solver.set("threshold", 2e-4); solver.set("max_iterations", 2); solver.set("convergence_type", SolveConvergenceCriteria::_residual); const Array & coordinates = mesh.getNodes(); Array & displacement = model.getDisplacement(); Array & boundary = model.getBlockedDOFs(); MeshUtils::buildFacets(mesh); mesh.createBoundaryGroupFromGeometry(); // Loop over (Sub)Boundar(ies) for (auto & group : mesh.iterateElementGroups()) { for (const auto & n : group.getNodeGroup()) { std::cout << "Node " << n << std::endl; for (UInt i = 0; i < dim; ++i) { displacement(n, i) = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); } boundary(n, i) = true; } } } + ++model.getDisplacementRelease(); /* ------------------------------------------------------------------------ */ /* Apply eigenstrain in each element */ /* ------------------------------------------------------------------------ */ Array & eigengradu_vect = model.getMaterial(0).getInternal("eigen_grad_u")(element_type); - auto eigengradu_it = eigengradu_vect.begin(dim, dim); - auto eigengradu_end = eigengradu_vect.end(dim, dim); - - for (; eigengradu_it != eigengradu_end; ++eigengradu_it) { - *eigengradu_it = prescribed_eigengradu; + for (auto & eigen_grad_u : make_view(eigengradu_vect, dim, dim)) { + eigen_grad_u = prescribed_eigengradu; } /* ------------------------------------------------------------------------ */ /* Static solve */ /* ------------------------------------------------------------------------ */ model.solveStep(); std::cout << "Converged in " << Int(solver.get("nb_iterations")) << " (" << Real(solver.get("error")) << ")" << std::endl; /* ------------------------------------------------------------------------ */ /* Checks */ /* ------------------------------------------------------------------------ */ const Array & stress_vect = model.getMaterial(0).getStress(element_type); auto stress_it = stress_vect.begin(dim, dim); auto stress_end = stress_vect.end(dim, dim); Matrix presc_stress; presc_stress = prescribed_stress(prescribed_eigengradu); Real stress_tolerance = 1e-13; for (; stress_it != stress_end; ++stress_it) { const auto & stress = *stress_it; Matrix diff(dim, dim); diff = stress; diff -= presc_stress; Real stress_error = diff.norm() / stress.norm(); if (stress_error > stress_tolerance) { std::cerr << "stress error: " << stress_error << " > " << stress_tolerance << std::endl; std::cerr << "stress: " << stress << std::endl << "prescribed stress: " << presc_stress << std::endl; return EXIT_FAILURE; } // else { // std::cerr << "stress error: " << stress_error // << " < " << stress_tolerance << std::endl; // } } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc index e31ab086a..6eeac6758 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc @@ -1,203 +1,205 @@ /** * @file test_solid_mechanics_model_reassign_material.cc * * @author Aurelia Isabel Cuba Ramos * * @date creation: Mon Feb 10 2014 * @date last modification: Tue Feb 20 2018 * * @brief test the function reassign material * * * Copyright (©) 2014-2018 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_grid_dynamic.hh" #include "communicator.hh" #include "material.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; class StraightInterfaceMaterialSelector : public MaterialSelector { public: StraightInterfaceMaterialSelector(SolidMechanicsModel & model, const std::string & mat_1_material, const std::string & mat_2_material, bool & horizontal, Real & pos_interface) : model(model), mat_1_material(mat_1_material), mat_2_material(mat_2_material), horizontal(horizontal), pos_interface(pos_interface) { Mesh & mesh = model.getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(); /// store barycenters of all elements barycenters.initialize(mesh, _spatial_dimension = spatial_dimension, _nb_component = spatial_dimension); for (auto ghost_type : ghost_types) { Element e; e.ghost_type = ghost_type; for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { UInt nb_element = mesh.getNbElement(type, ghost_type); e.type = type; Array & barycenter = barycenters(type, ghost_type); barycenter.resize(nb_element); Array::iterator> bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { e.element = elem; mesh.getBarycenter(e, *bary_it); ++bary_it; } } } } UInt operator()(const Element & elem) { UInt spatial_dimension = model.getSpatialDimension(); const Vector & bary = barycenters(elem.type, elem.ghost_type) .begin(spatial_dimension)[elem.element]; /// check for a given element on which side of the material interface plane /// the bary center lies and assign corresponding material if (bary(!horizontal) < pos_interface) { return model.getMaterialIndex(mat_1_material); ; } return model.getMaterialIndex(mat_2_material); ; } bool isConditonVerified() { /// check if material has been (re)-assigned correctly Mesh & mesh = model.getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(); GhostType ghost_type = _not_ghost; for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { Array & mat_indexes = model.getMaterialByElement(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); Array::iterator> bary = barycenters(type, ghost_type).begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem, ++bary) { /// compare element_index_by material to material index that should be /// assigned due to the geometry of the interface UInt mat_index; if ((*bary)(!horizontal) < pos_interface) mat_index = model.getMaterialIndex(mat_1_material); else mat_index = model.getMaterialIndex(mat_2_material); if (mat_indexes(elem) != mat_index) /// wrong material index, make test fail return false; } } return true; } void moveInterface(Real & pos_new, bool & horizontal_new) { /// update position and orientation of material interface plane pos_interface = pos_new; horizontal = horizontal_new; model.reassignMaterial(); } protected: SolidMechanicsModel & model; ElementTypeMapArray barycenters; std::string mat_1_material; std::string mat_2_material; bool horizontal; Real pos_interface; }; /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { bool test_passed; debug::setDebugLevel(dblWarning); initialize("two_materials.dat", argc, argv); /// specify position and orientation of material interface plane bool horizontal = true; Real pos_interface = 0.; UInt spatial_dimension = 3; const auto & comm = Communicator::getWorldCommunicator(); Int prank = comm.whoAmI(); Mesh mesh(spatial_dimension); - if (prank == 0) + if (prank == 0) { mesh.read("cube_two_materials.msh"); + } mesh.distribute(); /// model creation SolidMechanicsModel model(mesh); /// assign the two different materials using the /// StraightInterfaceMaterialSelector auto && mat_selector = std::make_shared( model, "mat_1", "mat_2", horizontal, pos_interface); model.setMaterialSelector(mat_selector); model.initFull(_analysis_method = _static); MeshUtils::buildFacets(mesh); /// check if different materials have been assigned correctly test_passed = mat_selector->isConditonVerified(); if (!test_passed) { AKANTU_ERROR("materials not correctly assigned"); return EXIT_FAILURE; } /// change orientation of material interface plane horizontal = false; mat_selector->moveInterface(pos_interface, horizontal); // model.dump(); /// test if material has been reassigned correctly test_passed = mat_selector->isConditonVerified(); if (!test_passed) { AKANTU_ERROR("materials not correctly reassigned"); return EXIT_FAILURE; } finalize(); - if (prank == 0) + if (prank == 0) { std::cout << "OK: test passed!" << std::endl; + } return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */