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 65a302b74..fbf0ccf69 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,196 @@ /** * @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" +#include "communicator.hh" +#include "solid_mechanics_model_RVE.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template +template MaterialFE2::MaterialFE2(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), Parent(model, id), - C("material_stiffness", *this) - { + const ID & 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 +MaterialFE2::~MaterialFE2() = default; /* -------------------------------------------------------------------------- */ -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::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 +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(); + /// create a Mesh and SolidMechanicsModel on each integration point of the + /// material + const auto & comm = this->model.getMesh().getCommunicator(); 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(); + auto C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size); + + for (auto && data : + enumerate(make_view(C(this->el_type), voigt_h::size, voigt_h::size))) { + auto q = std::get<0>(data); + auto & C = std::get<1>(data); + + meshes.emplace_back(std::make_unique( + spatial_dimension, "RVE_mesh_" + std::to_string(prank), q + 1)); + + auto & mesh = *meshes.back(); + mesh.read(mesh_file); + + RVEs.emplace_back(std::make_unique( + mesh, true, this->nb_gel_pockets, _all_dimensions, + "SMM_RVE_" + std::to_string(prank), q + 1)); + + auto & RVE = *RVEs.back(); + RVE.initFull(); + /// compute intial stiffness of the RVE - (*RVE_it)->homogenizeStiffness(*C_it); + RVE.homogenizeStiffness(C); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void MaterialFE2::computeStress(ElementType el_type, - GhostType ghost_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(); + 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); + 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) { + 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) { + 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; + 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) { +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); + 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(); +template +void MaterialFE2::advanceASR( + const Matrix & prestrain) { + AKANTU_DEBUG_IN(); - 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 (auto && data : + zip(RVEs, make_view(this->gradu(this->el_type), spatial_dimension, + spatial_dimension), + make_view(this->eigengradu(this->el_type), spatial_dimension, + spatial_dimension), + make_view(this->C(this->el_type), voigt_h::size, voigt_h::size))) { + auto & RVE = *(std::get<0>(data)); - 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); + /// apply boundary conditions based on the current macroscopic displ. + /// gradient + RVE.applyBoundaryConditions(std::get<1>(data)); /// advance the ASR in every RVE - (*RVE_it)->advanceASR(prestrain); + RVE.advanceASR(prestrain); /// compute the average eigen_grad_u - (*RVE_it)->homogenizeEigenGradU(*eigen_gradu_it); + RVE.homogenizeEigenGradU(std::get<2>(data)); /// compute the new effective stiffness of the RVE - (*RVE_it)->homogenizeStiffness(*C_it); - + RVE.homogenizeStiffness(std::get<3>(data)); } + AKANTU_DEBUG_OUT(); } - -INSTANTIATE_MATERIAL(MaterialFE2); - +INSTANTIATE_MATERIAL(material_FE2, MaterialFE2); } // 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 0ddd75ba7..64d6c4fa4 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,122 @@ /** * @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. + * 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__ +namespace akantu { +class SolidMechanicsModelRVE; +} + namespace akantu { /* -------------------------------------------------------------------------- */ -/// /!\ This material works ONLY for meshes with a single element type!!!!! +/// /!\ This material works ONLY for meshes with a single element type!!!!! /* -------------------------------------------------------------------------- */ - /** * MaterialFE2 * * parameters in the material files : - * - mesh_file + * - 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 */ - /* ------------------------------------------------------------------------ */ +// template +// class MaterialFE2 : public virtual Material { +// /* ------------------------------------------------------------------------ +// */ +// /* Constructors/Destructors */ +// /* ------------------------------------------------------------------------ +// */ +template class MaterialFE2 : public MaterialThermal { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ private: typedef MaterialThermal Parent; -// Emil - 27.01.2018 - public: - + // 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); + 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); + 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; + std::vector> RVEs; /// Meshes for all RVEs - std::vector meshes; + std::vector> meshes; - /// the element type of the associated mesh (this material handles only one type!!) + /// 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) + /// 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" } // 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 aefff6a96..4eb79aec2 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,552 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" +#include "element_group.hh" #include "material_damage_iterative.hh" -#ifdef AKANTU_USE_MUMPS -#include "solver_mumps.hh" -#endif -#ifdef AKANTU_USE_PETSC -#include "solver_petsc.hh" -#include "petsc_matrix.hh" -#endif +#include "node_group.hh" +#include "non_linear_solver.hh" +#include "parser.hh" /* -------------------------------------------------------------------------- */ -__BEGIN_AKANTU__ + +namespace akantu { /* -------------------------------------------------------------------------- */ -SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, - UInt nb_gel_pockets, UInt dim, - const ID & id, const MemoryID & memory_id) : - SolidMechanicsModel(mesh, dim, id, memory_id), - volume(0.), - use_RVE_mat_selector(use_RVE_mat_selector), - static_communicator_dummy(StaticCommunicator::getStaticCommunicatorDummy()), - nb_gel_pockets(nb_gel_pockets), - nb_dumps(0) { +SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, + bool use_RVE_mat_selector, + UInt nb_gel_pockets, UInt dim, + const ID & id, + const MemoryID & memory_id) + : SolidMechanicsModel(mesh, dim, id, memory_id), volume(0.), + use_RVE_mat_selector(use_RVE_mat_selector), + nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// 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 auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); + bottom_nodes.insert(bottom.begin(), bottom.end()); + + const auto & left = mesh.getElementGroup("left").getNodeGroup(); + left_nodes.insert(left.begin(), left.end()); - const ElementGroup & left = mesh.getElementGroup("left"); - left_nodes.insert( left.node_begin(), left.node_end() ); + // /// enforce periodicity on the displacement fluctuations + // auto surface_pair_1 = std::make_pair("top", "bottom"); + // auto surface_pair_2 = std::make_pair("right", "left"); + // SurfacePairList surface_pairs_list; + // surface_pairs_list.push_back(surface_pair_1); + // surface_pairs_list.push_back(surface_pair_2); + // TODO: To Nicolas correct the PBCs + // this->setPBC(surface_pairs_list); - /// 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(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -SolidMechanicsModelRVE::~SolidMechanicsModelRVE() { - delete static_communicator_dummy; -} +SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ -void SolidMechanicsModelRVE::initFull(const ModelOptions & options) { +void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); - SolidMechanicsModel::initFull(options); + + auto options_cp(options); + options_cp.analysis_method = AnalysisMethod::_static; + + SolidMechanicsModel::initFullImpl(options); this->initMaterials(); + auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); + /// 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); + for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { + Array Volume(this->mesh.getNbElement(element_type) * + fem.getNbIntegrationPoints(element_type), + 1, 1.); + this->volume = fem.integrate(Volume, element_type); } - std::cout << "The volume of the RVE is " << this->volume << std::endl; + 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->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->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("internal_forces"); this->dump(); - this->nb_dumps +=1; - AKANTU_DEBUG_OUT(); + this->nb_dumps += 1; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModelRVE::applyBoundaryConditions(const Matrix & displacement_gradient) { +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 + /// 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); + 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); + 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); + 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); + 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(); + const auto & position = mesh.getNodes(); + const auto & lower_bounds = mesh.getLowerBounds(); + const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); - for (UInt n = 0; n < mesh.getNbNodes(); ++n) { + for (auto && data : enumerate(make_view(position, spatial_dimension))) { + auto node = std::get<0>(data); + const auto & X = std::get<1>(data); + + auto distance = X.distance(lower_bounds); // node 1 - if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && - Math::are_float_equal(position(n, 1), lower_bounds(1))) { - corner_nodes(0) = n; + if (Math::are_float_equal(distance, 0)) { + corner_nodes(0) = node; } // 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; + else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && + Math::are_float_equal(X(_y), lower_bounds(_y))) { + corner_nodes(1) = node; } // node 3 - else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) && - Math::are_float_equal(position(n, 1), upper_bounds(1))) { - corner_nodes(2) = n; + else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && + Math::are_float_equal(X(_y), upper_bounds(_y))) { + corner_nodes(2) = node; } // node 4 - else if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && - Math::are_float_equal(position(n, 1), upper_bounds(1))) { - corner_nodes(3) = n; + else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && + Math::are_float_equal(X(_y), upper_bounds(_y))) { + corner_nodes(3) = node; } } - for (UInt i = 0; i < corner_nodes.getSize(); ++i) { + for (UInt i = 0; i < corner_nodes.size(); ++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) + for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { + Array & prestrain_vect = + const_cast &>(this->getMaterial("gel").getInternal( + "eigen_grad_u")(element_type)); + auto prestrain_it = + prestrain_vect.begin(spatial_dimension, spatial_dimension); + auto prestrain_end = + prestrain_vect.end(spatial_dimension, spatial_dimension); + + for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } - /// advance the damage - MaterialDamageIterative<2> & mat_paste = dynamic_cast & >(*this->materials[1]); - MaterialDamageIterative<2> & mat_aggregate = dynamic_cast & >(*this->materials[0]); + 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 + + auto & solver = this->getNonLinearSolver(); + solver.set("max_iterations", 2); + solver.set("threshold", 1e-6); + solver.set("convergence_type", _scc_solution); + + do { + this->solveStep(); + + /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); - + nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) nb_damaged_elements = mat_aggregate.updateDamage(); else if (max_eq_stress_aggregate < max_eq_stress_paste) nb_damaged_elements = mat_paste.updateDamage(); else - nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); + nb_damaged_elements = + (mat_paste.updateDamage() + mat_aggregate.updateDamage()); - std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; + 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) { +Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, + const ID & field_type) { AKANTU_DEBUG_IN(); - FEEngine * fem = this->fems["SolidMechanicsFEEngine"]; + auto & fem = this->getFEEngine("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; + for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { 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) + for (UInt m = 0; m < this->materials.size(); ++m) { + const auto & stress_vec = this->materials[m]->getStress(element_type); + const auto & elem_filter = + this->materials[m]->getElementFilter(element_type); + Array int_stress_vec(elem_filter.size(), + spatial_dimension * spatial_dimension, + "int_of_stress"); + + fem.integrate(stress_vec, int_stress_vec, + spatial_dimension * spatial_dimension, element_type, + _not_ghost, elem_filter); + + for (UInt k = 0; k < elem_filter.size(); ++k) + average += int_stress_vec( + k, row_index * spatial_dimension + col_index); // 3 is the value + // for the yy (in + // 3D, the value is + // 4) } - } - - else if (field_type == "strain") { - for(UInt m = 0; m < this->materials.size(); ++m) { - const 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 == "strain") { + for (UInt m = 0; m < this->materials.size(); ++m) { + const auto & gradu_vec = this->materials[m]->getGradU(element_type); + const auto & elem_filter = + this->materials[m]->getElementFilter(element_type); + Array int_gradu_vec(elem_filter.size(), + spatial_dimension * spatial_dimension, + "int_of_gradu"); + + fem.integrate(gradu_vec, int_gradu_vec, + spatial_dimension * spatial_dimension, element_type, + _not_ghost, elem_filter); + + for (UInt k = 0; k < elem_filter.size(); ++k) + /// averaging is done only for normal components, so stress and strain + /// are equal + average += + 0.5 * + (int_gradu_vec(k, row_index * spatial_dimension + col_index) + + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } - } - - else if (field_type == "eigen_grad_u") { - for(UInt m = 0; m < this->materials.size(); ++m) { - const 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 if (field_type == "eigen_grad_u") { + for (UInt m = 0; m < this->materials.size(); ++m) { + const auto & eigen_gradu_vec = + this->materials[m]->getInternal("eigen_grad_u")(element_type); + const auto & elem_filter = + this->materials[m]->getElementFilter(element_type); + Array int_eigen_gradu_vec(elem_filter.size(), + spatial_dimension * spatial_dimension, + "int_of_gradu"); + + fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, + spatial_dimension * spatial_dimension, element_type, + _not_ghost, elem_filter); + + for (UInt k = 0; k < elem_filter.size(); ++k) + /// averaging is done only for normal components, so stress and strain + /// are equal + average += + int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } - } - - else + } else { AKANTU_DEBUG_ERROR("Averaging not implemented for this field!!!"); + } } - return average/this->volume; + + 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!!!"); + 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) + for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { + auto & prestrain_vect = + const_cast &>(this->getMaterial("gel").getInternal( + "eigen_grad_u")(element_type)); + auto prestrain_it = + prestrain_vect.begin(spatial_dimension, spatial_dimension); + auto prestrain_end = + prestrain_vect.end(spatial_dimension, spatial_dimension); + + for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); - /// save the damage state before 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); - + /// save the damage state before filling up cracks + // ElementTypeMapReal saved_damage("saved_damage"); + // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = + // 0); + // this->fillCracks(saved_damage); + /// virtual test 1: - H(0,0) = 0.01; + H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); - H(1,1) = 0.01; + H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); - H(0,1) = 0.01; + H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks - //this->drainCracks(saved_damage); + // 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) { +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"); + auto & solver = this->getNonLinearSolver("static"); + solver.set("max_iterations", 2); + solver.set("threshold", 1e-6); + solver.set("convergence_type", _scc_solution); + this->solveStep("static"); /// 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"); + 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) { +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"); + 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 (!are_materials_instantiated) + instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); - Real bottom = lowerBounds(1); + Real bottom = lowerBounds(1); Real top = upperBounds(1); - Real box_size = std::abs(top-bottom); + 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; + auto tmp = std::make_shared(*this, box_size, "gel", + this->nb_gel_pockets, eps); + tmp->setFallback(material_selector); + material_selector = tmp; } 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 - } +void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { + const auto & mat_gel = this->getMaterial("gel"); + Real E_gel = mat_gel.get("E"); + Real E_homogenized = 0.; - 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(); -} + for (auto && mat : materials) { + if (mat->getName() == "gel" || mat->getName() == "FE2_mat") + continue; -/* -------------------------------------------------------------------------- */ -void SolidMechanicsModelRVE::initArrays() { - AKANTU_DEBUG_IN(); + Real E = mat->get("E"); + auto & damage = mat->getInternal("damage"); - 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(); + for (auto && type : damage.elementTypes()) { + const auto & elem_filter = mat->getElementFilter(type); + auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); - AKANTU_DEBUG_OUT(); -} + auto sav_dam_it = + make_view(saved_damage(type), nb_integration_point).begin(); + for (auto && data : + zip(elem_filter, make_view(damage(type), nb_integration_point))) { + auto el = std::get<0>(data); + auto & dam = std::get<1>(data); + Vector sav_dam = sav_dam_it[el]; -/* -------------------------------------------------------------------------- */ -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); + sav_dam = dam; + + for (auto q : arange(dam.size())) { + E_homogenized = (E_gel - E) * dam(q) + E; + dam(q) = 1. - (E_homogenized / E); + } } } } } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModelRVE::drainCracks(const ElementTypeMapReal & saved_damage) { - 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") +void SolidMechanicsModelRVE::drainCracks( + const ElementTypeMapReal & saved_damage) { + for (auto && mat : materials) { + 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)); - } + auto & damage = mat->getInternal("damage"); + + for (auto && type : damage.elementTypes()) { + const auto & elem_filter = mat->getElementFilter(type); + auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); + + auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); + for (auto && data : + zip(elem_filter, make_view(damage(type), nb_integration_point))) { + auto el = std::get<0>(data); + auto & dam = std::get<1>(data); + Vector sav_dam = sav_dam_it[el]; + + dam = sav_dam; } } } } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh index d83557c30..c38055a1a 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,230 @@ /** * @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 "solid_mechanics_model.hh" #include /* -------------------------------------------------------------------------- */ -__BEGIN_AKANTU__ - +namespace 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); + 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)); +protected: + void initFullImpl(const ModelOptions & option) override; /// initialize the materials - virtual void initMaterials(); + void initMaterials() override; +public: /// apply boundary contions based on macroscopic deformation gradient - virtual void applyBoundaryConditions(const Matrix & displacement_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); + 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); + inline void unpackData(CommunicationBuffer & buffer, + const Array & index, + const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ - /* Accessors */ + /* 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 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 */ + /* 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) { + const Array & index, + const 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)]); + auto disp_it = displacement->begin(spatial_dimension); + + for (auto node : index) { + Vector current_disp(disp_it[node]); + + // if node is at the bottom, u_bottom = u_top +u_2 -u_3 + if (bottom_nodes.count(node)) { + 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(node)) { + current_disp += Vector(disp_it[corner_nodes(3)]); + current_disp -= Vector(disp_it[corner_nodes(2)]); + } } } } /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ -class GelMaterialSelector : public MeshDataMaterialSelector { +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) { + 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); + Element el{_triangle_3, 0, _not_ghost}; + UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); + Array barycenter(nb_element, spatial_dimension); + + for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { + el.element = std::get<0>(data); + auto & bary = std::get<1>(data); + mesh.getBarycenter(el, bary); } /// 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; + continue; checked_baries.insert(bary_id); el.element = bary_id; if (MeshDataMaterialSelector::operator()(el) == 1) - continue; /// element belongs to paste + 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) { + 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 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; }; - } // 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_damage_iterative.cc b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.cc index e7387097b..99f1295b2 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,259 @@ /** * @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" +#include "communicator.hh" +#include "data_accessor.hh" +/* -------------------------------------------------------------------------- */ 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) { +template +MaterialDamageIterative::MaterialDamageIterative( + SolidMechanicsModel & model, const ID & id) + : MaterialDamage(model, id), + Sc("Sc", *this), reduction_step("damage_step", *this), + equivalent_stress("equivalent_stress", *this), max_reductions(0), + norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); - this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold"); - this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "prescribed damage" ); - this->registerParam("dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1" ); - this->registerParam("dam_tolerance", dam_tolerance, 0.01, _pat_parsable | _pat_modifiable, "damage tolerance to decide if quadrature point will be damageed" ); - this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" ); - this->registerParam("max_reductions", max_reductions, UInt(10), _pat_parsable | _pat_modifiable, "max reductions"); - - this->use_previous_stress = true; - this->use_previous_gradu = true; + this->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) { +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::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); + 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) { + for (; grad_u_it != grad_u_end; ++grad_u_it) { sigma.clear(); - MaterialElastic::computeStressOnQuad(*grad_u_it, sigma, 0.); - computeDamageAndStressOnQuad(sigma,*dam); + 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); + *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) { +template +void MaterialDamageIterative::computeAllStresses( + GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress - if(ghost_type==_not_ghost) + 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); + auto 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); + const auto & comm = this->model.getMesh().getCommunicator(); + comm.allReduce(norm_max_equivalent_stress, SynchronizerOperation::_max); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialDamageIterative::findMaxNormalizedEquivalentStress(ElementType el_type, - GhostType ghost_type) { +template +void MaterialDamageIterative:: + findMaxNormalizedEquivalentStress(ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); - if(ghost_type==_not_ghost) { + 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) + // 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; + 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) { +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); + computeDamageAndStressOnQuad(sigma, *dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type); + computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, + ghost_type); norm_max_equivalent_stress = 0; findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template UInt MaterialDamageIterative::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., - "Your prescribed damage must be greater than zero"); + "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;; + /// 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) { + 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; - } + 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); + auto * rve_model = + dynamic_cast(&this->model); if (rve_model == NULL) { - StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); - comm.allReduce(&nb_damaged_elements, 1, _so_sum); + const auto & comm = this->model.getMesh().getCommunicator(); + comm.allReduce(nb_damaged_elements, SynchronizerOperation::_sum); } + AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ -template -void MaterialDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { +template +void MaterialDamageIterative::updateEnergiesAfterDamage( + ElementType el_type, GhostType ghost_type) { MaterialDamage::updateEnergies(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ -INSTANTIATE_MATERIAL(MaterialDamageIterative); - +INSTANTIATE_MATERIAL(damage_iterative, MaterialDamageIterative); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh index 88d620503..acd902544 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" +#include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ #define __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ namespace akantu { /** * Material damage iterative * * parameters in the material files : - * - Sc + * - Sc */ -template +template class MaterialDamageIterative : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MaterialDamageIterative(SolidMechanicsModel & model, const ID & id = ""); - virtual ~MaterialDamageIterative() {}; + 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); + 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) { }; + virtual void updateEnergiesAfterDamage(ElementType el_type, + GhostType ghost_typ); - ///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); + /// 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); + 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); + virtual void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost); - inline void computeDamageAndStressOnQuad(Matrix & sigma, - Real & dam); + inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ - inline UInt getNbDataForElements(const Array & elements, - SynchronizationTag tag) const; + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; - inline void packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const; + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; - inline void unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag); + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get max normalized equivalent stress AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real); - /// get a non-const reference to the max normalized equivalent stress - AKANTU_GET_MACRO_NOT_CONST(NormMaxEquivalentStress, norm_max_equivalent_stress, Real &); + /// 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 + /// 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; + /// define damage threshold at which damage will be set to 1 + Real dam_threshold; /// maximum damage value Real max_damage; }; +} // namespace akantu + /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative_inline_impl.cc" -} // namespace akantu - #endif /* __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_inline_impl.cc b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_inline_impl.cc index 14430ee19..283b95a47 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative_inline_impl.cc @@ -1,80 +1,92 @@ /** * @file material_damage_iterative_inline_impl.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Implementation of inline functions of the material damage iterative * * @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" +/* -------------------------------------------------------------------------- */ + +namespace akantu { /* -------------------------------------------------------------------------- */ -template +template inline void -MaterialDamageIterative::computeDamageAndStressOnQuad(Matrix & sigma, Real & dam) { - sigma *= 1-dam; +MaterialDamageIterative::computeDamageAndStressOnQuad( + Matrix & sigma, Real & dam) { + sigma *= 1 - dam; } /* -------------------------------------------------------------------------- */ -template -UInt MaterialDamageIterative::updateDamage(UInt quad_index, const Real eq_stress, const ElementType & el_type, const GhostType & ghost_type) { +template +UInt MaterialDamageIterative::updateDamage( + UInt quad_index, const Real /*eq_stress*/, const ElementType & el_type, + const GhostType & ghost_type) { AKANTU_DEBUG_ASSERT(prescribed_dam > 0., - "Your prescribed damage must be greater than zero"); - + "Your prescribed damage must be greater than zero"); Array & dam = this->damage(el_type, ghost_type); - Real & dam_on_quad = dam(quad_index); + Real & dam_on_quad = dam(quad_index); /// check if damage occurs - if (equivalent_stress(el_type, ghost_type)(quad_index) >= (1-dam_tolerance) * norm_max_equivalent_stress) { + if (equivalent_stress(el_type, ghost_type)(quad_index) >= + (1 - dam_tolerance) * norm_max_equivalent_stress) { if (dam_on_quad < dam_threshold) - dam_on_quad +=prescribed_dam; - else dam_on_quad = max_damage; + dam_on_quad += prescribed_dam; + else + dam_on_quad = max_damage; return 1; } return 0; } /* -------------------------------------------------------------------------- */ -template -inline UInt MaterialDamageIterative::getNbDataForElements(const Array & elements, - SynchronizationTag tag) const { +template +inline UInt MaterialDamageIterative::getNbData( + const Array & elements, const SynchronizationTag & tag) const { if (tag == _gst_user_2) { return sizeof(Real) * this->getModel().getNbIntegrationPoints(elements); } - return MaterialDamage::getNbDataForElements(elements, tag); + return MaterialDamage::getNbData(elements, tag); } /* -------------------------------------------------------------------------- */ -template -inline void MaterialDamageIterative::packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const { +template +inline void MaterialDamageIterative::packData( + CommunicationBuffer & buffer, const Array & elements, + const SynchronizationTag & tag) const { if (tag == _gst_user_2) { - this->packElementDataHelper(this->damage, buffer, elements); + DataAccessor::packElementalDataHelper( + this->damage, buffer, elements, true, this->damage.getFEEngine()); } - return MaterialDamage::packElementData(buffer, elements, tag); + return MaterialDamage::packData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ -template -inline void MaterialDamageIterative::unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) { +template +inline void MaterialDamageIterative::unpackData( + CommunicationBuffer & buffer, const Array & elements, + const SynchronizationTag & tag) { if (tag == _gst_user_2) { - this->unpackElementDataHelper(this->damage, buffer, elements); + DataAccessor::unpackElementalDataHelper( + this->damage, buffer, elements, true, this->damage.getFEEngine()); } - - return MaterialDamage::unpackElementData(buffer, elements, tag); + return MaterialDamage::unpackData(buffer, elements, tag); } +} // akantu + /* -------------------------------------------------------------------------- */ 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 63d0ba181..551e85f48 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,75 @@ /** * @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" namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialDamageLinear::MaterialDamageLinear(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), - MaterialDamage(model, id), - K("K", *this) { +template +MaterialDamageLinear::MaterialDamageLinear( + SolidMechanicsModel & model, const ID & 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->registerParam("Gc", Gc, 2., _pat_parsable, "Gc"); this->K.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void MaterialDamageLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); Epsmin = Sigc / this->E; - Epsmax = 2 * Gc/ Sigc + Epsmin; + Epsmax = 2 * Gc / Sigc + Epsmin; this->K.setDefaultValue(Epsmin); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialDamageLinear::computeStress(ElementType el_type, - GhostType ghost_type) { +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); +INSTANTIATE_MATERIAL(damage_linear, MaterialDamageLinear); } // namespace akantu 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 d73dc78b1..35fdd0c96 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,211 @@ /** * @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 "communicator.hh" #include "solid_mechanics_model_RVE.hh" -#include +/* -------------------------------------------------------------------------- */ 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.) { +template +MaterialIterativeStiffnessReduction:: + MaterialIterativeStiffnessReduction(SolidMechanicsModel & model, + const ID & 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->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 +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); + for (auto ghost_type : ghost_types) { /// loop over all types in the filter - for(; it != last_type; ++it) { - ElementType el_type = *it; + for (auto & el_type : + this->element_filter.elementTypes(_ghost_type = ghost_type)) { /// get the stiffness on each quad point - Array::const_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); + auto 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(); + auto D_it = this->D(el_type, ghost_type).begin(); + auto D_end = this->D(el_type, ghost_type).end(); + /// get the ultimate strain on each quad + auto eps_u_it = this->eps_u(el_type, ghost_type).begin(); // compute the tangent and the ultimate strain for each quad for (; D_it != D_end; ++Sc_it, ++D_it, ++eps_u_it) { - *eps_u_it = ( (2. * this->Gf) / (*Sc_it * this->crack_band_width) ); - *D_it = *(Sc_it) / ((*eps_u_it) - ( (*Sc_it) / this->E ) ); + *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) { +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); + auto Sc_it = this->Sc(el_type, ghost_type).begin(); + auto dam_it = this->damage(el_type, ghost_type).begin(); + auto equivalent_stress_it = + this->equivalent_stress(el_type, ghost_type).begin(); + auto grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); + auto 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) { + 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); + 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); + *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 +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;; + /// 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) { + 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(); + auto equivalent_stress_it = + this->equivalent_stress(el_type, ghost_type).begin(); + auto equivalent_stress_end = + this->equivalent_stress(el_type, ghost_type).end(); + auto dam_it = this->damage(el_type, ghost_type).begin(); + auto reduction_it = this->reduction_step(el_type, ghost_type).begin(); + auto eps_u_it = this->eps_u(el_type, ghost_type).begin(); + auto Sc_it = this->Sc(el_type, ghost_type).begin(); + auto 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; - } + 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); + auto rve_model = dynamic_cast(&this->model); if (rve_model == NULL) { - StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); - comm.allReduce(&nb_damaged_elements, 1, _so_sum); + const auto & comm = this->model.getMesh().getCommunicator(); + comm.allReduce(nb_damaged_elements, SynchronizerOperation::_sum); } AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ - - -INSTANTIATE_MATERIAL(MaterialIterativeStiffnessReduction); - +INSTANTIATE_MATERIAL(iterative_stiffness_reduction, + MaterialIterativeStiffnessReduction); } // namespace akantu 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 63e3d00f5..21b7c819a 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,379 @@ /** * @file material_damage_iterative.cc * * @author Aurelia Isabel Cuba Ramos * * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage_iterative.hh" +#include "communicator.hh" #include "solid_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template +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) { + MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, + const ID & id) + : MaterialOrthotropicDamage(model, id), Sc("Sc", *this), + equivalent_stress("equivalent_stress", *this), + stress_dir("equiv_stress_dir", *this), norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); - this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold"); - this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "increase of damage in every step" ); - this->registerParam("dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1" ); - + this->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->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) { +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::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::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); + 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) + /// 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); + 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 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) { + 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.); + 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)) ) { + /// 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(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); + 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); - + 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) { +template +void MaterialOrthotropicDamageIterative::computeAllStresses( + GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress - if(ghost_type==_not_ghost) + 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); + this->model.getMesh().getCommunicator().allReduce( + norm_max_equivalent_stress, SynchronizerOperation::_max); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialOrthotropicDamageIterative::findMaxNormalizedEquivalentStress(ElementType el_type, - GhostType ghost_type) { +template +void MaterialOrthotropicDamageIterative:: + findMaxNormalizedEquivalentStress(ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); - if(ghost_type==_not_ghost) { + 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(); + auto equivalent_stress_it = e_stress.begin(); + auto 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); - + auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); + auto damage_directions_it = + this->damage_dir_vecs(el_type, _not_ghost) + .begin(this->spatial_dimension, this->spatial_dimension); + auto stress_dir_it = this->stress_dir(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); /// initialize the matrices for damage rotation results Matrix tmp(spatial_dimension, spatial_dimension); Matrix dam_in_computation_frame(spatial_dimension, spatial_dimension); Matrix dam_in_stress_frame(spatial_dimension, spatial_dimension); - for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it, ++damage_directions_it, ++stress_dir_it ) { - /// check if max equivalent stress for this element type is greater than the current norm_max_eq_stress + 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(); - } - } + (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) + } /// 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) { +template +void MaterialOrthotropicDamageIterative::computeStress( + ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); + MaterialOrthotropicDamage::computeStress(el_type, + ghost_type); - 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); + auto damage_iterator = + this->damage(el_type, ghost_type) + .begin(this->spatial_dimension, this->spatial_dimension); + auto damage_dir_it = + this->damage_dir_vecs(el_type, ghost_type) + .begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); - Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); - Matrix one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); - Matrix sqrt_one_minus_D_rotated(this->spatial_dimension, this->spatial_dimension); + Matrix 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)) ) { + /// 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(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); + 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); + 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); + 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 +template UInt MaterialOrthotropicDamageIterative::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., - "Your prescribed damage must be greater than zero"); + "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 + /// 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); + auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); + auto damage_directions_it = + this->damage_dir_vecs(el_type, _not_ghost) + .begin(this->spatial_dimension, this->spatial_dimension); + auto stress_dir_it = this->stress_dir(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); /// initialize the matrices for damage rotation results Matrix tmp(spatial_dimension, spatial_dimension); Matrix dam_in_computation_frame(spatial_dimension, spatial_dimension); Matrix dam_in_stress_frame(spatial_dimension, spatial_dimension); /// references to damage and directions of highest Gauss point Matrix q_dam = dam_it[q_global_num]; Matrix q_dam_dir = damage_directions_it[q_global_num]; Matrix q_stress_dir = stress_dir_it[q_global_num]; /// increment damage /// find the damage increment on this Gauss point /// rotate damage into stress frame - this->rotateIntoComputationFrame(q_dam, dam_in_computation_frame, q_dam_dir, tmp); - this->rotateIntoNewFrame(dam_in_computation_frame, dam_in_stress_frame, q_stress_dir, tmp); + 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; + 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); + this->model.getMesh().getCommunicator().allReduce( + nb_damaged_elements, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ -template -void MaterialOrthotropicDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { - MaterialOrthotropicDamage::updateEnergies(el_type, ghost_type); +template +void MaterialOrthotropicDamageIterative< + spatial_dimension>::updateEnergiesAfterDamage(ElementType el_type, + GhostType ghost_type) { + MaterialOrthotropicDamage::updateEnergies(el_type, + ghost_type); } /* -------------------------------------------------------------------------- */ - -INSTANTIATE_MATERIAL(MaterialOrthotropicDamageIterative); - +INSTANTIATE_MATERIAL(orthotropic_damage_iterative, + MaterialOrthotropicDamageIterative); } // 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 56668d59e..55837e443 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,129 @@ /** * @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" +#include "material_orthotropic_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ namespace akantu { /** * Material damage iterative * * parameters in the material files : - * - Sc + * - Sc */ -template -class MaterialOrthotropicDamageIterative : public MaterialOrthotropicDamage { +template +class MaterialOrthotropicDamageIterative + : public MaterialOrthotropicDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: + MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, + const ID & id = ""); - MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, const ID & id = ""); - - virtual ~MaterialOrthotropicDamageIterative() {}; + 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 updateEnergiesAfterDamage(ElementType el_type, + GhostType ghost_typ); - 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); + 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); + /// 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); + 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; + /// 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" } // namespace akantu #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_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 a03c0efee..91a98795a 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,322 @@ /** * @file material_orthotropic_damage_tmpl.hh * @author Aurelia Isabel Cuba Ramos * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the material class for the orthotropic * damage material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -template class Parent> -MaterialOrthotropicDamage::MaterialOrthotropicDamage(SolidMechanicsModel & model, - const ID & id) : - 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) { +template class Parent> +MaterialOrthotropicDamage::MaterialOrthotropicDamage( + SolidMechanicsModel & model, const ID & id) + : Parent(model, id), damage("damage", *this), + dissipated_energy("damage dissipated energy", *this), + int_sigma("integral of sigma", *this), + damage_dir_vecs("damage_principal_directions", *this) { AKANTU_DEBUG_IN(); - this->registerParam("eta" , eta , 2. , _pat_parsable | _pat_modifiable, "Damage sensitivity parameter" ); - this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" ); - + this->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->damage.initialize(spatial_dimension * spatial_dimension); this->dissipated_energy.initialize(1); - this->int_sigma .initialize(1); + this->int_sigma.initialize(1); this->damage_dir_vecs.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template class Parent> +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$ + * @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) { +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); + auto epsilon_p = this->gradu.previous(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); + auto 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(); + auto 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(); + auto ints = this->int_sigma(el_type, ghost_type).begin(); + auto 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) { +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); + 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); + auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); - /// get the directions of damage for the current element type + /// 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); + auto damage_directions_it = + dam_dirs.begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); - Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); + Matrix 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); + 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 +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)); + tangent *= (1 - dam(0, 0)); if (spatial_dimension == 2) { - Real min_val = std::min( (this->eta / spatial_dimension * trace_D), - this->max_damage ); + 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(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(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(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(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.)); + 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.)); -// /// 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)); + /// second row + tangent(1, 0) = + (C(0, 0) * S(0, 1) * S(0, 1) + C(1, 0) * S(1, 1) * S(1, 1) - + (min_val / 2. - 1. / 2) * (C(0, 0) + C(1, 0)) + + (O(1, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.)); + tangent(1, 1) = + (C(0, 1) * S(0, 1) * S(0, 1) + C(1, 1) * S(1, 1) * S(1, 1) - + (min_val / 2. - 1. / 2) * (C(0, 1) + C(1, 1)) + + (O(1, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.)); + tangent(1, 2) = (2. * C(2, 2) * S(0, 1) * S(1, 1) + + (2. * C(2, 2) * O(0, 1) * O(1, 1)) / (trace_D - 2.)); - - - - - - + /// third row + tangent(2, 0) = + ((O(0, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.) + + C(0, 0) * S(0, 0) * S(0, 1) + C(1, 0) * S(0, 1) * S(1, 1)); + + tangent(2, 1) = + ((O(0, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.) + + C(0, 1) * S(0, 0) * S(0, 1) + C(1, 1) * S(0, 1) * S(1, 1)); + tangent(2, 2) = ((2. * C(2, 2) * O(0, 1) * O(0, 1)) / (trace_D - 2.) + + C(2, 2) * S(0, 1) * S(0, 1) + C(2, 2) * S(0, 0) * S(1, 1)); } if (spatial_dimension == 3) { - - } - - } /* -------------------------------------------------------------------------- */ -template class Parent> -inline void -MaterialOrthotropicDamage::computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & sqrt_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term) { - - // /// 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.; - +template class Parent> +inline void MaterialOrthotropicDamage:: + computeDamageAndStressOnQuad(Matrix & sigma, + Matrix & one_minus_D, + Matrix & sqrt_one_minus_D, + Matrix & damage, + Matrix & first_term, + Matrix & third_term) { /// Definition of Cauchy stress based on second order damage tensor: /// "Anisotropic damage modelling of biaxial behaviour and rupture /// of concrete strucutres", Ragueneau et al., 2008, Eq. 7 first_term.mul(sqrt_one_minus_D, sigma); first_term *= sqrt_one_minus_D; Real second_term = 0; for (UInt i = 0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) - second_term += sigma(i,j) * one_minus_D(i,j); - } + second_term += 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())); + 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); - +template class Parent> +inline void MaterialOrthotropicDamage:: + rotateIntoComputationFrame(const Matrix & to_rotate, + Matrix & rotated, + const Matrix & damage_directions, + Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ -template class Parent> -inline void MaterialOrthotropicDamage -::rotateIntoNewFrame(const Matrix & to_rotate, - Matrix & rotated, - const Matrix & damage_directions, - Matrix & rotation_tmp) { - - /// copy matrix - // rotated.copy(to_rotate); - // to_rotate.mul(rotated, damage_directions); - // rotated.mul(damage_directions, to_rotate); - +template class Parent> +inline void +MaterialOrthotropicDamage::rotateIntoNewFrame( + const Matrix & to_rotate, Matrix & rotated, + const Matrix & damage_directions, Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ -template class Parent> -inline void MaterialOrthotropicDamage -::computeOneMinusD(Matrix & one_minus_D, const Matrix & damage) { +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> +template class Parent> inline void -MaterialOrthotropicDamage::computeSqrtOneMinusD(const Matrix & one_minus_D, Matrix & sqrt_one_minus_D) { +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 +/// 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"); + 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 - +#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)); + 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(); - -// /// 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, -// _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); -// } - -// /* -------------------------------------------------------------------------- */ - } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.cc index cc4a6e5c6..48b476fb9 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,108 @@ /** * @file material_viscoplastic.cc * * @author Ramin Aghababaei * * - * @brief Specialization of the material class for isotropic viscoplastic (small deformation) + * @brief Specialization of the material class for isotropic viscoplastic + * (small deformation) * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ - /* -------------------------------------------------------------------------- */ #include "material_viscoplastic.hh" #include "solid_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialViscoPlastic::MaterialViscoPlastic(SolidMechanicsModel & model, const ID & id) : - Material(model, id), - MaterialPlastic(model, id) { +template +MaterialViscoPlastic::MaterialViscoPlastic(SolidMechanicsModel & model, + const ID & id) + : MaterialPlastic(model, id) { AKANTU_DEBUG_IN(); - this->registerParam("rate", rate, 0., _pat_parsable | _pat_modifiable, "Rate sensitivity component"); - this->registerParam("edot0", edot0, 0., _pat_parsable | _pat_modifiable, "Reference strain rate"); - this->registerParam("ts", ts, 0., _pat_parsable | _pat_modifiable, "Time Step"); + 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(); +template +void MaterialViscoPlastic::computeStress(ElementType el_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); - Real * iso_hardening = this->iso_hardening(el_type, ghost_type).storage(); + 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); + auto previous_grad_u_it = + this->gradu.previous(el_type, ghost_type).begin(dim, dim); - Array::matrix_iterator previous_sigma_it = - this->stress.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); + auto previous_sigma_it = + this->stress.previous(el_type, ghost_type).begin(dim, dim); - Array::matrix_iterator inelastic_strain_it = - this->inelastic_strain(el_type, ghost_type).begin(spatial_dimension,spatial_dimension); + auto inelastic_strain_it = + this->inelastic_strain(el_type, ghost_type).begin(dim, dim); - Array::matrix_iterator previous_inelastic_strain_it = - this->inelastic_strain.previous(el_type, ghost_type).begin(spatial_dimension,spatial_dimension); + auto previous_inelastic_strain_it = + this->inelastic_strain.previous(el_type, ghost_type).begin(dim, dim); - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + 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); + 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; + ++inelastic_strain_it; + ++iso_hardening; + ++previous_grad_u_it; - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - AKANTU_DEBUG_OUT(); + 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); +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_strain_it = - this->gradu.previous(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); + auto previous_sigma_it = + this->stress.previous(el_type, ghost_type).begin(dim, dim); - Real * iso_hardening= this->iso_hardening(el_type, ghost_type).storage(); + auto previous_strain_it = + this->gradu.previous(el_type, ghost_type).begin(dim, dim); + Real * iso_hardening = this->iso_hardening(el_type, ghost_type).storage(); - MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); + MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); - Matrix & previous_grad_u = *previous_strain_it; - Matrix & previous_sigma_tensor = *previous_sigma_it; + 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; + 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(); + AKANTU_DEBUG_OUT(); } -INSTANTIATE_MATERIAL(MaterialViscoPlastic); +INSTANTIATE_MATERIAL(visco_plastic, MaterialViscoPlastic); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic_inline_impl.cc b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic_inline_impl.cc index 1c790f454..fcb55190a 100644 --- a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic_inline_impl.cc +++ b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic_inline_impl.cc @@ -1,92 +1,92 @@ /** * @file material_viscoplastic_inline_impl.cc * * @author Ramin Aghababaei * * * @brief Implementation of the inline functions of the material viscoplastic * * @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 #include "material_viscoplastic.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void MaterialViscoPlastic::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 { // Infinitesimal plastic // Compute stress magnitude Real s = sigma.doubleDot(sigma); Real sigma_mag=sqrt(s); // Compute plastic strain increment Real factor = (this->ts * this->edot0 * pow(sigma_mag, (this->rate - 1.)) / pow(this->sigma_y + iso_hardening, this->rate)); Matrix delta_inelastic_strain(sigma); delta_inelastic_strain *= factor; // Compute plastic strain increment magnitude s = delta_inelastic_strain.doubleDot(delta_inelastic_strain); Real dep_mag = std::sqrt(s); Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; // Update stress and plastic strain Matrix grad_u_elastic(dim, dim); grad_u_elastic = grad_delta_u; grad_u_elastic -= delta_inelastic_strain; Matrix sigma_elastic(dim, dim); MaterialElastic::computeStressOnQuad(grad_u_elastic, sigma_elastic); sigma += sigma_elastic; inelastic_strain += delta_inelastic_strain; //Update resistance stress iso_hardening = iso_hardening + this->h * dep_mag; MaterialPlastic::computeStressAndInelasticStrainOnQuad(grad_delta_u, sigma, previous_sigma, inelastic_strain, previous_inelastic_strain, delta_inelastic_strain); } /* -------------------------------------------------------------------------- */ template inline void MaterialViscoPlastic::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 { + const Matrix & /*grad_u*/, + const Matrix & /*previous_grad_u*/, + const Matrix & /*sigma_tensor*/, + const Matrix & /*previous_sigma_tensor*/, + const Real & /*iso_hardening*/) const { UInt cols = tangent.cols(); UInt rows = tangent.rows(); for (UInt m = 0; m < rows; ++m) { UInt i = VoigtHelper::vec[m][0]; UInt j = VoigtHelper::vec[m][1]; for (UInt n = 0; n < cols; ++n) { UInt k = VoigtHelper::vec[n][0]; UInt l = VoigtHelper::vec[n][1]; tangent(m,n) = (i==k) * (j==l) * 2. * this->mu + (i==j) * (k==l) * this->lambda; tangent(m,n) -= (m==n) * (m>=dim) * this->mu; } } } 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 d28ca2c11..797fbb275 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,125 @@ /** * @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" 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) { +template +MaterialStiffnessProportional::MaterialStiffnessProportional( + SolidMechanicsModel & model, const ID & 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->registerParam("Alpha", alpha, 0., _pat_parsable | _pat_modifiable, + "Artificial viscous ratio"); this->stress_viscosity.initialize(spatial_dimension * spatial_dimension); - this->stress_elastic .initialize(spatial_dimension * spatial_dimension); - + this->stress_elastic.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void MaterialStiffnessProportional::initMaterial() { AKANTU_DEBUG_IN(); MaterialElastic::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialStiffnessProportional::computeStress(ElementType el_type, - GhostType ghost_type) { +template +void MaterialStiffnessProportional::computeStress( + ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array & elem_filter = this->element_filter (el_type, ghost_type); + 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); + Array & stress_el = stress_elastic(el_type, ghost_type); MaterialElastic::computeStress(el_type, ghost_type); Array & velocity = this->model.getVelocity(); Array strain_rate(0, spatial_dimension * spatial_dimension, - "strain_rate"); + "strain_rate"); - this->model.getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate, - spatial_dimension, - el_type, ghost_type, elem_filter); + 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); + strain_rate.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_visc_it = - stress_visc.begin(spatial_dimension, spatial_dimension); + stress_visc.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_el_it = - stress_el.begin(spatial_dimension, spatial_dimension); + stress_el.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); - Matrix & grad_v = *strain_rate_it; + Matrix & grad_v = *strain_rate_it; Matrix & sigma_visc = *stress_visc_it; - Matrix & sigma_el = *stress_el_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) { +template +void MaterialStiffnessProportional::computePotentialEnergy( + ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - if(ghost_type != _not_ghost) return; + 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); + 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); + 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); - +INSTANTIATE_MATERIAL(ve_stiffness_prop, MaterialStiffnessProportional); } // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index 7c65113e6..bc1b06b62 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,340 +1,340 @@ /** * @file heat_transfer_model.hh * * @author Guillaume Anciaux * @author Lucas Frerot * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Dec 08 2015 * * @brief Model of Heat Transfer * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_HEAT_TRANSFER_MODEL_HH__ #define __AKANTU_HEAT_TRANSFER_MODEL_HH__ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; } namespace akantu { class HeatTransferModel : public Model, public DataAccessor, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; void predictor() override; /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// compute the internal heat flux void assembleInternalHeatRate(); /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// assemble the conductivity matrix void assembleConductivityMatrix(); /// assemble the conductivity matrix void assembleCapacity(); /// assemble the conductivity matrix void assembleCapacity(GhostType ghost_type); /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); private: /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(const GhostType & ghost_type); /// assemble the conductivity matrix (w/ ghost type) template void assembleConductivityMatrix(const GhostType & ghost_type); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point void computeKgradT(const GhostType & ghost_type); /// compute the thermal energy Real computeThermalEnergyByNode(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; inline UInt getNbData(const Array & indexes, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); - virtual void dump(); + void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Density, density, Real); AKANTU_GET_MACRO(Capacity, capacity, Real); /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux AKANTU_GET_MACRO(InternalHeatRate, *internal_heat_rate, Array &); /// get the lumped capacity AKANTU_GET_MACRO(CapacityLumped, *capacity_lumped, Array &); /// get the boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the external heat rate vector AKANTU_GET_MACRO(ExternalHeatRate, *external_heat_rate, Array &); /// get the temperature gradient AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); /// internal variables AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KgradT, k_gradt_on_qpoints, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Array &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Array &); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id); /// get the thermal energy for a given element Real getThermalEnergy(const ElementType & type, UInt index); /// get the thermal energy for a given element Real getThermalEnergy(); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ----------------------------------------------------------------------- */ template void getThermalEnergy(iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const; template void allocNodalField(Array *& array, const ID & name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// number of iterations UInt n_iter; /// time step Real time_step; /// temperatures array Array * temperature{nullptr}; /// temperatures derivatives array Array * temperature_rate{nullptr}; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Array * increment{nullptr}; /// the density Real density; /// the speed of the changing temperature ElementTypeMapArray temperature_gradient; /// temperature field on quadrature points ElementTypeMapArray temperature_on_qpoints; /// conductivity tensor on quadrature points ElementTypeMapArray conductivity_on_qpoints; /// vector k \grad T on quad points ElementTypeMapArray k_gradt_on_qpoints; /// external flux vector Array * external_heat_rate{nullptr}; /// residuals array Array * internal_heat_rate{nullptr}; // lumped vector Array * capacity_lumped{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; // realtime Real time; /// capacity Real capacity; // conductivity matrix Matrix conductivity; // linear variation of the conductivity (for temperature dependent // conductivity) Real conductivity_variation; // reference temperature for the interpretation of temperature variation Real T_ref; // the biggest parameter of conductivity matrix Real conductivitymax; UInt temperature_release{0}; UInt conductivity_matrix_release{0}; std::unordered_map initial_conductivity{{_not_ghost, true}, {_ghost, true}}; std::unordered_map conductivity_release{{_not_ghost, 0}, {_ghost, 0}}; }; } // akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model_inline_impl.cc" #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */