diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc index 0f840ef5d..e97fdf316 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,604 +1,604 @@ /** * @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 "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), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// find the four corner nodes of the RVE findCornerNodes(); /// 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()); // /// enforce periodicity on the displacement fluctuations // auto surface_pair_1 = std::make_pair("top", "bottom"); // auto surface_pair_2 = std::make_pair("right", "left"); // SurfacePairList surface_pairs_list; // surface_pairs_list.push_back(surface_pair_1); // surface_pairs_list.push_back(surface_pair_2); // TODO: To Nicolas correct the PBCs // this->setPBC(surface_pairs_list); 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); // this->initMaterials(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); /// compute the volume of the RVE GhostType gt = _not_ghost; for (auto element_type : this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { Array Volume(this->mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem.integrate(Volume, element_type); } 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->dump(); this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will /// be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyHomogeneousTemperature( const Real & temperature) { for (UInt m = 0; m < this->getNbMaterials(); ++m) { Material & mat = this->getMaterial(m); const ElementTypeMapArray & filter_map = mat.getElementFilter(); - Mesh::type_iterator type_it = filter_map.firstType(spatial_dimension); - Mesh::type_iterator type_end = filter_map.lastType(spatial_dimension); // Loop over all element types - for (; type_it != type_end; ++type_it) { - const Array & filter = filter_map(*type_it); + for (auto && type : filter_map.elementTypes(spatial_dimension)) { + const Array & filter = filter_map(type); if (filter.size() == 0) continue; - Array & delta_T = mat.getArray("delta_T", *type_it); + Array & delta_T = mat.getArray("delta_T", type); Array::scalar_iterator delta_T_it = delta_T.begin(); Array::scalar_iterator delta_T_end = delta_T.end(); for (; delta_T_it != delta_T_end; ++delta_T_it) { *delta_T_it = temperature; } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { AKANTU_DEBUG_IN(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, spatial_dimension))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } for (UInt i = 0; i < corner_nodes.size(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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 MaterialDamageIterative<2> & mat_paste = dynamic_cast &>(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast &>(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); do { 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()); std::cout << "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(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); Real average = 0; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & stress_vec = this->materials[m]->getStress(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec( k, row_index * spatial_dimension + col_index); // 3 is the value // for the yy (in // 3D, the value is // 4) } } else if (field_type == "strain") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & gradu_vec = this->materials[m]->getGradU(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else { AKANTU_ERROR("Averaging not implemented for this field!!!"); } } return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// clear the eigenstrain Matrix zero_eigengradu(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & 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) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = // 0); // this->fillCracks(saved_damage); /// virtual test 1: H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); this->solveStep(); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU( Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u"); 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) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top - bottom); Real eps = box_size * 1e-6; auto tmp = std::make_shared(*this, box_size, "gel", this->nb_gel_pockets, eps); 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(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = this->getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; Real E = mat->get("E"); auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks( const ElementTypeMapReal & saved_damage) { for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.cc b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.cc index 76c1e14b8..58612f76f 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.cc +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.cc @@ -1,255 +1,255 @@ /** * @file material_damage_iterative.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative.hh" #include "communicator.hh" #include "data_accessor.hh" #include "solid_mechanics_model_RVE.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialDamageIterative::MaterialDamageIterative( SolidMechanicsModel & model, const ID & id) : MaterialDamage(model, id), Sc("Sc", *this), reduction_step("damage_step", *this), equivalent_stress("equivalent_stress", *this), max_reductions(0), norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold"); this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "prescribed damage"); this->registerParam( "dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1"); this->registerParam( "dam_tolerance", dam_tolerance, 0.01, _pat_parsable | _pat_modifiable, "damage tolerance to decide if quadrature point will be damageed"); this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value"); this->registerParam("max_reductions", max_reductions, UInt(10), _pat_parsable | _pat_modifiable, "max reductions"); this->use_previous_stress = true; this->use_previous_gradu = true; this->Sc.initialize(1); this->equivalent_stress.initialize(1); this->reduction_step.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative:: computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); auto Sc_it = Sc(el_type, ghost_type).begin(); auto equivalent_stress_it = equivalent_stress(el_type, ghost_type).begin(); Array::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); Real * dam = this->damage(el_type, ghost_type).storage(); Matrix sigma(spatial_dimension, spatial_dimension); for (; grad_u_it != grad_u_end; ++grad_u_it) { sigma.clear(); MaterialElastic::computeStressOnQuad(*grad_u_it, sigma, 0.); computeDamageAndStressOnQuad(sigma, *dam); /// compute eigenvalues sigma.eig(eigenvalues); /// find max eigenvalue and normalize by tensile strength *equivalent_stress_it = *(std::max_element(eigenvalues.storage(), eigenvalues.storage() + spatial_dimension)) / *(Sc_it); ++Sc_it; ++equivalent_stress_it; ++dam; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::computeAllStresses( GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress if (ghost_type == _not_ghost) norm_max_equivalent_stress = 0; MaterialDamage::computeAllStresses(ghost_type); /// find global Gauss point with highest stress auto rve_model = dynamic_cast(&this->model); if (rve_model == NULL) { /// is no RVE model const auto & comm = this->model.getMesh().getCommunicator(); comm.allReduce(norm_max_equivalent_stress, SynchronizerOperation::_max); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative:: findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (ghost_type == _not_ghost) { // const Array & e_stress = equivalent_stress(el_type); // if (e_stress.begin() != e_stress.end() ) { // auto equivalent_stress_it_max = // std::max_element(e_stress.begin(),e_stress.end()); // /// check if max equivalent stress for this element type is greater // than the current norm_max_eq_stress // if (*equivalent_stress_it_max > norm_max_equivalent_stress) // norm_max_equivalent_stress = *equivalent_stress_it_max; // } const Array & e_stress = equivalent_stress(el_type); auto equivalent_stress_it = e_stress.begin(); auto equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); auto dam_it = dam.begin(); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it) { /// check if max equivalent stress for this element type is greater than /// the current norm_max_eq_stress and if the element is not already fully /// damaged if (*equivalent_stress_it > norm_max_equivalent_stress && *dam_it < max_damage) { norm_max_equivalent_stress = *equivalent_stress_it; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialDamage::computeStress(el_type, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeDamageAndStressOnQuad(sigma, *dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type); norm_max_equivalent_stress = 0; findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialDamageIterative::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., "Your prescribed damage must be greater than zero"); if (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; ; Mesh::type_iterator it = this->model.getFEEngine().getMesh().firstType( spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model.getFEEngine().getMesh().lastType(spatial_dimension, ghost_type); for (; it != last_type; ++it) { ElementType el_type = *it; const Array & e_stress = equivalent_stress(el_type); auto equivalent_stress_it = e_stress.begin(); auto equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); auto dam_it = dam.begin(); auto reduction_it = this->reduction_step(el_type, ghost_type).begin(); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it, ++reduction_it) { /// check if damage occurs if (*equivalent_stress_it >= (1 - dam_tolerance) * norm_max_equivalent_stress) { /// check if this element can still be damaged if (*reduction_it == this->max_reductions) continue; *reduction_it += 1; if (*reduction_it == this->max_reductions) { *dam_it = max_damage; } else { *dam_it += prescribed_dam; } 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; } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::updateEnergiesAfterDamage( - ElementType el_type, GhostType ghost_type) { - MaterialDamage::updateEnergies(el_type, ghost_type); + ElementType el_type) { + MaterialDamage::updateEnergies(el_type); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(damage_iterative, MaterialDamageIterative); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh index acd902544..b5ef98296 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh @@ -1,147 +1,146 @@ /** * @file material_damage_iterative.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ #define __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ namespace akantu { /** * Material damage iterative * * parameters in the material files : * - Sc */ template class MaterialDamageIterative : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamageIterative(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialDamageIterative(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// virtual void updateInternalParameters(); virtual void computeAllStresses(GhostType ghost_type = _not_ghost); /// update internal field damage virtual UInt updateDamage(); UInt updateDamage(UInt quad_index, const Real eq_stress, const ElementType & el_type, const GhostType & ghost_type); /// update energies after damage has been updated - virtual void updateEnergiesAfterDamage(ElementType el_type, - GhostType ghost_typ); + virtual void updateEnergiesAfterDamage(ElementType el_type); /// compute the equivalent stress on each Gauss point (i.e. the max prinicpal /// stress) and normalize it by the tensile strength virtual void computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type = _not_ghost); /// find max normalized equivalent stress void findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ 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; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get max normalized equivalent stress AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real); /// get a non-const reference to the max normalized equivalent stress AKANTU_GET_MACRO_NOT_CONST(NormMaxEquivalentStress, norm_max_equivalent_stress, Real &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// resistance to damage RandomInternalField Sc; /// the reduction InternalField reduction_step; /// internal field to store equivalent stress on each Gauss point InternalField equivalent_stress; /// the number of total reductions steps until complete failure UInt max_reductions; /// damage increment Real prescribed_dam; /// maximum equivalent stress Real norm_max_equivalent_stress; /// deviation from max stress at which Gauss point will still get damaged Real dam_tolerance; /// define damage threshold at which damage will be set to 1 Real dam_threshold; /// maximum damage value Real max_damage; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative_inline_impl.cc" #endif /* __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh index afd4555ce..69261f49a 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh @@ -1,147 +1,147 @@ /** * @file material_orthotropic_damage.hh * @author Aurelia Cuba Ramos * @date Sun Mar 8 12:49:56 2015 * * @brief Material for orthotropic damage * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH__ #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH__ namespace akantu { template class Parent = MaterialElastic> class MaterialOrthotropicDamage : public Parent { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialOrthotropicDamage(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialOrthotropicDamage(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// compute the tangent stiffness matrix for an element type virtual void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); protected: /// update the dissipated energy, must be called after the stress have been /// computed - virtual void updateEnergies(ElementType el_type, GhostType ghost_type); + void updateEnergies(ElementType el_type) override; /// compute the tangent stiffness matrix for a given quadrature point inline void computeTangentModuliOnQuad( Matrix & tangent, const Matrix C, const Matrix & dam, const Matrix & dam_directions, Matrix & O_prime, Matrix & S_prime, Matrix & O, Matrix & S, Matrix & rotation_tmp); inline void computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & root_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term); /// rotate a Matrix of size dim*dim into the coordinate system of the FE /// computation inline void rotateIntoComputationFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp); /// rotate a Matrix of size dim*dim into the coordinate system of the damage inline void rotateIntoNewFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp); /// compute (1-D) inline void computeOneMinusD(Matrix & one_minus_D, const Matrix & damage); /// compute (1-D)^(1/2) inline void computeSqrtOneMinusD(const Matrix & one_minus_D, Matrix & sqrt_one_minus_D); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; // virtual Real getEnergy(std::string type); // virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) // { // return Parent::getEnergy(energy_id, type, index); // }; AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ElementTypeMapArray &); AKANTU_GET_MACRO(Damage, damage, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real) /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// damage internal variable InternalField damage; /// dissipated energy InternalField dissipated_energy; /// contain the current value of @f$ \int_0^{\epsilon}\sigma(\omega)d\omega /// @f$ the dissipated energy InternalField int_sigma; /// direction vectors for the damage frame InternalField damage_dir_vecs; Real eta; /// maximum damage value Real max_damage; }; } // namespace akantu #include "material_orthotropic_damage_tmpl.hh" #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.cc b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.cc index 73d3cc282..2f67f7d40 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.cc +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.cc @@ -1,378 +1,376 @@ /** * @file material_damage_iterative.cc * * @author Aurelia Isabel Cuba Ramos * * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage_iterative.hh" #include "communicator.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialOrthotropicDamageIterative:: MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, const ID & id) : MaterialOrthotropicDamage(model, id), Sc("Sc", *this), equivalent_stress("equivalent_stress", *this), stress_dir("equiv_stress_dir", *this), norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold"); this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "increase of damage in every step"); this->registerParam( "dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1"); this->use_previous_stress = true; this->use_previous_gradu = true; this->Sc.initialize(1); this->equivalent_stress.initialize(1); this->stress_dir.initialize(spatial_dimension * spatial_dimension); /// the Gauss point with the highest stress can only be of type _not_ghost q_max.ghost_type = _not_ghost; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative:: computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); auto Sc_it = Sc(el_type).begin(); auto equivalent_stress_it = equivalent_stress(el_type).begin(); Array::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_dir_it = this->stress_dir(el_type).begin(spatial_dimension, spatial_dimension); /// initialize matrix to store the stress for the criterion (only different in /// non-local) Matrix sigma(spatial_dimension, spatial_dimension); Array::matrix_iterator damage_iterator = this->damage(el_type, ghost_type) .begin(this->spatial_dimension, this->spatial_dimension); Array::matrix_iterator damage_dir_it = this->damage_dir_vecs(el_type, ghost_type) .begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); Matrix rotation_tmp(this->spatial_dimension, this->spatial_dimension); /// create matrix to store the first term of the computation of the /// Cauchy stress Matrix first_term(this->spatial_dimension, this->spatial_dimension); Matrix third_term(this->spatial_dimension, this->spatial_dimension); for (; grad_u_it != grad_u_end; ++Sc_it, ++equivalent_stress_it, ++stress_dir_it, ++grad_u_it) { sigma.clear(); MaterialOrthotropicDamage::computeStressOnQuad( *grad_u_it, sigma, 0.); /// rotate the tensors from the damage principal coordinate system to the CS /// of the computation if (!(Math::are_float_equal((*damage_iterator).trace(), 0))) { /// compute (1-D) and (1-D)^1/2 this->computeOneMinusD(one_minus_D, *damage_iterator); this->computeSqrtOneMinusD(one_minus_D, sqrt_one_minus_D); this->rotateIntoComputationFrame(one_minus_D, one_minus_D_rotated, *damage_dir_it, rotation_tmp); this->rotateIntoComputationFrame(sqrt_one_minus_D, sqrt_one_minus_D_rotated, *damage_dir_it, rotation_tmp); } else { this->computeOneMinusD(one_minus_D_rotated, *damage_iterator); this->computeSqrtOneMinusD(one_minus_D_rotated, sqrt_one_minus_D_rotated); } computeDamageAndStressOnQuad(sigma, one_minus_D_rotated, sqrt_one_minus_D_rotated, *damage_iterator, first_term, third_term); /// compute the maximum principal stresses and their directions sigma.eig(eigenvalues, *stress_dir_it); *equivalent_stress_it = eigenvalues(0) / *(Sc_it); ++damage_dir_it; ++damage_iterator; } // for(;sigma_it != sigma_end; ++sigma_it, // ++Sc_it, ++equivalent_stress_it, ++stress_dir_it) { // /// compute the maximum principal stresses and their directions // (*sigma_it).eig(eigenvalues, *stress_dir_it); // *equivalent_stress_it = eigenvalues(0) / *(Sc_it); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative::computeAllStresses( GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress if (ghost_type == _not_ghost) norm_max_equivalent_stress = 0; MaterialOrthotropicDamage::computeAllStresses(ghost_type); /// find global Gauss point with highest stress this->model.getMesh().getCommunicator().allReduce( norm_max_equivalent_stress, SynchronizerOperation::_max); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative:: findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (ghost_type == _not_ghost) { /// initialize the iterators for the equivalent stress and the damage const Array & e_stress = equivalent_stress(el_type); auto equivalent_stress_it = e_stress.begin(); auto equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); auto damage_directions_it = this->damage_dir_vecs(el_type, _not_ghost) .begin(this->spatial_dimension, this->spatial_dimension); auto stress_dir_it = this->stress_dir(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); /// initialize the matrices for damage rotation results Matrix tmp(spatial_dimension, spatial_dimension); Matrix dam_in_computation_frame(spatial_dimension, spatial_dimension); Matrix dam_in_stress_frame(spatial_dimension, spatial_dimension); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it, ++damage_directions_it, ++stress_dir_it) { /// check if max equivalent stress for this element type is greater than /// the current norm_max_eq_stress if (*equivalent_stress_it > norm_max_equivalent_stress && (spatial_dimension * this->max_damage - (*dam_it).trace() > Math::getTolerance())) { if (Math::are_float_equal((*dam_it).trace(), 0)) { /// gauss point has not been damaged yet norm_max_equivalent_stress = *equivalent_stress_it; q_max.type = el_type; q_max.global_num = equivalent_stress_it - e_stress.begin(); } else { /// find the damage increment on this Gauss point /// rotate damage into stress frame this->rotateIntoComputationFrame(*dam_it, dam_in_computation_frame, *damage_directions_it, tmp); this->rotateIntoNewFrame(dam_in_computation_frame, dam_in_stress_frame, *stress_dir_it, tmp); /// add damage increment dam_in_stress_frame(0, 0) += prescribed_dam; /// find new principal directions of damage Vector dam_eigenvalues(spatial_dimension); dam_in_stress_frame.eig(dam_eigenvalues); bool limit_reached = false; for (UInt i = 0; i < spatial_dimension; ++i) { if (dam_eigenvalues(i) + Math::getTolerance() > this->max_damage) limit_reached = true; } if (!limit_reached) { norm_max_equivalent_stress = *equivalent_stress_it; q_max.type = el_type; q_max.global_num = equivalent_stress_it - e_stress.begin(); } } } /// end if equiv_stress > max_equiv_stress } /// end loop over all gauss points of this element type } // end if(_not_ghost) AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialOrthotropicDamage::computeStress(el_type, ghost_type); auto damage_iterator = this->damage(el_type, ghost_type) .begin(this->spatial_dimension, this->spatial_dimension); auto damage_dir_it = this->damage_dir_vecs(el_type, ghost_type) .begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); Matrix rotation_tmp(this->spatial_dimension, this->spatial_dimension); /// create matrix to store the first term of the computation of the /// Cauchy stress Matrix first_term(this->spatial_dimension, this->spatial_dimension); Matrix third_term(this->spatial_dimension, this->spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); /// rotate the tensors from the damage principal coordinate system to the CS /// of the computation if (!(Math::are_float_equal((*damage_iterator).trace(), 0))) { /// compute (1-D) and (1-D)^1/2 this->computeOneMinusD(one_minus_D, *damage_iterator); this->computeSqrtOneMinusD(one_minus_D, sqrt_one_minus_D); this->rotateIntoComputationFrame(one_minus_D, one_minus_D_rotated, *damage_dir_it, rotation_tmp); this->rotateIntoComputationFrame(sqrt_one_minus_D, sqrt_one_minus_D_rotated, *damage_dir_it, rotation_tmp); } else { this->computeOneMinusD(one_minus_D_rotated, *damage_iterator); this->computeSqrtOneMinusD(one_minus_D_rotated, sqrt_one_minus_D_rotated); } computeDamageAndStressOnQuad(sigma, one_minus_D_rotated, sqrt_one_minus_D_rotated, *damage_iterator, first_term, third_term); ++damage_dir_it; ++damage_iterator; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type); norm_max_equivalent_stress = 0; findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialOrthotropicDamageIterative::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., "Your prescribed damage must be greater than zero"); if (norm_max_equivalent_stress >= 1.) { AKANTU_DEBUG_IN(); /// get the arrays and iterators for the element_type of the highest /// quadrature point ElementType el_type = q_max.type; UInt q_global_num = q_max.global_num; Array & dam = this->damage(el_type, _not_ghost); auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); auto damage_directions_it = this->damage_dir_vecs(el_type, _not_ghost) .begin(this->spatial_dimension, this->spatial_dimension); auto stress_dir_it = this->stress_dir(el_type, _not_ghost) .begin(spatial_dimension, spatial_dimension); /// initialize the matrices for damage rotation results Matrix tmp(spatial_dimension, spatial_dimension); Matrix dam_in_computation_frame(spatial_dimension, spatial_dimension); Matrix dam_in_stress_frame(spatial_dimension, spatial_dimension); /// references to damage and directions of highest Gauss point Matrix q_dam = dam_it[q_global_num]; Matrix q_dam_dir = damage_directions_it[q_global_num]; Matrix q_stress_dir = stress_dir_it[q_global_num]; /// increment damage /// find the damage increment on this Gauss point /// rotate damage into stress frame this->rotateIntoComputationFrame(q_dam, dam_in_computation_frame, q_dam_dir, tmp); this->rotateIntoNewFrame(dam_in_computation_frame, dam_in_stress_frame, q_stress_dir, tmp); /// add damage increment dam_in_stress_frame(0, 0) += prescribed_dam; /// find new principal directions of damage Vector dam_eigenvalues(spatial_dimension); dam_in_stress_frame.eig(dam_eigenvalues, q_dam_dir); for (UInt i = 0; i < spatial_dimension; ++i) { q_dam(i, i) = dam_eigenvalues(i); if (q_dam(i, i) + Math::getTolerance() >= dam_threshold) q_dam(i, i) = this->max_damage; } nb_damaged_elements += 1; } this->model.getMesh().getCommunicator().allReduce( nb_damaged_elements, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative< - spatial_dimension>::updateEnergiesAfterDamage(ElementType el_type, - GhostType ghost_type) { - MaterialOrthotropicDamage::updateEnergies(el_type, - ghost_type); + spatial_dimension>::updateEnergiesAfterDamage(ElementType el_type) { + MaterialOrthotropicDamage::updateEnergies(el_type); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(orthotropic_damage_iterative, MaterialOrthotropicDamageIterative); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh index 91a98795a..7ab01bafe 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh @@ -1,322 +1,322 @@ /** * @file material_orthotropic_damage_tmpl.hh * @author Aurelia Isabel Cuba Ramos * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the material class for the orthotropic * damage material * * @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 "material_orthotropic_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> MaterialOrthotropicDamage::MaterialOrthotropicDamage( SolidMechanicsModel & model, const ID & id) : Parent(model, id), damage("damage", *this), dissipated_energy("damage dissipated energy", *this), int_sigma("integral of sigma", *this), damage_dir_vecs("damage_principal_directions", *this) { AKANTU_DEBUG_IN(); this->registerParam("eta", eta, 2., _pat_parsable | _pat_modifiable, "Damage sensitivity parameter"); this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value"); this->is_non_local = false; this->use_previous_stress = true; this->use_previous_gradu = true; /// use second order tensor for description of damage state this->damage.initialize(spatial_dimension * spatial_dimension); this->dissipated_energy.initialize(1); this->int_sigma.initialize(1); this->damage_dir_vecs.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the dissipated energy in each element by a trapezoidal approximation * of * @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - * \frac{1}{2}\sigma:\epsilon@f$ */ template class Parent> void MaterialOrthotropicDamage::updateEnergies( - ElementType el_type, GhostType ghost_type) { - Parent::updateEnergies(el_type, ghost_type); + ElementType el_type) { + Parent::updateEnergies(el_type); - this->computePotentialEnergy(el_type, ghost_type); + this->computePotentialEnergy(el_type); - auto epsilon_p = this->gradu.previous(el_type, ghost_type) + auto epsilon_p = this->gradu.previous(el_type) .begin(spatial_dimension, spatial_dimension); - auto sigma_p = this->stress.previous(el_type, ghost_type) + auto sigma_p = this->stress.previous(el_type) .begin(spatial_dimension, spatial_dimension); - auto epot = this->potential_energy(el_type, ghost_type).begin(); + auto epot = this->potential_energy(el_type).begin(); - auto ints = this->int_sigma(el_type, ghost_type).begin(); - auto ed = this->dissipated_energy(el_type, ghost_type).begin(); + auto ints = this->int_sigma(el_type).begin(); + auto ed = this->dissipated_energy(el_type).begin(); - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); - Matrix delta_gradu_it(*gradu_it); + Matrix delta_gradu_it(grad_u); delta_gradu_it -= *epsilon_p; Matrix sigma_h(sigma); sigma_h += *sigma_p; Real dint = .5 * sigma_h.doubleDot(delta_gradu_it); *ints += dint; *ed = *ints - *epot; ++epsilon_p; ++sigma_p; ++epot; ++ints; ++ed; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeTangentModuli(el_type, tangent_matrix, ghost_type); /// get the damage array for current element type Array & dam = this->damage(el_type); auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); /// get the directions of damage for the current element type Array & dam_dirs = this->damage_dir_vecs(el_type); auto damage_directions_it = dam_dirs.begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix one_minus_D_rot(spatial_dimension, spatial_dimension); Matrix sqrt_one_minus_D_rot(spatial_dimension, spatial_dimension); Matrix rotation_tmp(spatial_dimension, spatial_dimension); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); if (!(Math::are_float_equal((*dam_it).trace(), 0))) computeTangentModuliOnQuad(tangent, tangent, *dam_it, *damage_directions_it, one_minus_D, sqrt_one_minus_D, one_minus_D_rot, sqrt_one_minus_D_rot, rotation_tmp); ++dam_it; ++damage_directions_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage:: computeTangentModuliOnQuad(Matrix & tangent, const Matrix C, const Matrix & dam, const Matrix & dam_directions, Matrix & O_prime, Matrix & S_prime, Matrix & O, Matrix & S, Matrix & rotation_tmp) { /// effect of damage on stiffness matrix: See Ragueneau et al. 2008, p. 423, /// ep. 7 Real trace_D = dam.trace(); this->computeOneMinusD(O_prime, dam); this->computeSqrtOneMinusD(O_prime, S_prime); this->rotateIntoComputationFrame(O_prime, O, dam_directions, rotation_tmp); this->rotateIntoComputationFrame(S_prime, S, dam_directions, rotation_tmp); /// compute new stiffness matrix in damage coordinate system if (spatial_dimension == 1) tangent *= (1 - dam(0, 0)); if (spatial_dimension == 2) { Real min_val = std::min((this->eta / spatial_dimension * trace_D), this->max_damage); /// first row tangent(0, 0) = (C(0, 0) * S(0, 0) * S(0, 0) + C(1, 0) * S(0, 1) * S(0, 1) - (min_val / 2. - 1. / 2) * (C(0, 0) + C(1, 0)) + (O(0, 0) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.)); tangent(0, 1) = (C(0, 1) * S(0, 0) * S(0, 0) + C(1, 1) * S(0, 1) * S(0, 1) - (min_val / 2. - 1. / 2) * (C(0, 1) + C(1, 1)) + (O(0, 0) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.)); tangent(0, 2) = (2. * C(2, 2) * S(0, 0) * S(0, 1) + (2. * C(2, 2) * O(0, 0) * O(0, 1)) / (trace_D - 2.)); /// second row tangent(1, 0) = (C(0, 0) * S(0, 1) * S(0, 1) + C(1, 0) * S(1, 1) * S(1, 1) - (min_val / 2. - 1. / 2) * (C(0, 0) + C(1, 0)) + (O(1, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.)); tangent(1, 1) = (C(0, 1) * S(0, 1) * S(0, 1) + C(1, 1) * S(1, 1) * S(1, 1) - (min_val / 2. - 1. / 2) * (C(0, 1) + C(1, 1)) + (O(1, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.)); tangent(1, 2) = (2. * C(2, 2) * S(0, 1) * S(1, 1) + (2. * C(2, 2) * O(0, 1) * O(1, 1)) / (trace_D - 2.)); /// third row tangent(2, 0) = ((O(0, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.) + C(0, 0) * S(0, 0) * S(0, 1) + C(1, 0) * S(0, 1) * S(1, 1)); tangent(2, 1) = ((O(0, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.) + C(0, 1) * S(0, 0) * S(0, 1) + C(1, 1) * S(0, 1) * S(1, 1)); tangent(2, 2) = ((2. * C(2, 2) * O(0, 1) * O(0, 1)) / (trace_D - 2.) + C(2, 2) * S(0, 1) * S(0, 1) + C(2, 2) * S(0, 0) * S(1, 1)); } if (spatial_dimension == 3) { } } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage:: computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & sqrt_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term) { /// Definition of Cauchy stress based on second order damage tensor: /// "Anisotropic damage modelling of biaxial behaviour and rupture /// of concrete strucutres", Ragueneau et al., 2008, Eq. 7 first_term.mul(sqrt_one_minus_D, sigma); first_term *= sqrt_one_minus_D; Real second_term = 0; for (UInt i = 0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) second_term += sigma(i, j) * one_minus_D(i, j); } second_term /= (this->spatial_dimension - damage.trace()); one_minus_D *= second_term; third_term.eye(1. / this->spatial_dimension * sigma.trace() * (1 - eta / (this->spatial_dimension) * damage.trace())); sigma.copy(first_term); sigma -= one_minus_D; sigma += third_term; } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage:: rotateIntoComputationFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::rotateIntoNewFrame( const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::computeOneMinusD( Matrix & one_minus_D, const Matrix & damage) { /// compute one minus one_minus_D.eye(); one_minus_D -= damage; } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::computeSqrtOneMinusD( const Matrix & one_minus_D, Matrix & sqrt_one_minus_D) { /// To compute (1-D)^1/2 we need to check that we are in the /// principal coordinate system of the damage #ifndef AKANTU_NDEBUG for (UInt i = 0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) { if (i != j) AKANTU_DEBUG_ASSERT(Math::are_float_equal(one_minus_D(i, j), 0), "The damage tensor has off-diagonal parts"); } } #endif // AKANTU_NDEBUG /// compute (1-D)^1/2 sqrt_one_minus_D.copy(one_minus_D); for (UInt i = 0; i < this->spatial_dimension; ++i) sqrt_one_minus_D(i, i) = std::sqrt(sqrt_one_minus_D(i, i)); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc index 48b476fb9..8d1b0a2bc 100644 --- a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc +++ b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc @@ -1,108 +1,108 @@ /** * @file material_viscoplastic.cc * * @author Ramin Aghababaei * * * @brief Specialization of the material class for isotropic viscoplastic * (small deformation) * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_viscoplastic.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialViscoPlastic::MaterialViscoPlastic(SolidMechanicsModel & model, const ID & id) : MaterialPlastic(model, id) { AKANTU_DEBUG_IN(); this->registerParam("rate", rate, 0., _pat_parsable | _pat_modifiable, "Rate sensitivity component"); this->registerParam("edot0", edot0, 0., _pat_parsable | _pat_modifiable, "Reference strain rate"); this->registerParam("ts", ts, 0., _pat_parsable | _pat_modifiable, "Time Step"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoPlastic::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * iso_hardening = this->iso_hardening(el_type, ghost_type).storage(); auto previous_grad_u_it = this->gradu.previous(el_type, ghost_type).begin(dim, dim); auto previous_sigma_it = this->stress.previous(el_type, ghost_type).begin(dim, dim); auto inelastic_strain_it = this->inelastic_strain(el_type, ghost_type).begin(dim, dim); auto previous_inelastic_strain_it = this->inelastic_strain.previous(el_type, ghost_type).begin(dim, dim); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeStressOnQuad(grad_u, *previous_grad_u_it, sigma, *previous_sigma_it, *inelastic_strain_it, *previous_inelastic_strain_it, *iso_hardening); ++inelastic_strain_it; ++iso_hardening; ++previous_grad_u_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoPlastic::computeTangentModuli( __attribute__((unused)) const ElementType & el_type, Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); auto previous_sigma_it = this->stress.previous(el_type, ghost_type).begin(dim, dim); auto previous_strain_it = this->gradu.previous(el_type, ghost_type).begin(dim, dim); Real * iso_hardening = this->iso_hardening(el_type, ghost_type).storage(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); Matrix & previous_grad_u = *previous_strain_it; Matrix & previous_sigma_tensor = *previous_sigma_it; - computeTangentModuliOnQuad(tangent, grad_u, previous_grad_u, sigma_tensor, + computeTangentModuliOnQuad(tangent, grad_u, previous_grad_u, sigma, previous_sigma_tensor, *iso_hardening); ++previous_sigma_it; ++previous_strain_it; ++iso_hardening; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } INSTANTIATE_MATERIAL(visco_plastic, MaterialViscoPlastic); } // namespace akantu