diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc index 84faf9648..65a302b74 100644 --- a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc +++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc @@ -1,199 +1,199 @@ /** * @file material_FE2.cc * * @author Aurelia Isabel Cuba Ramos * * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * @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_FE2.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialFE2::MaterialFE2(SolidMechanicsModel & model, const ID & id) : Material(model, id), Parent(model, id), C("material_stiffness", *this) { AKANTU_DEBUG_IN(); this->C.initialize(voigt_h::size * voigt_h::size); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialFE2::~MaterialFE2() { AKANTU_DEBUG_IN(); for (UInt i = 0; i < RVEs.size(); ++i) { delete meshes[i]; delete RVEs[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initialize() { this->registerParam("element_type" ,el_type, _triangle_3 , _pat_parsable | _pat_modifiable, "element type in RVE mesh" ); this->registerParam("mesh_file" ,mesh_file , _pat_parsable | _pat_modifiable, "the mesh file for the RVE"); this->registerParam("nb_gel_pockets" ,nb_gel_pockets , _pat_parsable | _pat_modifiable, "the number of gel pockets in each RVE"); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// compute the number of integration points in this material and resize the RVE vector UInt nb_integration_points = this->element_filter(this->el_type, _not_ghost).getSize() * this->fem->getNbIntegrationPoints(this->el_type); RVEs.resize(nb_integration_points); meshes.resize(nb_integration_points); /// create a Mesh and SolidMechanicsModel on each integration point of the material std::vector::iterator RVE_it = RVEs.begin(); std::vector::iterator mesh_it = meshes.begin(); StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); UInt prank = comm.whoAmI(); Array::matrix_iterator C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size); for (UInt i = 1; i < nb_integration_points+1; ++RVE_it, ++mesh_it, ++i, ++C_it) { std::stringstream mesh_name; mesh_name << "RVE_mesh_" << prank; std::stringstream rve_name; rve_name << "SMM_RVE_" << prank; *mesh_it = new Mesh(spatial_dimension, mesh_name.str(), i); (*mesh_it)->read(mesh_file); *RVE_it = new SolidMechanicsModelRVE(*(*(mesh_it)), true, this->nb_gel_pockets, _all_dimensions, rve_name.str(), i); (*RVE_it)->initFull(); /// compute intial stiffness of the RVE (*RVE_it)->homogenizeStiffness(*C_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Compute thermal stresses first Parent::computeStress(el_type, ghost_type); Array::const_scalar_iterator sigma_th_it = this->sigma_th(el_type, ghost_type).begin(); // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); // create vectors to store stress and strain in Voigt notation // for efficient computation of stress Vector voigt_strain(voigt_h::size); Vector voigt_stress(voigt_h::size); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); const Matrix & C_mat = *C_it; const Real & sigma_th = *sigma_th_it; /// copy strains in Voigt notation for(UInt I = 0; I < voigt_h::size; ++I) { /// copy stress in Real voigt_factor = voigt_h::factors[I]; UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; } // compute stresses in Voigt notation voigt_stress.mul(C_mat, voigt_strain); /// copy stresses back in full vectorised notation for(UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; sigma(i, j) = sigma(j, i) = voigt_stress(I)+ (i == j) * sigma_th; } ++C_it; ++sigma_th_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); tangent.copy(*C_it); ++C_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); std::vector::iterator RVE_it = RVEs.begin(); std::vector::iterator RVE_end = RVEs.end(); Array::matrix_iterator C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size); Array::matrix_iterator gradu_it = this->gradu(this->el_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator eigen_gradu_it = this->eigengradu(this->el_type).begin(spatial_dimension, spatial_dimension); for (; RVE_it != RVE_end; ++RVE_it, ++C_it, ++gradu_it, ++eigen_gradu_it) { /// apply boundary conditions based on the current macroscopic displ. gradient (*RVE_it)->applyBoundaryConditions(*gradu_it); /// advance the ASR in every RVE (*RVE_it)->advanceASR(prestrain); /// compute the average eigen_grad_u (*RVE_it)->homogenizeEigenGradU(*eigen_gradu_it); /// compute the new effective stiffness of the RVE (*RVE_it)->homogenizeStiffness(*C_it); } AKANTU_DEBUG_OUT(); } INSTANTIATE_MATERIAL(MaterialFE2); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh index 02f0843d0..0ddd75ba7 100644 --- a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh +++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh @@ -1,123 +1,123 @@ /** * @file material_FE2.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * @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_thermal.hh" #include "solid_mechanics_model_RVE.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_FE_2_HH__ #define __AKANTU_MATERIAL_FE_2_HH__ -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ /// /!\ This material works ONLY for meshes with a single element type!!!!! /* -------------------------------------------------------------------------- */ /** * MaterialFE2 * * parameters in the material files : * - mesh_file */ // Emil - 27.01.2018 - re-inheriting from MaterialThermal and PlaneStressToolBox //template //class MaterialFE2 : public virtual Material { // /* ------------------------------------------------------------------------ */ // /* Constructors/Destructors */ // /* ------------------------------------------------------------------------ */ template class MaterialFE2 : public MaterialThermal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: typedef MaterialThermal Parent; // Emil - 27.01.2018 public: MaterialFE2(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialFE2(); typedef VoigtHelper voigt_h; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); /// advance alkali-silica reaction void advanceASR(const Matrix & prestrain); private: void initialize(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying RVE at each integration point std::vector RVEs; /// Meshes for all RVEs std::vector meshes; /// the element type of the associated mesh (this material handles only one type!!) ElementType el_type; /// the name of RVE mesh file ID mesh_file; /// Elastic stiffness tensor at each Gauss point (in voigt notation) InternalField C; /// number of gel pockets in each underlying RVE UInt nb_gel_pockets; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_FE2_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_FE_2_HH__ */ 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 f830f6cba..aefff6a96 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,600 +1,600 @@ /** * @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 "material_damage_iterative.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_PETSC #include "solver_petsc.hh" #include "petsc_matrix.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_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), static_communicator_dummy(StaticCommunicator::getStaticCommunicatorDummy()), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// create node groups for PBCs mesh.createGroupsFromMeshData("physical_names"); /// 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 ElementGroup & bottom = mesh.getElementGroup("bottom"); bottom_nodes.insert( bottom.node_begin(), bottom.node_end() ); const ElementGroup & left = mesh.getElementGroup("left"); left_nodes.insert( left.node_begin(), left.node_end() ); /// enforce periodicity on the displacement fluctuations SurfacePair surface_pair_1 = std::make_pair("top","bottom"); SurfacePair 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); this->setPBC(surface_pairs_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() { delete static_communicator_dummy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); SolidMechanicsModel::initFull(options); this->initMaterials(); /// compute the volume of the RVE FEEngine * fem = this->fems["SolidMechanicsFEEngine"]; GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; Array Volume(this->mesh.getNbElement(element_type) * this->fems["SolidMechanicsFEEngine"]->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 ("force"); this->addDumpField ("equivalent_stress"); this->addDumpField ("residual"); 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::findCornerNodes() { AKANTU_DEBUG_IN(); mesh.computeBoundingBox(); // find corner nodes const Array & position = mesh.getNodes(); const Vector & lower_bounds = mesh.getLowerBounds(); const Vector & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { // node 1 if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && Math::are_float_equal(position(n, 1), lower_bounds(1))) { corner_nodes(0) = n; } // node 2 else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) && Math::are_float_equal(position(n, 1), lower_bounds(1))) { corner_nodes(1) = n; } // node 3 else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) && Math::are_float_equal(position(n, 1), upper_bounds(1))) { corner_nodes(2) = n; } // node 4 else if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && Math::are_float_equal(position(n, 1), upper_bounds(1))) { corner_nodes(3) = n; } } for (UInt i = 0; i < corner_nodes.getSize(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_DEBUG_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 GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal("eigen_grad_u")(element_type, gt)); Array::iterator< Matrix > prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > 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; Real error = 0; bool converged = false; do { converged = this->solveStep<_scm_newton_raphson_tangent, _scc_increment>(1e-6, error, 2, false, *static_communicator_dummy); AKANTU_DEBUG_ASSERT(converged, "Did not converge"); std::cout << "the error is " << error << std::endl; /// 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(); FEEngine * fem = this->fems["SolidMechanicsFEEngine"]; Real average = 0; GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; if (field_type == "stress") { for(UInt m = 0; m < this->materials.size(); ++m) { const Array & stress_vec = this->materials[m]->getStress(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.getSize(), 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.getSize(); ++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 Array & gradu_vec = this->materials[m]->getGradU(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.getSize(), 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.getSize(); ++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 Array & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.getSize(), 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.getSize(); ++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_DEBUG_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 GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(dim, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(dim, gt, _ek_not_defined); Matrix zero_eigengradu(dim, dim, 0.); for(; it != end; ++it) { const ElementType element_type = *it; Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal("eigen_grad_u")(element_type, gt)); Array::iterator< Matrix > prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > 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 fillin up cracks ElementTypeMapReal saved_damage("saved_damage"); mesh.initElementTypeMapArray(saved_damage, 1, spatial_dimension, false, _ek_regular, true); saved_damage.clear(); /// fill the cracks for the tension test // 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); /// solve system this->assembleStiffnessMatrix(); Real error = 0; bool converged= this->solveStep<_scm_newton_raphson_tangent_not_computed, _scc_increment>(1e-6, error, 2, false, *static_communicator_dummy); std::cout << "error in tension test " << error << std::endl; AKANTU_DEBUG_ASSERT(converged, "Did not converge"); /// 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; if(is_default_material_selector) delete material_selector; material_selector = new GelMaterialSelector(*this, box_size, "gel", this->nb_gel_pockets, eps); is_default_material_selector = false; } SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initSolver(__attribute__((unused)) SolverOptions & options) { AKANTU_DEBUG_IN(); #if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #else jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _unsymmetric, sstr.str(), memory_id, 1); #endif //AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id, 1); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #endif //AKANTU_USE_PETSC } delete solver; std::stringstream sstr_solv; sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str(), memory_id); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str(), memory_id); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS SolverMumpsOptions opt(SolverMumpsOptions::_serial_split); if(solver) solver->initialize(opt); #endif //AKANTU_HAS_SOLVER AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension, *static_communicator_dummy); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const Material & mat_gel = this->getMaterial("gel"); const Real E_gel = mat_gel.getParam("E"); Real E_homogenized = 0.; GhostType gt = _not_ghost; for (UInt m = 0; m < this->getNbMaterials(); ++m) { Material & mat = this->getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; const Real E = mat.getParam("E"); InternalField & damage = mat.getInternal("damage"); Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); for (; it != end; ++it) { const ElementType element_type = *it; const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.getSize()) continue; Array & damage_vec = damage(element_type, gt); Array & saved_damage_vec = saved_damage(element_type, gt); for (UInt i = 0; i < damage_vec.getSize(); ++i) { saved_damage_vec(elem_filter(i)) = damage_vec(i); E_homogenized = (E_gel - E) * damage_vec(i) + E; damage_vec(i) = 1. -(E_homogenized/E); } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks(const ElementTypeMapReal & saved_damage) { GhostType gt = _not_ghost; for (UInt m = 0; m < this->getNbMaterials(); ++m) { Material & mat = this->getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; else { InternalField & damage = mat.getInternal("damage"); Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); for (; it != end; ++it) { const ElementType element_type = *it; const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.getSize()) continue; Array & damage_vec = damage(element_type, gt); const Array & saved_damage_vec = saved_damage(element_type, gt); for (UInt i = 0; i < damage_vec.getSize(); ++i) { damage_vec(i) = saved_damage_vec(elem_filter(i)); } } } } } -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh index 61a9467bb..d83557c30 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh @@ -1,249 +1,249 @@ /** * @file solid_mechanics_model_RVE.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 14:54:18 2016 * * @brief SMM for RVE computations in FE2 simulations * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "aka_grid_dynamic.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class SolidMechanicsModelRVE : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector = true, UInt nb_gel_pockets = 400, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModelRVE(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(const ModelOptions & options = SolidMechanicsModelOptions(_static, true)); /// initialize the materials virtual void initMaterials(); /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditions(const Matrix & displacement_gradient); /// advance the reactions -> grow gel and apply homogenized properties void advanceASR(const Matrix & prestrain); /// compute average stress or strain in the model Real averageTensorField(UInt row_index, UInt col_index, const ID & field_type); /// compute effective stiffness of the RVE void homogenizeStiffness(Matrix & C_macro); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /// initialize the solver and the jacobian_matrix (called by initImplicit) virtual void initSolver(SolverOptions & options = _solver_no_options); /// allocate all vectors virtual void initArrays(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(Volume, volume, Real); private: /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); void fillCracks(ElementTypeMapReal & saved_damage); void drainCracks(const ElementTypeMapReal & saved_damage); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ /// volume of the RVE Real volume; /// corner nodes 1, 2, 3, 4 (see Leonardo's thesis, page 98) Array corner_nodes; /// bottom nodes std::unordered_set bottom_nodes; /// left nodes std::unordered_set left_nodes; /// standard mat selector or user one bool use_RVE_mat_selector; StaticCommunicator * static_communicator_dummy; /// the number of gel pockets inside the RVE UInt nb_gel_pockets; /// dump counter UInt nb_dumps; }; inline void SolidMechanicsModelRVE::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) { SolidMechanicsModel::unpackData(buffer, index, tag); if (tag == _gst_smm_uv) { Array::vector_iterator disp_it = displacement->begin(spatial_dimension); Vector current_disp(disp_it[index]); // if node is at the bottom, u_bottom = u_top +u_2 -u_3 if ( bottom_nodes.count(index) ) { current_disp += Vector(disp_it[corner_nodes(1)]); current_disp -= Vector(disp_it[corner_nodes(2)]); } // if node is at the left, u_left = u_right +u_4 -u_3 else if ( left_nodes.count(index) ) { current_disp += Vector(disp_it[corner_nodes(3)]); current_disp -= Vector(disp_it[corner_nodes(2)]); } } } /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ class GelMaterialSelector : public MeshDataMaterialSelector { public: GelMaterialSelector(SolidMechanicsModel & model, const Real box_size, const std::string & gel_material, const UInt nb_gel_pockets, Real tolerance = 0.) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), box_size(box_size) { Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = model.getSpatialDimension(); ElementType type = _triangle_3; GhostType ghost_type = _not_ghost; UInt nb_element = mesh.getNbElement(type, ghost_type); Element el; el.type = type; el.ghost_type = ghost_type; Array barycenter(0,spatial_dimension); barycenter.resize(nb_element); Array::vector_iterator bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++bary_it, ++elem) { mesh.getBarycenter(elem, type, bary_it->storage(), ghost_type); } /// generate the gel pockets srand(0.); Vector center(spatial_dimension); UInt placed_gel_pockets = 0; std::set checked_baries; while (placed_gel_pockets != nb_gel_pockets) { /// get a random bary center UInt bary_id = rand() % nb_element; if (checked_baries.find(bary_id) != checked_baries.end()) continue; checked_baries.insert(bary_id); el.element = bary_id; if (MeshDataMaterialSelector::operator()(el) == 1) continue; /// element belongs to paste gel_pockets.push_back(el); placed_gel_pockets += 1; } } UInt operator()(const Element & elem) { UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index == 1) return temp_index; std::vector::const_iterator iit = gel_pockets.begin(); std::vector::const_iterator eit = gel_pockets.end(); if(std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material);; } return 0; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; }; -__END_AKANTU__ +} // namespace akantu ///#include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle.cc b/extra_packages/extra-materials/src/material_damage/material_brittle.cc index 8e3c3d964..f3d03636f 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle.cc +++ b/extra_packages/extra-materials/src/material_damage/material_brittle.cc @@ -1,122 +1,121 @@ /** * @file material_brittle.cc * * @author Aranda Ruiz Josue * @author Daniel Pino Muñoz * * * @brief Specialization of the material class for the brittle material * * @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_brittle.hh" #include "solid_mechanics_model.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialBrittle::MaterialBrittle(SolidMechanicsModel & model, const ID & id): - Material(model, id), MaterialDamage(model, id), strain_rate_brittle("strain_rate_brittle", *this){ AKANTU_DEBUG_IN(); - this->registerParam("S_0", S_0, 157e6, ParamAccessType(_pat_parsable | _pat_modifiable)); + this->registerParam("S_0", S_0, 157e6, _pat_parsable | _pat_modifiable); this->registerParam("E_0", E_0, 27e3, _pat_parsable, "Strain rate threshold"); this->registerParam("A", A, 1.622e-5, _pat_parsable, "Polynome cubic constant"); this->registerParam("B", B, -1.3274, _pat_parsable, "Polynome quadratic constant"); this->registerParam("C", C, 3.6544e4, _pat_parsable, "Polynome linear constant"); this->registerParam("D", D, -181.38e6, _pat_parsable, "Polynome constant"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); this->strain_rate_brittle.initialize(spatial_dimension * spatial_dimension); updateInternalParameters(); //this->Yd.resize(); - // const Mesh & mesh = this->model->getFEEngine().getMesh(); + // const Mesh & mesh = this->model.getFEEngine().getMesh(); // Mesh::type_iterator it = mesh.firstType(spatial_dimension); // Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); // for(; it != last_type; ++it) { // UInt nb_element = this->element_filter(*it).getSize(); - // UInt nb_quad = this->model->getFEEngine().getNbQuadraturePoints(*it); + // UInt nb_quad = this->model.getFEEngine().getNbQuadraturePoints(*it); // Array & Yd_rand_vec = Yd_rand(*it); // for(UInt e = 0; e < nb_element; ++e) { // Real rand_part = (2 * drand48()-1) * Yd_randomness * Yd; // for(UInt q = 0; q < nb_quad; ++q) // Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::updateInternalParameters() { MaterialDamage::updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); - Array & velocity = this->model->getVelocity(); + Array & velocity = this->model.getVelocity(); Array & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type); Array & elem_filter = this->element_filter(el_type, ghost_type); - this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle, + this->model.getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle, spatial_dimension, el_type, ghost_type, elem_filter); Array::iterator< Matrix > strain_rate_it = this->strain_rate_brittle(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Real sigma_equivalent = 0; Real fracture_stress = 0; Matrix & grad_v = *strain_rate_it; computeStressOnQuad(grad_u, grad_v, sigma, *dam, sigma_equivalent, fracture_stress); ++strain_rate_it; ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -INSTANTIATE_MATERIAL(MaterialBrittle); +INSTANTIATE_MATERIAL(brittle, MaterialBrittle); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle.hh b/extra_packages/extra-materials/src/material_damage/material_brittle.hh index 5c64e998d..e86165c63 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle.hh +++ b/extra_packages/extra-materials/src/material_damage/material_brittle.hh @@ -1,125 +1,111 @@ /** * @file material_brittle.hh * * @author Aranda Ruiz Josue * @author Daniel Pino Muñoz * * * @brief Brittle damage law * * @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_damage.hh" #include "material.hh" +#include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_BRITTLE_HH__ #define __AKANTU_MATERIAL_BRITTLE_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material brittle * * parameters in the material files : * - S_0 : Critical stress at low strain rate (default: 157e6) * - E_0 : Low strain rate threshold (default: 27e3) - * - A,B,C,D : Fitting parameters for the critical stress at high strain rates + * - A,B,C,D : Fitting parameters for the critical stress at high strain + * rates * (default: 1.622e-11, -1.3274e-6, 3.6544e-2, -181.38) */ -template +template class MaterialBrittle : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MaterialBrittle(SolidMechanicsModel & model, const ID & id = ""); - virtual ~MaterialBrittle() {}; - /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - void initMaterial(); virtual void updateInternalParameters(); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for a given quadrature point - inline void computeStressOnQuad(Matrix & grad_u, - Matrix & grad_v, - Matrix & sigma, - Real & dam, - Real & sigma_equivalent, + inline void computeStressOnQuad(Matrix & grad_u, Matrix & grad_v, + Matrix & sigma, Real & dam, + Real & sigma_equivalent, Real & fracture_stress); - inline void computeDamageAndStressOnQuad(Matrix & sigma, - Real & dam, - Real & sigma_c, - Real & fracture_stress); + inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam, + Real & sigma_c, + Real & fracture_stress); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; - inline virtual UInt getNbDataForElements(const Array & elements, - SynchronizationTag tag) const; - - inline virtual void packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const; + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; - inline virtual void unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag); - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ -public: + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// strain rate arrays ordered by element types InternalField strain_rate_brittle; - //polynome constants for critical stress value + // polynome constants for critical stress value Real A; Real B; Real C; Real D; - //minimum strain rate + // minimum strain rate Real E_0; - //Critical stress at low strain rates + // Critical stress at low strain rates Real S_0; - }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_brittle_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_brittle_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_brittle_inline_impl.cc index dc9dd21b0..d88ce1ac9 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_brittle_inline_impl.cc @@ -1,139 +1,101 @@ /** * @file material_brittle_inline_impl.cc * * @author Aranda Ruiz Josue * @author Daniel Pino Muñoz * * * @brief Implementation of the inline functions of the material brittle * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ -template -inline void - MaterialBrittle::computeStressOnQuad(Matrix & grad_u, - Matrix & grad_v, - Matrix & sigma, - Real & dam, - Real & sigma_equivalent, - Real & fracture_stress) { +template +inline void MaterialBrittle::computeStressOnQuad( + Matrix & grad_u, Matrix & grad_v, Matrix & sigma, + Real & dam, Real & sigma_equivalent, Real & fracture_stress) { MaterialElastic::computeStressOnQuad(grad_u, sigma); - Real equiv_strain_rate=0.; - Real volume_change_rate=grad_v.trace(); - if(spatial_dimension==2){ -// if(this->plane_stress){ -// Real e_dot_33 = this->lambda*(volume_change_rate)/(2. * this->mu - this->lambda); -// volume_change_rate+=e_dot_33; -// equiv_strain_rate+=2./3.*pow( e_dot_33 - volume_change_rate/3. , 2. ); -// } -// else - equiv_strain_rate+=2./3.*pow( volume_change_rate/3. , 2. ); + Real equiv_strain_rate = 0.; + Real volume_change_rate = grad_v.trace(); + if (spatial_dimension == 2) { + equiv_strain_rate += 2. / 3. * pow(volume_change_rate / 3., 2.); } - for(UInt i=0; iE_0) - //fracture_stress=A*pow(equiv_strain_rate,3)+B*pow(equiv_strain_rate,2)+C*equiv_strain_rate+D; - fracture_stress=A; + fracture_stress = S_0; + if (equiv_strain_rate > E_0) + fracture_stress = A; Vector principal_stress(spatial_dimension); sigma.eig(principal_stress); sigma_equivalent = principal_stress(0); - for(UInt i=1; iis_non_local) { + if (!this->is_non_local) { computeDamageAndStressOnQuad(sigma, dam, sigma_equivalent, fracture_stress); } } /* -------------------------------------------------------------------------- */ -template -inline void -MaterialBrittle::computeDamageAndStressOnQuad(Matrix & sigma, - Real & dam, - Real & sigma_c, - Real & fracture_stress) { - if ( sigma_c > fracture_stress ) +template +inline void MaterialBrittle::computeDamageAndStressOnQuad( + Matrix & sigma, Real & dam, Real & sigma_c, Real & fracture_stress) { + if (sigma_c > fracture_stress) dam = 1.; - dam = std::min(dam,1.); + dam = std::min(dam, 1.); - sigma *= 1-dam; + sigma *= 1 - dam; } /* -------------------------------------------------------------------------- */ -template -inline UInt MaterialBrittle::getNbDataForElements(const Array & elements, - SynchronizationTag tag) const { +template +inline UInt MaterialBrittle::getNbData( + const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); - /*UInt size = 0; - if(tag == _gst_smm_init_mat) { - size += sizeof(Real) * this->getModel().getNbQuadraturePoints(elements); - } - */ - UInt size = MaterialDamage::getNbDataForElements(elements, tag); + UInt size = MaterialDamage::getNbData(elements, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ -template -inline void MaterialBrittle::packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const { +template +inline void MaterialBrittle::packData( + CommunicationBuffer & buffer, const Array & elements, + const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); - /*if(tag == _gst_smm_init_mat) { - this->packElementDataHelper(Yd, buffer, elements); - }*/ - - MaterialDamage::packElementData(buffer, elements, tag); + MaterialDamage::packData(buffer, elements, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -inline void MaterialBrittle::unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) { +template +inline void +MaterialBrittle::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); - /*if(tag == _gst_smm_init_mat) { - this->unpackElementDataHelper(Yd, buffer, elements); - }*/ - - MaterialDamage::unpackElementData(buffer, elements, tag); + MaterialDamage::unpackData(buffer, elements, tag); AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local.hh b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local.hh index b8e6dd731..b77c270fe 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local.hh +++ b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local.hh @@ -1,85 +1,85 @@ /** * @file material_brittle_non_local.hh * * @author Daniel Pino Muñoz * * * @brief MaterialBrittleNonLocal header for non-local damage * * @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_brittle.hh" #include "material_damage_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_BRITTLE_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_BRITTLE_NON_LOCAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material Brittle Non local * * parameters in the material files : */ template class MaterialBrittleNonLocal : public MaterialDamageNonLocal > { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialDamageNonLocal > MaterialBrittleNonLocalParent; MaterialBrittleNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialBrittleNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); protected: /// constitutive law void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost); /// associate the non-local variables of the material to their neighborhoods virtual void nonLocalVariableToNeighborhood(); private: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Sigma_max, Sigma_max, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: InternalField Sigma_max; InternalField Sigma_maxnl; InternalField Sigma_fracture; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_brittle_non_local_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_BRITTLE_NON_LOCAL_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.cc index 216bc78d6..8c5261e5d 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.cc @@ -1,111 +1,111 @@ /** * @file material_brittle_non_local_inline_impl.cc * * @author Daniel Pino Muñoz * * * @brief MaterialBrittleNonLocal inline function implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ -__END_AKANTU__ +} // namespace akantu #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialBrittleNonLocal::MaterialBrittleNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialBrittleNonLocalParent(model, id), Sigma_max("Sigma max", *this), Sigma_maxnl("Sigma_max non local", *this), Sigma_fracture("Sigma_fracture", *this){ AKANTU_DEBUG_IN(); this->is_non_local = true; this->Sigma_max.initialize(1); this->Sigma_maxnl.initialize(1); this->Sigma_fracture.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::initMaterial() { AKANTU_DEBUG_IN(); - this->model->getNonLocalManager().registerNonLocalVariable(this->Sigma_max.getName(), Sigma_maxnl.getName(),1); + this->model.getNonLocalManager().registerNonLocalVariable(this->Sigma_max.getName(), Sigma_maxnl.getName(),1); MaterialBrittleNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Sigma_maxt = this->Sigma_max(el_type, ghost_type).storage(); Real * fracture_stress = this->Sigma_fracture(el_type, ghost_type).storage(); - Array & velocity = this->model->getVelocity(); + Array & velocity = this->model.getVelocity(); Array & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type); Array & elem_filter = this->element_filter(el_type, ghost_type); - this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle, + this->model.getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle, spatial_dimension, el_type, ghost_type, elem_filter); Array::iterator< Matrix > strain_rate_it = this->strain_rate_brittle(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & grad_v = *strain_rate_it; MaterialBrittle::computeStressOnQuad(grad_u, grad_v, sigma, *dam, *Sigma_maxt, *fracture_stress); ++dam; ++strain_rate_it; ++Sigma_maxt; ++fracture_stress; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::computeNonLocalStress(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(type, ghost_type).storage(); Real * Sigma_maxnlt = this->Sigma_maxnl(type, ghost_type).storage(); Real * fracture_stress = this->Sigma_fracture(type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(type, ghost_type); this->computeDamageAndStressOnQuad(sigma, *dam, *Sigma_maxnlt, *fracture_stress); ++dam; ++Sigma_maxnlt; ++fracture_stress; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::nonLocalVariableToNeighborhood() { - this->model->getNonLocalManager().nonLocalVariableToNeighborhood(Sigma_maxnl.getName(), this->name); + this->model.getNonLocalManager().nonLocalVariableToNeighborhood(Sigma_maxnl.getName(), this->name); } 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 039847eca..e7387097b 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,234 +1,234 @@ /** * @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 "solid_mechanics_model_RVE.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialDamageIterative::MaterialDamageIterative(SolidMechanicsModel & model, const ID & id) : Material(model, 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); Array::const_iterator Sc_it = Sc(el_type, ghost_type).begin(); Array::iterator 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 const SolidMechanicsModelRVE * rve_model = dynamic_cast(this->model); if (rve_model == NULL) { /// is no RVE model StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&norm_max_equivalent_stress, 1, _so_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() ) { // 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; // } const Array & e_stress = equivalent_stress(el_type); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator 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); + 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); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator dam_it = dam.begin(); Array::scalar_iterator 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; } } } } const SolidMechanicsModelRVE * rve_model = dynamic_cast(this->model); if (rve_model == NULL) { StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); } AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { MaterialDamage::updateEnergies(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialDamageIterative); -__END_AKANTU__ +} // 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 8287af9ea..88d620503 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,147 @@ /** * @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_damage.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ #define __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ -__BEGIN_AKANTU__ +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 onBeginningSolveStep(const AnalysisMethod & method) { }; virtual void onEndSolveStep(const AnalysisMethod & method) { }; ///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 getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* 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; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.cc b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.cc index 063462e9f..5c0102369 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.cc +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.cc @@ -1,40 +1,40 @@ /** * @file material_damage_iterative.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Implementation of MaterialDamageIterativeNonLocal * * @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_non_local.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeNonLocal::computeNonLocalStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress if(ghost_type==_not_ghost) this->norm_max_equivalent_stress = 0; MaterialDamageIterativeNonLocalParent::computeNonLocalStresses(ghost_type); /// find global Gauss point with highest stress StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&(this->norm_max_equivalent_stress), 1, _so_max); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialDamageIterativeNonLocal); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.hh b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.hh index e7dbf8a34..3a8aee989 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.hh +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local.hh @@ -1,80 +1,80 @@ /** * @file material_damage_iterative_non_local.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief MaterialDamageIterativeNonLocal header for non-local damage * * @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_damage_iterative.hh" #include "material_damage_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_ITERATIVE_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_DAMAGE_ITERATIVE_NON_LOCAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material Damage Iterative Non local * * parameters in the material files : */ template class MaterialDamageIterativeNonLocal : public MaterialDamageNonLocal > { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialDamageNonLocal > MaterialDamageIterativeNonLocalParent; MaterialDamageIterativeNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialDamageIterativeNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual void computeNonLocalStresses(GhostType ghost_type); protected: void computeStress(ElementType type, GhostType ghost_type); void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost); private: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: InternalField grad_u_nl; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative_non_local_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_DAMAGE_ITERATIVE_NON_LOCAL_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local_inline_impl.cc index c264999b4..dda04de6c 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_non_local_inline_impl.cc @@ -1,87 +1,87 @@ /** * @file material_damage_iterative_non_local_inline_impl.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief MaterialDamageIterativeNonLocal inline function implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ -__END_AKANTU__ +} // namespace akantu #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialDamageIterativeNonLocal::MaterialDamageIterativeNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamageIterativeNonLocalParent(model, id), grad_u_nl("grad_u non local", *this) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->grad_u_nl.initialize(spatial_dimension*spatial_dimension); - this->model->getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), spatial_dimension*spatial_dimension); + this->model.getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), spatial_dimension*spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeNonLocal::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamageIterativeNonLocalParent::initMaterial(); - this->model->getNonLocalManager().nonLocalVariableToNeighborhood(grad_u_nl.getName(), this->name); + this->model.getNonLocalManager().nonLocalVariableToNeighborhood(grad_u_nl.getName(), this->name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeNonLocal::computeStress(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterativeNonLocal::computeNonLocalStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// compute the stress (based on the elastic law) MaterialDamage::computeStress(el_type, ghost_type); /// multiply the stress by (1-d) to get the effective stress Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeDamageAndStressOnQuad(sigma,*dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; /// compute the normalized equivalent stress this->computeNormalizedEquivalentStress(this->grad_u_nl(el_type, ghost_type), el_type, ghost_type); /// find the maximum this->norm_max_equivalent_stress = 0; this->findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_linear.cc b/extra_packages/extra-materials/src/material_damage/material_damage_linear.cc index 645a21946..63d0ba181 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_linear.cc +++ b/extra_packages/extra-materials/src/material_damage/material_damage_linear.cc @@ -1,76 +1,76 @@ /** * @file material_damage_linear.cc * * @author Marion Estelle Chambart * * * @brief Specialization of the material class for the damage material * * @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_linear.hh" #include "solid_mechanics_model.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialDamageLinear::MaterialDamageLinear(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), K("K", *this) { AKANTU_DEBUG_IN(); this->registerParam("Sigc", Sigc, 1e5, _pat_parsable, "Sigma Critique"); this->registerParam("Gc" , Gc , 2. , _pat_parsable, "Gc"); this->K.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); Epsmin = Sigc / this->E; Epsmax = 2 * Gc/ Sigc + Epsmin; this->K.setDefaultValue(Epsmin); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageLinear::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * K = this->K(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeStressOnQuad(grad_u, sigma, *dam, *K); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialDamageLinear); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_linear.hh b/extra_packages/extra-materials/src/material_damage/material_damage_linear.hh index 7f8d625d1..9f1143f56 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_linear.hh +++ b/extra_packages/extra-materials/src/material_damage/material_damage_linear.hh @@ -1,92 +1,92 @@ /** * @file material_damage_linear.hh * * @author Nicolas Richart * @author Marion Estelle Chambart * * * @brief Material isotropic elastic + linear softening * * @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_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_LINEAR_HH__ #define __AKANTU_MATERIAL_DAMAGE_LINEAR_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material liner damage * * parameters in the material files : * - Sigc : (default: 1e5) * - Gc : (default: 2) */ template class MaterialDamageLinear : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamageLinear(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialDamageLinear() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(Matrix & F, Matrix & sigma, Real & damage, Real &K); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// kind of toughness Real Gc; /// critical stress Real Sigc; /// damage internal variable InternalField K; Real Epsmin, Epsmax; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_linear_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_DAMAGE_LINEAR_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.cc b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.cc index b7f608b0f..d73dc78b1 100644 --- a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.cc +++ b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.cc @@ -1,203 +1,203 @@ /** * @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 "material_iterative_stiffness_reduction.hh" #include "solid_mechanics_model_RVE.hh" #include -__BEGIN_AKANTU__ +namespace akantu { template /* -------------------------------------------------------------------------- */ MaterialIterativeStiffnessReduction::MaterialIterativeStiffnessReduction(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamageIterative(model, id), eps_u("ultimate_strain", *this), D("tangent", *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->D.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIterativeStiffnessReduction::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamageIterative::initMaterial(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType ghost_type = *g; /// loop over all types in the material typedef ElementTypeMapArray:: type_iterator iterator; iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_regular); iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_regular); /// loop over all types in the filter for(; it != last_type; ++it) { ElementType el_type = *it; /// get the stiffness on each quad point Array::const_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); /// get the tangent of the tensile softening on each quad point Array::iterator D_it = this->D(el_type, ghost_type).begin(); Array::iterator D_end = this->D(el_type, ghost_type).end(); /// get the ultimate strain on each quad Array::iterator 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 void MaterialIterativeStiffnessReduction::computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// storage for the current stress Matrix sigma(spatial_dimension, spatial_dimension); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); /// iterators on the needed internal fields Array::const_scalar_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); Array::scalar_iterator dam_it = this->damage(el_type, ghost_type).begin(); Array::scalar_iterator equivalent_stress_it = this->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); /// loop over all the quadrature points and compute the equivalent stress for(;grad_u_it != grad_u_end; ++ grad_u_it) { /// compute the stress sigma.clear(); MaterialElastic::computeStressOnQuad(*grad_u_it, sigma, 0.); MaterialDamageIterative::computeDamageAndStressOnQuad(sigma, *dam_it); /// 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_it; } 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;; - 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); + 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); /// loop over all the elements for(; it != last_type; ++it) { ElementType el_type = *it; /// get iterators on the needed internal fields Array::const_scalar_iterator equivalent_stress_it = this->equivalent_stress(el_type, ghost_type).begin(); Array::const_scalar_iterator equivalent_stress_end = this->equivalent_stress(el_type, ghost_type).end(); Array::scalar_iterator dam_it = this->damage(el_type, ghost_type).begin(); Array::scalar_iterator reduction_it = this->reduction_step(el_type, ghost_type).begin(); Array::const_scalar_iterator eps_u_it = this->eps_u(el_type, ghost_type).begin(); Array::scalar_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); Array::const_scalar_iterator D_it = this->D(el_type, ghost_type).begin(); /// 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, ++D_it) { /// 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) ) * this->E * (*D_it)/ ((1. - (*dam_it) ) * this->E + (*D_it)); } nb_damaged_elements += 1; } } } } const SolidMechanicsModelRVE * rve_model = dynamic_cast(this->model); if (rve_model == NULL) { StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); } AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialIterativeStiffnessReduction); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.hh b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.hh index 0495115b4..8a3bdb015 100644 --- a/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.hh +++ b/extra_packages/extra-materials/src/material_damage/material_iterative_stiffness_reduction.hh @@ -1,105 +1,105 @@ /** * @file material_iterative_stiffness_reduction.hh * @author Aurelia Isabel Cuba Ramos * @date Thu Feb 18 15:25:05 2016 * * @brief Damage material with constant 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 "material_damage_iterative.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ITERATIVE_STIFFNESS_REDUCTION_HH__ #define __AKANTU_MATERIAL_ITERATIVE_STIFFNESS_REDUCTION_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material damage iterative * * parameters in the material files : * - Gfx * - h * - Sc */ /// Proposed by Rots and Invernizzi, 2004: Regularized sequentially linear // saw-tooth softening model (section 4.2) template class MaterialIterativeStiffnessReduction : public MaterialDamageIterative { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialIterativeStiffnessReduction(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialIterativeStiffnessReduction() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// init the material virtual void initMaterial(); ///compute the equivalent stress on each Gauss point (i.e. the max prinicpal stress) and normalize it by the tensile stiffness virtual void computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type = _not_ghost); /// update internal field damage virtual UInt updateDamage(); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the ultimate strain InternalField eps_u; /// the tangent of the tensile stress-strain softening InternalField D; /// fracture energy Real Gf; /// crack_band_width for normalization of fracture energy Real crack_band_width; /// the reduction constant (denoated by a in the paper of rots) Real reduction_constant; }; -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_ITERATIVE_STIFFNESS_REDUCTION_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 41669db1b..4084cd3bb 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,143 +1,143 @@ /** * @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__ -__BEGIN_AKANTU__ +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); /// 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; }; -__END_AKANTU__ +} // 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 e03617963..63e3d00f5 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,367 +1,367 @@ /** * @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 "solid_mechanics_model.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialOrthotropicDamageIterative:: MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, const ID & id) : Material(model, 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); 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); 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 StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&norm_max_equivalent_stress, 1, _so_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); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::const_matrix_iterator dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); Array::matrix_iterator damage_directions_it = this->damage_dir_vecs(el_type, _not_ghost).begin(this->spatial_dimension, this->spatial_dimension); Array::matrix_iterator 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); 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); 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); Array::matrix_iterator dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); Array::matrix_iterator damage_directions_it = this->damage_dir_vecs(el_type, _not_ghost).begin(this->spatial_dimension, this->spatial_dimension); Array::matrix_iterator 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; } StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { MaterialOrthotropicDamage::updateEnergies(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialOrthotropicDamageIterative); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.hh index 0938a1faa..56668d59e 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative.hh @@ -1,123 +1,123 @@ /** * @file material_orthtropic_damage_iterative.hh * * @author Aurelia Isabel Cuba Ramos * * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the class material orthotropic 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_orthotropic_damage.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material damage iterative * * parameters in the material files : * - Sc */ template class MaterialOrthotropicDamageIterative : public MaterialOrthotropicDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialOrthotropicDamageIterative() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// virtual void updateInternalParameters(); virtual void computeAllStresses(GhostType ghost_type = _not_ghost); /// update internal field damage UInt updateDamage(); /// update energies after damage has been updated virtual void updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_typ); virtual void onBeginningSolveStep(const AnalysisMethod & method) { }; virtual void onEndSolveStep(const AnalysisMethod & method) { }; protected: /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); ///compute the equivalent stress on each Gauss point (i.e. the max prinicpal stress) and normalize it by the tensile strength 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); inline void computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & root_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get max normalized equivalent stress AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// resistance to damage RandomInternalField Sc; /// internal field to store equivalent stress on each Gauss point InternalField equivalent_stress; /// internal field to store the direction of the current damage frame InternalField stress_dir; /// damage increment Real prescribed_dam; /// maximum equivalent stress Real norm_max_equivalent_stress; /// define damage threshold at which damage will be set to 1 Real dam_threshold; /// quadrature point with highest equivalent Stress IntegrationPoint q_max; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage_iterative_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local.hh index d2e4772b8..d69e6c8b8 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local.hh @@ -1,81 +1,81 @@ /** * @file material_orthotropic_damage_iterative_non_local.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief MaterialOrthotropicDamageIterativeNonLocal header for non-local damage * * @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_orthotropic_damage_iterative.hh" #include "material_damage_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_NON_LOCAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material Damage Iterative Non local * * parameters in the material files : */ template class MaterialOrthotropicDamageIterativeNonLocal : public MaterialDamageNonLocal > { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialDamageNonLocal > MaterialOrthotropicDamageIterativeNonLocalParent; MaterialOrthotropicDamageIterativeNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialOrthotropicDamageIterativeNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); protected: void computeStress(ElementType type, GhostType ghost_type); void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost); /// associate the non-local variables of the material to their neighborhoods virtual void nonLocalVariableToNeighborhood(); private: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: InternalField grad_u_nl; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage_iterative_non_local_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_NON_LOCAL_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local_inline_impl.cc index fc3d688aa..73e149644 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_iterative_non_local_inline_impl.cc @@ -1,129 +1,129 @@ /** * @file material_orthotropic_damage_iterative_non_local_inline_impl.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief MaterialOrthotropicDamageIterativeNonLocal inline function * implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ -__END_AKANTU__ +} // namespace akantu #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialOrthotropicDamageIterativeNonLocal::MaterialOrthotropicDamageIterativeNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialOrthotropicDamageIterativeNonLocalParent(model, id), grad_u_nl("grad_u non local", *this) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->grad_u_nl.initialize(spatial_dimension*spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterativeNonLocal::initMaterial() { AKANTU_DEBUG_IN(); - this->model->getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), spatial_dimension*spatial_dimension); + this->model.getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), spatial_dimension*spatial_dimension); MaterialOrthotropicDamageIterativeNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterativeNonLocal::computeStress(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterativeNonLocal::computeNonLocalStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialOrthotropicDamage::computeStress(el_type, ghost_type); 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); 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 { MaterialOrthotropicDamage::computeOneMinusD(one_minus_D_rotated, *damage_iterator); MaterialOrthotropicDamage::computeSqrtOneMinusD(one_minus_D_rotated, sqrt_one_minus_D_rotated); } this->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; this->computeNormalizedEquivalentStress(this->grad_u_nl(el_type, ghost_type), el_type, ghost_type); this->norm_max_equivalent_stress = 0; this->findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialOrthotropicDamageIterativeNonLocal::nonLocalVariableToNeighborhood() { - this->model->getNonLocalManager().nonLocalVariableToNeighborhood(grad_u_nl.getName(), this->name); + this->model.getNonLocalManager().nonLocalVariableToNeighborhood(grad_u_nl.getName(), this->name); } diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_non_local.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_non_local.hh index 745c5a7a6..7fd3e780b 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_non_local.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_non_local.hh @@ -1,100 +1,100 @@ /** * @file material_orthotropic_damage_non_local.hh * @author Aurelia Isabel Cuba Ramos * @date Sun Mar 22 21:10:27 2015 * * @brief interface for non local 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 "aka_common.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_NON_LOCAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { template class MaterialOrthotropicDamageNonLocal : public MaterialOrthotropicDamageLocal, public MaterialNonLocal { public: typedef MaterialNonLocal MaterialNonLocalParent; typedef MaterialOrthotropicDamageLocal MaterialOrthotropicDamageParent; MaterialOrthotropicDamageNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialOrthotropicDamageParent(model, id), MaterialNonLocalParent(model, id) { }; /* ------------------------------------------------------------------------ */ virtual void initMaterial() { MaterialOrthotropicDamageParent::initMaterial(); MaterialNonLocalParent::initMaterial(); } protected: /* -------------------------------------------------------------------------- */ virtual void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost) = 0; /* ------------------------------------------------------------------------ */ void computeNonLocalStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); - 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); + 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) { computeNonLocalStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } public: /* ------------------------------------------------------------------------ */ virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const { return MaterialNonLocalParent::getNbDataForElements(elements, tag) + MaterialOrthotropicDamageParent::getNbDataForElements(elements, tag); } virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { MaterialNonLocalParent::packElementData(buffer, elements, tag); MaterialOrthotropicDamageParent::packElementData(buffer, elements, tag); } virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { MaterialNonLocalParent::unpackElementData(buffer, elements, tag); MaterialOrthotropicDamageParent::unpackElementData(buffer, elements, tag); } }; -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_NON_LOCAL_HH__ */ 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 08e384041..a03c0efee 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,420 +1,420 @@ /** * @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" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> MaterialOrthotropicDamage::MaterialOrthotropicDamage(SolidMechanicsModel & model, const ID & id) : Material(model, 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); this->computePotentialEnergy(el_type, ghost_type); Array::matrix_iterator epsilon_p = this->gradu.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator sigma_p = this->stress.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::const_scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); Array::scalar_iterator ints = this->int_sigma(el_type, ghost_type).begin(); Array::scalar_iterator ed = this->dissipated_energy(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix delta_gradu_it(*gradu_it); 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); Array::matrix_iterator 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); Array::matrix_iterator 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)); // /// first row // tangent(0, 0) = (C(0,0)*S(0,0)*S(0,0) + C(1,0)*S(0,1)*S(0,1) - ((eta_effective*trace_D)/4. - 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) - ((eta_effective*trace_D)/4. - 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) - ((eta_effective*trace_D)/4. - 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) - ((eta_effective*trace_D)/4. - 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) { // /// set all terms to zero // first_term.clear(); // third_term.clear(); // /// hydrostatic sensitivity parameter // Real eta = 2.; // /// 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()); // for (UInt i = 0; i < this->spatial_dimension; ++i) { // for (UInt j = 0; j < this->spatial_dimension; ++j) // one_minus_D(i,j) *= second_term; // } // third_term.eye(1./this->spatial_dimension * sigma.trace() * (1 - eta/(this->spatial_dimension) * damage.trace())); // sigma = first_term - one_minus_D + third_term; /// hydrostatic sensitivity parameter // Real eta = 2.; /// 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()); // for (UInt i = 0; i < this->spatial_dimension; ++i) { // for (UInt j = 0; j < this->spatial_dimension; ++j) // one_minus_D(i,j) *= second_term; // } 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) { /// copy matrix // rotated.copy(to_rotate); // to_rotate.mul(rotated, damage_directions); // rotated.mul(damage_directions, to_rotate); 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) { /// copy matrix // rotated.copy(to_rotate); // to_rotate.mul(rotated, damage_directions); // rotated.mul(damage_directions, to_rotate); 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)); } /* -------------------------------------------------------------------------- */ // /* -------------------------------------------------------------------------- */ // template class Parent> // Real MaterialOrthotropicDamage::getDissipatedEnergy() const { // AKANTU_DEBUG_IN(); // Real de = 0.; -// const Mesh & mesh = this->model->getFEEngine().getMesh(); +// const Mesh & mesh = this->model.getFEEngine().getMesh(); // /// integrate the dissipated energy for each type of elements // Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost); // Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost); // for(; it != end; ++it) { -// de += this->model->getFEEngine().integrate(dissipated_energy(*it, _not_ghost), *it, +// de += this->model.getFEEngine().integrate(dissipated_energy(*it, _not_ghost), *it, // _not_ghost, this->element_filter(*it, _not_ghost)); // } // AKANTU_DEBUG_OUT(); // return de; // } // /* -------------------------------------------------------------------------- */ // template class Parent> // Real MaterialOrthotropicDamage::getEnergy(std::string type) { // if(type == "dissipated") return getDissipatedEnergy(); // else return Parent::getEnergy(type); // } // /* -------------------------------------------------------------------------- */ -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings.hh b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings.hh index 5bf4976ec..38b28c615 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings.hh +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings.hh @@ -1,146 +1,146 @@ /** * @file material_vreepeerlings.hh * * @author Cyprien Wolff * * * @brief Specialization of the material class for the VreePeerlings material * * @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_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VREEPEERLINGS_HH__ #define __AKANTU_MATERIAL_VREEPEERLINGS_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material vreepeerlings * * parameters in the material files : * - Kapaoi : (default: 0.0001) Initial threshold (of the equivalent strain) >= Crate * - Kapac : (default: 0.0002) Final threshold (of the equivalent strain) * - Arate : (default: 1.) Fitting parameter (must be close to 1 to do tend to 0 the stress in the damaged element) * - Brate : (default: 1.) This parameter determines the rate at which the damage grows * - Crate : (default: 0.0001) This parameter determines the rate at which the damage grows * - Kct : (default: 1.) Ratio between compressive and tensile strength */ template class MatParent = MaterialElastic> class MaterialVreePeerlings : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialDamage MaterialVreePeerlingsParent; MaterialVreePeerlings(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialVreePeerlings() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(Matrix & F, Matrix & sigma, Real & dam, Real & Equistrain_rate, Real & Equistrain, Real & Kapaq, Real dt, Matrix & strain_rate_vrpgls, Real & FullDam_ValStrain, Real & FullDam_ValStrain_rate, Real & Nb_damage); inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam, Real & Equistrain_rate, Real & Equistrain, Real & Kapaq, Real dt, Real & FullDam_Valstrain, Real & FullDam_Valstrain_rate, Real & Nb_damage); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Initial threshold (of the equivalent strain) (used in the initial step) Real Kapaoi; /// Final threshold (of the equivalent strain) (used in the initial step) Real Kapac; /// This parameter determines the rate at which the damage grows Real Arate; /// This parameter determines the rate at which the damage grows Real Brate; /// This parameter determines the rate at which the damage grows Real Crate; /// Ratio between compressive and tensile strength Real Kct; /// Randomness on Kapaoi Real Kapao_randomness; /// Kapa vector which contains the initial damage threshold RandomInternalField Kapa; /// Strain rate tensor to compute the rate dependent damage law InternalField strain_rate_vreepeerlings; /// Value of the equivalent strain when damage = 1 InternalField Full_dam_value_strain; /// Value of the equivalent strain rate when damage = 1 InternalField Full_dam_value_strain_rate; /// Count the number of times that the material is damaged to damage = 0 until damage = 1 InternalField Number_damage; /// Equivalent strain used to compute the damage evolution InternalField equi_strain; /// Equivalent strain rate used to compute the damage evolution InternalField equi_strain_rate; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_vreepeerlings_inline_impl.cc" #include "material_vreepeerlings_tmpl.hh" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_VREEPEERLINGS_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.cc b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.cc index 9688deaa9..72f5a0ba5 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.cc +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.cc @@ -1,28 +1,28 @@ /** * @file material_vreepeerlings_non_local.cc * * @author Nicolas Richart * @author Cyprien Wolff * * * @brief Specialization of the material class for the non-local Vree-Peerlings material * * @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_vreepeerlings_non_local.hh" #include "solid_mechanics_model.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ //INSTANTIATE_MATERIAL(MaterialVreePeerlingsNonLocal); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.hh b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.hh index 8f63c7fcc..d0979b9ce 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.hh +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local.hh @@ -1,91 +1,91 @@ /** * @file material_vreepeerlings_non_local.hh * * @author Nicolas Richart * @author Cyprien Wolff * * * @brief MaterialVreePeerlings header for non-local damage * * @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_vreepeerlings.hh" #include "material_damage_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material VreePeerlings Non local * * parameters in the material files : */ template class MatParent = MaterialElastic> class MaterialVreePeerlingsNonLocal : public MaterialDamageNonLocal > { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialVreePeerlings Parent; typedef MaterialDamageNonLocal MaterialVreePeerlingsNonLocalParent; MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialVreePeerlingsNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); /// constitutive law for all element of a type //void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// non local version of equivalent strain InternalField equi_strain_non_local; /// non local version of equivalent strain rate InternalField equi_strain_rate_non_local; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_vreepeerlings_non_local_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local_inline_impl.cc index 790b798e4..7a5f1f8ee 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_non_local_inline_impl.cc @@ -1,145 +1,145 @@ /** * @file material_vreepeerlings_non_local_inline_impl.cc * * @author Nicolas Richart * @author Cyprien Wolff * * * @brief Specialization of the material class for the non-local Vree-Peerlings material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class MatParent> MaterialVreePeerlingsNonLocal::MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialVreePeerlingsNonLocalParent(model, id), equi_strain_non_local("equi-strain_non_local", *this), equi_strain_rate_non_local("equi-strain-rate_non_local", *this) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->equi_strain_non_local.initialize(1); this->equi_strain_rate_non_local.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlingsNonLocal::initMaterial() { AKANTU_DEBUG_IN(); this->registerNonLocalVariable(this->equi_strain, this->equi_strain_non_local, 1); this->registerNonLocalVariable(this->equi_strain_rate, this->equi_strain_rate_non_local, 1); MaterialVreePeerlingsNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ //template class MatParent> //void MaterialVreePeerlingsNonLocal::computeStress(ElementType el_type, // GhostType ghost_type) { // AKANTU_DEBUG_IN(); // // Real * dam = this->damage(el_type, ghost_type).storage(); // Real * equi_straint = equi_strain(el_type, ghost_type).storage(); // Real * equi_straint_rate = equi_strain_rate(el_type, ghost_type).storage(); // Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); // Real * crit_strain = this->critical_strain(el_type, ghost_type).storage(); // Real * crit_strain_rate = this->critical_strain_rate(el_type, ghost_type).storage(); // Real * rdr_damage = this->recorder_damage(el_type, ghost_type).storage(); // Real * nb_damage = this->number_damage(el_type, ghost_type).storage(); -// Real dt = this->model->getTimeStep(); +// Real dt = this->model.getTimeStep(); // // Vector & elem_filter = this->element_filter(el_type, ghost_type); -// Vector & velocity = this->model->getVelocity(); +// Vector & velocity = this->model.getVelocity(); // Vector & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type); // // -// this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate_vrplgs, +// this->model.getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate_vrplgs, // spatial_dimension, // el_type, ghost_type, &elem_filter); // // Vector::iterator strain_rate_vrplgs_it = // strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension); // // // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); // // types::RMatrix & strain_rate = *strain_rate_vrplgs_it; // // // // MaterialVreePeerlings::computeStressOnQuad(grad_u, sigma, // *dam, // *equi_straint, // *equi_straint_rate, // *Kapaq, // dt, // strain_rate, // *crit_strain, // *crit_strain_rate, // *rdr_damage, // *nb_damage); // ++dam; // ++equi_straint; // ++equi_straint_rate; // ++Kapaq; // ++strain_rate_vrplgs_it; // ++crit_strain; // ++crit_strain_rate; // ++rdr_damage; // ++nb_damage; // // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; // // AKANTU_DEBUG_OUT(); //} // /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlingsNonLocal::computeNonLocalStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); Real * equi_strain_nl = this->equi_strain_non_local(el_type, ghost_type).storage(); Real * equi_strain_rate_nl = this->equi_strain_rate_non_local(el_type, ghost_type).storage(); //Real * equi_strain_rate_nl = this->equi_strain_rate(el_type, ghost_type).storage(); - Real dt = this->model->getTimeStep(); + Real dt = this->model.getTimeStep(); Real * FullDam_Valstrain = this->Full_dam_value_strain(el_type, ghost_type).storage(); Real * FullDam_Valstrain_rate = this->Full_dam_value_strain_rate(el_type, ghost_type).storage(); Real * Nb_damage = this->Number_damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeDamageAndStressOnQuad(sigma, *dam, *equi_strain_nl, *equi_strain_rate_nl, *Kapaq, dt, *FullDam_Valstrain, *FullDam_Valstrain_rate, *Nb_damage); ++dam; ++equi_strain_nl; ++equi_strain_rate_nl; ++Kapaq; ++FullDam_Valstrain; ++FullDam_Valstrain_rate; ++Nb_damage; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh index f0406a03a..b211f6077 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh @@ -1,112 +1,112 @@ /** * @file material_vreepeerlings_tmpl.hh * * @author Cyprien Wolff * * * @brief Specialization of the material class for the VreePeerlings material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ template class MatParent> MaterialVreePeerlings::MaterialVreePeerlings(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialVreePeerlingsParent(model, id), Kapa("Kapa", *this), strain_rate_vreepeerlings("strain-rate-vreepeerlings", *this), Full_dam_value_strain("fulldam-valstrain", *this), Full_dam_value_strain_rate("fulldam-valstrain-rate", *this), Number_damage("number-damage", *this), equi_strain("equi-strain", *this), equi_strain_rate("equi-strain-rate", *this) { AKANTU_DEBUG_IN(); this->registerParam("Kapaoi" , Kapaoi , 0.0001, _pat_parsable); this->registerParam("Kapac" , Kapac , 0.0002, _pat_parsable); this->registerParam("Arate" , Arate , 0. , _pat_parsable); this->registerParam("Brate" , Brate , 1. , _pat_parsable); this->registerParam("Crate" , Brate , 1. , _pat_parsable); this->registerParam("Kct" , Kct , 1. , _pat_parsable); this->registerParam("Kapao_randomness", Kapao_randomness, 0. , _pat_parsable); this->Kapa.initialize(1); this->equi_strain.initialize(1); this->equi_strain_rate.initialize(1); this->Full_dam_value_strain.initialize(1); this->Full_dam_value_strain_rate.initialize(1); this->Number_damage.initialize(1); this->strain_rate_vreepeerlings.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlings::initMaterial() { AKANTU_DEBUG_IN(); MaterialVreePeerlingsParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlings::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialVreePeerlingsParent::computeStress(el_type, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); Real * equi_straint = equi_strain(el_type, ghost_type).storage(); Real * equi_straint_rate = equi_strain_rate(el_type, ghost_type).storage(); Real * Kapaq = Kapa(el_type, ghost_type).storage(); Real * FullDam_Valstrain = Full_dam_value_strain(el_type, ghost_type).storage(); Real * FullDam_Valstrain_rate = Full_dam_value_strain_rate(el_type, ghost_type).storage(); Real * Nb_damage = Number_damage(el_type, ghost_type).storage(); - Real dt = this->model->getTimeStep(); + Real dt = this->model.getTimeStep(); Array & elem_filter = this->element_filter(el_type, ghost_type); - Array & velocity = this->model->getVelocity(); + Array & velocity = this->model.getVelocity(); Array & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type); - this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_vrplgs, + this->model.getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_vrplgs, spatial_dimension, el_type, ghost_type, elem_filter); Array::matrix_iterator strain_rate_vrplgs_it = strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & strain_rate = *strain_rate_vrplgs_it; computeStressOnQuad(grad_u, sigma, *dam, *equi_straint, *equi_straint_rate, *Kapaq, dt, strain_rate, *FullDam_Valstrain, *FullDam_Valstrain_rate, *Nb_damage); ++dam; ++equi_straint; ++equi_straint_rate; ++Kapaq; ++strain_rate_vrplgs_it; ++FullDam_Valstrain; ++FullDam_Valstrain_rate; ++Nb_damage; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } 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 418d4643f..cc4a6e5c6 100644 --- a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc +++ b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc @@ -1,104 +1,104 @@ /** * @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" -__BEGIN_AKANTU__ +namespace 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->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->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(); } INSTANTIATE_MATERIAL(MaterialViscoPlastic); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh index 8f2aa293e..4559397a4 100644 --- a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh +++ b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh @@ -1,104 +1,104 @@ /** * @file material_viscoplastic.hh * * @author Ramin Aghababaei * * * @brief Specialization of the material class for * MaterialLinearIsotropicHardening to include viscous effects (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 "aka_common.hh" #include "material_plastic.hh" #include "aka_voigthelper.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VISCOPLASTIC_HH__ #define __AKANTU_MATERIAL_VISCOPLASTIC_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material plastic isotropic * * parameters in the material files : * - h : Hardening parameter (default: 0) * - sigmay : Yield stress * - rate : Rate sensitivity * - edot0 : Reference strain rate * * - ts: Time step */ template class MaterialViscoPlastic : public MaterialPlastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialViscoPlastic(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); protected: inline void computeStressOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, Real & iso_hardening) const; inline void computeTangentModuliOnQuad(Matrix & tangent, const Matrix & grad_u, const Matrix & previous_grad_u, const Matrix & sigma_tensor, const Matrix & previous_sigma_tensor, const Real & iso_hardening) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Rate sensitivity component (rate) Real rate; /// Reference strain rate (edot0) Real edot0; /// Time step (ts) Real ts; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_viscoplastic_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_VISCOPLASTIC_HH__ */ diff --git a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc index 738cfbdcd..d28ca2c11 100644 --- a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc +++ b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc @@ -1,128 +1,128 @@ /** * @file material_stiffness_proportional.cc * * @author David Simon Kammer * @author Nicolas Richart * * * @brief Special. of the material class for the caughey viscoelastic material * * @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_stiffness_proportional.hh" #include "solid_mechanics_model.hh" -__BEGIN_AKANTU__ +namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialStiffnessProportional::MaterialStiffnessProportional(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), stress_viscosity("stress_viscosity", *this), stress_elastic("stress_elastic", *this) { AKANTU_DEBUG_IN(); this->registerParam("Alpha", alpha, 0., ParamAccessType(_pat_parsable | _pat_modifiable), "Artificial viscous ratio"); this->stress_viscosity.initialize(spatial_dimension * spatial_dimension); this->stress_elastic .initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::initMaterial() { AKANTU_DEBUG_IN(); MaterialElastic::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = this->element_filter (el_type, ghost_type); Array & stress_visc = stress_viscosity(el_type, ghost_type); Array & stress_el = stress_elastic (el_type, ghost_type); MaterialElastic::computeStress(el_type, ghost_type); - Array & velocity = this->model->getVelocity(); + Array & velocity = this->model.getVelocity(); Array strain_rate(0, spatial_dimension * spatial_dimension, "strain_rate"); - this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate, + this->model.getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate, spatial_dimension, el_type, ghost_type, elem_filter); Array::matrix_iterator strain_rate_it = strain_rate.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_visc_it = stress_visc.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & grad_v = *strain_rate_it; Matrix & sigma_visc = *stress_visc_it; Matrix & sigma_el = *stress_el_it; MaterialElastic::computeStressOnQuad(grad_v, sigma_visc); sigma_visc *= alpha; sigma_el.copy(sigma); sigma += sigma_visc; ++strain_rate_it; ++stress_visc_it; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Array & stress_el = stress_elastic(el_type, ghost_type); Array::matrix_iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); Real * epot = this->potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & sigma_el = *stress_el_it; MaterialElastic::computePotentialEnergyOnQuad(grad_u, sigma_el, *epot); epot++; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialStiffnessProportional); -__END_AKANTU__ +} // namespace akantu diff --git a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh index 3d5699b0d..43779d0b5 100644 --- a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh +++ b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh @@ -1,102 +1,102 @@ /** * @file material_stiffness_proportional.hh * * @author David Simon Kammer * @author Nicolas Richart * * * @brief Material isotropic visco-elastic with viscosity proportional to the * stiffness * * @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_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH__ #define __AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH__ -__BEGIN_AKANTU__ +namespace akantu { /** * Material visco-elastic @f[\sigma = E\epsilon + \alpha E* \frac{d\epsilon}{dt}@f] * it can be seen as a Kelvin-Voigt solid with @f[\eta = \alpha E @f] * * The material satisfies the Caughey condition, the visco-elastic solid has the * same eigen-modes as the elastic one. (T.K. Caughey 1960 - Journal of Applied * Mechanics 27, 269-271. Classical normal modes in damped linear systems.) * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) * - alpha : viscous ratio */ template class MaterialStiffnessProportional : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialStiffnessProportional(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialStiffnessProportional() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for all elements virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for a given quadrature point //inline void computeStress(Real * F, Real * sigma); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// stress due to viscosity InternalField stress_viscosity; /// stress due to elasticity InternalField stress_elastic; /// viscous ratio Real alpha; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_elastic_caughey_inline_impl.cc" -__END_AKANTU__ +} // namespace akantu #endif /* __AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH__ */