diff --git a/src/material_FE2/material_FE2.cc b/src/material_FE2/material_FE2.cc index 42d84c365..84faf9648 100644 --- a/src/material_FE2/material_FE2.cc +++ b/src/material_FE2/material_FE2.cc @@ -1,189 +1,199 @@ /** * @file material_FE2.cc * * @author Aurelia Isabel Cuba Ramos * * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_FE2.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialFE2::MaterialFE2(SolidMechanicsModel & model, const ID & id) : - Material(model, id), + Material(model, id), Parent(model, id), C("material_stiffness", *this) { AKANTU_DEBUG_IN(); this->C.initialize(voigt_h::size * voigt_h::size); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialFE2::~MaterialFE2() { AKANTU_DEBUG_IN(); for (UInt i = 0; i < RVEs.size(); ++i) { delete meshes[i]; delete RVEs[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initialize() { this->registerParam("element_type" ,el_type, _triangle_3 , _pat_parsable | _pat_modifiable, "element type in RVE mesh" ); this->registerParam("mesh_file" ,mesh_file , _pat_parsable | _pat_modifiable, "the mesh file for the RVE"); this->registerParam("nb_gel_pockets" ,nb_gel_pockets , _pat_parsable | _pat_modifiable, "the number of gel pockets in each RVE"); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initMaterial() { AKANTU_DEBUG_IN(); - Material::initMaterial(); + Parent::initMaterial(); /// compute the number of integration points in this material and resize the RVE vector UInt nb_integration_points = this->element_filter(this->el_type, _not_ghost).getSize() * this->fem->getNbIntegrationPoints(this->el_type); RVEs.resize(nb_integration_points); meshes.resize(nb_integration_points); /// create a Mesh and SolidMechanicsModel on each integration point of the material std::vector::iterator RVE_it = RVEs.begin(); std::vector::iterator mesh_it = meshes.begin(); StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); UInt prank = comm.whoAmI(); Array::matrix_iterator C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size); for (UInt i = 1; i < nb_integration_points+1; ++RVE_it, ++mesh_it, ++i, ++C_it) { std::stringstream mesh_name; mesh_name << "RVE_mesh_" << prank; std::stringstream rve_name; rve_name << "SMM_RVE_" << prank; *mesh_it = new Mesh(spatial_dimension, mesh_name.str(), i); (*mesh_it)->read(mesh_file); *RVE_it = new SolidMechanicsModelRVE(*(*(mesh_it)), true, this->nb_gel_pockets, _all_dimensions, rve_name.str(), i); (*RVE_it)->initFull(); /// compute intial stiffness of the RVE (*RVE_it)->homogenizeStiffness(*C_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); + + // Compute thermal stresses first + + Parent::computeStress(el_type, ghost_type); + Array::const_scalar_iterator sigma_th_it = + this->sigma_th(el_type, ghost_type).begin(); + // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation - Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, + Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); // create vectors to store stress and strain in Voigt notation // for efficient computation of stress Vector voigt_strain(voigt_h::size); Vector voigt_stress(voigt_h::size); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); const Matrix & C_mat = *C_it; + const Real & sigma_th = *sigma_th_it; /// copy strains in Voigt notation for(UInt I = 0; I < voigt_h::size; ++I) { /// copy stress in Real voigt_factor = voigt_h::factors[I]; UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; } // compute stresses in Voigt notation voigt_stress.mul(C_mat, voigt_strain); /// copy stresses back in full vectorised notation for(UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; - sigma(i, j) = sigma(j, i) = voigt_stress(I); + sigma(i, j) = sigma(j, i) = voigt_stress(I)+ (i == j) * sigma_th; } - + ++C_it; - + ++sigma_th_it; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); tangent.copy(*C_it); ++C_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); std::vector::iterator RVE_it = RVEs.begin(); std::vector::iterator RVE_end = RVEs.end(); Array::matrix_iterator C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size); Array::matrix_iterator gradu_it = this->gradu(this->el_type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator eigen_gradu_it = this->eigengradu(this->el_type).begin(spatial_dimension, spatial_dimension); for (; RVE_it != RVE_end; ++RVE_it, ++C_it, ++gradu_it, ++eigen_gradu_it) { /// apply boundary conditions based on the current macroscopic displ. gradient (*RVE_it)->applyBoundaryConditions(*gradu_it); /// advance the ASR in every RVE (*RVE_it)->advanceASR(prestrain); /// compute the average eigen_grad_u (*RVE_it)->homogenizeEigenGradU(*eigen_gradu_it); /// compute the new effective stiffness of the RVE (*RVE_it)->homogenizeStiffness(*C_it); } AKANTU_DEBUG_OUT(); } + INSTANTIATE_MATERIAL(MaterialFE2); __END_AKANTU__ diff --git a/src/material_FE2/material_FE2.hh b/src/material_FE2/material_FE2.hh index b39d13d1c..02f0843d0 100644 --- a/src/material_FE2/material_FE2.hh +++ b/src/material_FE2/material_FE2.hh @@ -1,111 +1,123 @@ /** * @file material_FE2.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" +#include "material_thermal.hh" #include "solid_mechanics_model_RVE.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_FE_2_HH__ #define __AKANTU_MATERIAL_FE_2_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// /!\ This material works ONLY for meshes with a single element type!!!!! /* -------------------------------------------------------------------------- */ /** * MaterialFE2 * * parameters in the material files : * - mesh_file */ +// Emil - 27.01.2018 - re-inheriting from MaterialThermal and PlaneStressToolBox + +//template +//class MaterialFE2 : public virtual Material { +// /* ------------------------------------------------------------------------ */ +// /* Constructors/Destructors */ +// /* ------------------------------------------------------------------------ */ template -class MaterialFE2 : public virtual Material { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: +class MaterialFE2 : public MaterialThermal { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +private: + typedef MaterialThermal Parent; + +// Emil - 27.01.2018 + public: MaterialFE2(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialFE2(); typedef VoigtHelper voigt_h; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); /// advance alkali-silica reaction void advanceASR(const Matrix & prestrain); private: void initialize(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying RVE at each integration point std::vector RVEs; /// Meshes for all RVEs std::vector meshes; /// the element type of the associated mesh (this material handles only one type!!) ElementType el_type; /// the name of RVE mesh file ID mesh_file; /// Elastic stiffness tensor at each Gauss point (in voigt notation) InternalField C; /// number of gel pockets in each underlying RVE UInt nb_gel_pockets; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_FE2_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_FE_2_HH__ */ diff --git a/src/material_damage/material_damage_iterative.cc b/src/material_damage/material_damage_iterative.cc index 6362837f1..039847eca 100644 --- a/src/material_damage/material_damage_iterative.cc +++ b/src/material_damage/material_damage_iterative.cc @@ -1,234 +1,234 @@ /** * @file material_damage_iterative.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative.hh" #include "solid_mechanics_model_RVE.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialDamageIterative::MaterialDamageIterative(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), Sc("Sc", *this), reduction_step("damage_step", *this), equivalent_stress("equivalent_stress", *this), max_reductions(0), norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold"); this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "prescribed damage" ); this->registerParam("dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1" ); this->registerParam("dam_tolerance", dam_tolerance, 0.01, _pat_parsable | _pat_modifiable, "damage tolerance to decide if quadrature point will be damageed" ); - this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" ); - this->registerParam("max_reductions", max_reductions, UInt(10), _pat_parsable | _pat_modifiable, "max reductions"); + this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" ); + this->registerParam("max_reductions", max_reductions, UInt(10), _pat_parsable | _pat_modifiable, "max reductions"); this->use_previous_stress = true; this->use_previous_gradu = true; this->Sc.initialize(1); this->equivalent_stress.initialize(1); this->reduction_step.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); Array::const_iterator Sc_it = Sc(el_type, ghost_type).begin(); Array::iterator equivalent_stress_it = equivalent_stress(el_type, ghost_type).begin(); Array::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); Real * dam = this->damage(el_type, ghost_type).storage(); Matrix sigma(spatial_dimension, spatial_dimension); for(;grad_u_it != grad_u_end; ++ grad_u_it) { sigma.clear(); MaterialElastic::computeStressOnQuad(*grad_u_it, sigma, 0.); computeDamageAndStressOnQuad(sigma,*dam); /// compute eigenvalues sigma.eig(eigenvalues); /// find max eigenvalue and normalize by tensile strength *equivalent_stress_it = *(std::max_element(eigenvalues.storage(), eigenvalues.storage() + spatial_dimension)) / *(Sc_it); ++Sc_it; ++equivalent_stress_it; ++dam; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress if(ghost_type==_not_ghost) norm_max_equivalent_stress = 0; MaterialDamage::computeAllStresses(ghost_type); /// find global Gauss point with highest stress const SolidMechanicsModelRVE * rve_model = dynamic_cast(this->model); if (rve_model == NULL) { /// is no RVE model StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&norm_max_equivalent_stress, 1, _so_max); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type==_not_ghost) { // const Array & e_stress = equivalent_stress(el_type); // if (e_stress.begin() != e_stress.end() ) { // Array::const_iterator equivalent_stress_it_max = std::max_element(e_stress.begin(),e_stress.end()); // /// check if max equivalent stress for this element type is greater than the current norm_max_eq_stress // if (*equivalent_stress_it_max > norm_max_equivalent_stress) // norm_max_equivalent_stress = *equivalent_stress_it_max; // } const Array & e_stress = equivalent_stress(el_type); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator dam_it = dam.begin(); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it ) { /// check if max equivalent stress for this element type is greater than the current norm_max_eq_stress and if the element is not already fully damaged if (*equivalent_stress_it > norm_max_equivalent_stress && *dam_it < max_damage) { norm_max_equivalent_stress = *equivalent_stress_it; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialDamage::computeStress(el_type, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeDamageAndStressOnQuad(sigma,*dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type); norm_max_equivalent_stress = 0; findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialDamageIterative::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., "Your prescribed damage must be greater than zero"); if (norm_max_equivalent_stress >= 1.) { AKANTU_DEBUG_IN(); /// update the damage only on non-ghosts elements! Doesn't make sense to update on ghost. GhostType ghost_type = _not_ghost;; Mesh::type_iterator it = this->model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { ElementType el_type = *it; const Array & e_stress = equivalent_stress(el_type); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator dam_it = dam.begin(); Array::scalar_iterator reduction_it = this->reduction_step(el_type, ghost_type).begin(); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it,++reduction_it ) { /// check if damage occurs if (*equivalent_stress_it >= (1-dam_tolerance)*norm_max_equivalent_stress) { /// check if this element can still be damaged if (*reduction_it == this->max_reductions) continue; *reduction_it += 1; if (*reduction_it == this->max_reductions) { *dam_it = max_damage; } else { *dam_it += prescribed_dam; } nb_damaged_elements += 1; } } } } const SolidMechanicsModelRVE * rve_model = dynamic_cast(this->model); if (rve_model == NULL) { StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); } AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialDamageIterative::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) { MaterialDamage::updateEnergies(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialDamageIterative); __END_AKANTU__