diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc index 3aa14eb6a..7aa4435a7 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc @@ -1,214 +1,214 @@ /** * @file material_damage_iterative.cc * @author Aurelia Cuba Ramos * @date Fri Feb 14 14:54:44 2014 * * @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-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_damage_iterative.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialDamageIterative::MaterialDamageIterative(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), Sc("Sc", *this), equivalent_stress("equivalent_stress", *this), norm_max_equivalent_stress(0), q_point() { 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_strain = true; + this->use_previous_stress = true; + this->use_previous_gradu = true; this->Sc.initialize(1); this->equivalent_stress.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); Array::const_iterator Sc_it = Sc(el_type).begin(); Array::iterator 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); 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 Gauss point to store the information where highest critical stress occurs if(ghost_type==_not_ghost) { norm_max_equivalent_stress = 0; q_point.type = _not_defined; q_point.global_num = 0; q_point.ghost_type = ghost_type; } MaterialDamage::computeAllStresses(ghost_type); /// find global Gauss point with highest stress StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); Real global_norm_max_equivalent_stress = norm_max_equivalent_stress; comm.allReduce(&global_norm_max_equivalent_stress, 1, _so_max); if (!Math::are_float_equal(global_norm_max_equivalent_stress, norm_max_equivalent_stress)) { q_point.type = _not_defined; q_point.global_num = 0; } 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() ) { Array::const_iterator 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; q_point.type = el_type; q_point.global_num = equivalent_stress_it_max-e_stress.begin(); } } } 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->strain(el_type, ghost_type), el_type, ghost_type); + 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 (q_point.type != _not_defined) { //const Array & Sc_vect = material.getInternal("Sc")(q_point.type); const Array & e_stress_vect = equivalent_stress(q_point.type); /// check if damage occurs if (e_stress_vect(q_point.global_num) >= 1.) { //Sc_vect(q_point.global_num)) { /// update internal field damage and increment # if damaged elements Array & dam_vect = this->damage(q_point.type); if (dam_vect(q_point.global_num) < dam_threshold) dam_vect(q_point.global_num) += prescribed_dam; else dam_vect(q_point.global_num) = 1.; nb_damaged_elements += 1; } } StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { MaterialDamage::updateEnergies(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialDamageIterative); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_plastic/material_viscoplastic.cc b/src/model/solid_mechanics/materials/material_plastic/material_viscoplastic.cc index a802921fc..af3920d4e 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_viscoplastic.cc +++ b/src/model/solid_mechanics/materials/material_plastic/material_viscoplastic.cc @@ -1,118 +1,118 @@ /** * @file material_viscoplastic.cc * * @author Ramin Aghababaei * * @date Tue Jul 09 18:15:37 20130 * * @brief Specialization of the material class for isotropic viscoplastic (small deformation) * * @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_viscoplastic.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialViscoPlastic::MaterialViscoPlastic(SolidMechanicsModel & model, const ID & id) : Material(model, 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(); Array::matrix_iterator previous_grad_u_it = - this->strain.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); + this->gradu.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator previous_sigma_it = this->stress.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator inelastic_strain_it = this->inelastic_strain(el_type, ghost_type).begin(spatial_dimension,spatial_dimension); Array::matrix_iterator previous_inelastic_strain_it = this->inelastic_strain.previous(el_type, ghost_type).begin(spatial_dimension,spatial_dimension); 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(); Array::matrix_iterator previous_sigma_it = this->stress.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator previous_strain_it = - this->strain.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); + this->gradu.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); 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, previous_sigma_tensor, *iso_hardening); ++previous_sigma_it; ++previous_strain_it; ++iso_hardening; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } INSTANSIATE_MATERIAL(MaterialViscoPlastic); __END_AKANTU__ diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit_anisotropic.cc b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit_anisotropic.cc index b4def1758..bc7e5329b 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit_anisotropic.cc +++ b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit_anisotropic.cc @@ -1,302 +1,302 @@ /** * @file patch_test_explicit.cc * * @author David Simon Kammer * @author Nicolas Richart * @author Cyprien Wolff * * @date Thu Feb 17 16:05:48 2011 * * @brief patch test for elastic material in solid mechanics model * * @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 #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 } }; //Stiffness tensor, rotated by hand Real C[3][3][3][3] = {{{{112.93753505, 1.85842452538e-10, -4.47654358027e-10}, {1.85847317471e-10, 54.2334345331, -3.69840984824}, {-4.4764768395e-10, -3.69840984824, 56.848605217}}, {{1.85847781609e-10, 25.429294233, -3.69840984816}, {25.429294233, 3.31613847493e-10, -8.38797920011e-11}, {-3.69840984816, -8.38804581349e-11, -1.97875715813e-10}}, {{-4.47654358027e-10, -3.69840984816, 28.044464917}, {-3.69840984816, 2.09374961813e-10, 9.4857455224e-12}, {28.044464917, 9.48308098714e-12, -2.1367885239e-10}}}, {{{1.85847781609e-10, 25.429294233, -3.69840984816}, {25.429294233, 3.31613847493e-10, -8.38793479119e-11}, {-3.69840984816, -8.38795699565e-11, -1.97876381947e-10}}, {{54.2334345331, 3.31617400207e-10, 2.09372075233e-10}, {3.3161562385e-10, 115.552705733, -3.15093728886e-10}, {2.09372075233e-10, -3.15090176173e-10, 54.2334345333}}, {{-3.69840984824, -8.38795699565e-11, 9.48219280872e-12}, {-8.38795699565e-11, -3.1509195253e-10, 25.4292942335}, {9.48441325477e-12, 25.4292942335, 3.69840984851}}}, {{{-4.47653469848e-10, -3.69840984816, 28.044464917}, {-3.69840984816, 2.09374073634e-10, 9.48752187924e-12}, {28.044464917, 9.48552347779e-12, -2.1367885239e-10}}, {{-3.69840984824, -8.3884899027e-11, 9.48219280872e-12}, {-8.3884899027e-11, -3.150972816e-10, 25.4292942335}, {9.48041645188e-12, 25.4292942335, 3.69840984851}}, {{56.848605217, -1.97875493768e-10, -2.13681516925e-10}, {-1.97877270125e-10, 54.2334345333, 3.69840984851}, {-2.13683293282e-10, 3.69840984851, 112.93753505}}}}; /* -------------------------------------------------------------------------- */ template static Matrix prescribed_grad_u() { UInt spatial_dimension = ElementClass::getSpatialDimension(); Matrix grad_u(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { grad_u(i, j) = alpha[i][j + 1]; } } return grad_u; } template static Matrix prescribed_stress() { 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_grad_u(); 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) strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); for (UInt i = 0 ; i < spatial_dimension ; ++i) { for (UInt j = 0 ; j < spatial_dimension ; ++j) { stress(i, j) = 0; for (UInt k = 0 ; k < spatial_dimension ; ++k) { for (UInt l = 0 ; l < spatial_dimension ; ++l) { stress(i, j) += C[i][j][k][l]*strain(k, l); } } } } return stress; } #define TYPE _tetrahedron_4 /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize("material_anisotropic.dat", argc, argv); UInt dim = ElementClass::getSpatialDimension(); const ElementType element_type = TYPE; UInt damping_steps = 600000; UInt damping_interval = 50; Real damping_ratio = 0.99; UInt additional_steps = 20000; UInt max_steps = damping_steps + additional_steps; /// load mesh Mesh my_mesh(dim); std::stringstream filename; filename << TYPE << ".msh"; my_mesh.read(filename.str()); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass)); std::cout << my_model.getMaterial(0) << std::endl; Real time_step = my_model.getStableTimeStep()/5.; my_model.setTimeStep(time_step); my_model.assembleMassLumped(); std::cout << "The number of time steps is: " << max_steps << " (" << time_step << "s)" << std::endl; // boundary conditions const Array & coordinates = my_mesh.getNodes(); Array & displacement = my_model.getDisplacement(); Array & boundary = my_model.getBlockedDOFs(); MeshUtils::buildFacets(my_mesh); my_mesh.createBoundaryGroupFromGeometry(); // Loop over (Sub)Boundar(ies) // Loop over (Sub)Boundar(ies) for(GroupManager::const_element_group_iterator it(my_mesh.element_group_begin()); it != my_mesh.element_group_end(); ++it) { for(ElementGroup::const_node_iterator nodes_it(it->second->node_begin()); nodes_it!= it->second->node_end(); ++nodes_it) { UInt n(*nodes_it); std::cout << "Node " << *nodes_it << 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; } } } // Actually, loop over all nodes, since I wanna test a static solution for (UInt n = 0 ; n < nb_nodes ; ++n) { 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); } } } Array & velocity = my_model.getVelocity(); std::ofstream energy; std::stringstream energy_filename; energy_filename << "energy_" << TYPE << ".csv"; energy.open(energy_filename.str().c_str()); energy << "id,time,ekin" << std::endl; Real ekin_mean = 0.; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ UInt s; for(s = 1; s <= max_steps; ++s) { if(s % 10000 == 0) std::cout << "passing step " << s << "/" << max_steps << " (" << s*time_step << "s)" < & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); - Array & strain_vect = const_cast &>(my_model.getMaterial(0).getStrain(element_type)); + Array & gradu_vect = const_cast &>(my_model.getMaterial(0).getGradU(element_type)); Array::iterator< Matrix > stress_it = stress_vect.begin(dim, dim); - Array::iterator< Matrix > strain_it = strain_vect.begin(dim, dim); + Array::iterator< Matrix > gradu_it = gradu_vect.begin(dim, dim); Matrix presc_stress; presc_stress = prescribed_stress(); - Matrix presc_strain; presc_strain = prescribed_grad_u(); + Matrix presc_gradu; presc_gradu = prescribed_grad_u(); UInt nb_element = my_mesh.getNbElement(TYPE); - Real strain_tolerance = 1e-9; + Real gradu_tolerance = 1e-9; Real stress_tolerance = 1e-2; if(s > max_steps) { stress_tolerance = 1e-4; - strain_tolerance = 1e-7; + gradu_tolerance = 1e-7; } for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & stress = *stress_it; - Matrix & strain = *strain_it; + Matrix & gradu = *gradu_it; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { - if(!(std::abs(strain(i, j) - presc_strain(i, j)) < strain_tolerance)) { - std::cerr << "strain[" << i << "," << j << "] = " << strain(i, j) << " but should be = " << presc_strain(i, j) << " (-" << std::abs(strain(i, j) - presc_strain(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; - std::cerr << strain << presc_strain << std::endl; + if(!(std::abs(gradu(i, j) - presc_gradu(i, j)) < gradu_tolerance)) { + std::cerr << "gradu[" << i << "," << j << "] = " << gradu(i, j) << " but should be = " << presc_gradu(i, j) << " (-" << std::abs(gradu(i, j) - presc_gradu(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; + std::cerr << gradu << presc_gradu << std::endl; return EXIT_FAILURE; } if(!(std::abs(stress(i, j) - presc_stress(i, j)) < stress_tolerance)) { std::cerr << "stress[" << i << "," << j << "] = " << stress(i, j) << " but should be = " << presc_stress(i, j) << " (-" << std::abs(stress(i, j) - presc_stress(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << stress << presc_stress << std::endl; return EXIT_FAILURE; } } } ++stress_it; - ++strain_it; + ++gradu_it; } } for (UInt n = 0; n < nb_nodes; ++n) { for (UInt i = 0; i < dim; ++i) { Real disp = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { disp += alpha[i][j + 1] * coordinates(n, j); } if(!(std::abs(displacement(n,i) - disp) < 1e-7)) { std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << "(" << displacement(n,i) - disp << ")" << std::endl; return EXIT_FAILURE; } } } finalize(); return EXIT_SUCCESS; }