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 2dcdbaab9..59cedd99d 100644 --- a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc +++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc @@ -1,196 +1,202 @@ /** * @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 MaterialFE2::MaterialFE2(SolidMechanicsModel & model, 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() = 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::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// create a Mesh and SolidMechanicsModel on each integration point of the /// material const auto & comm = this->model.getMesh().getCommunicator(); UInt prank = comm.whoAmI(); 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)); + spatial_dimension, "RVE_mesh_" + std::to_string(q), 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)); + "SMM_RVE_" + std::to_string(q), q + 1)); auto & RVE = *RVEs.back(); RVE.initFull(_analysis_method = _static); /// compute intial stiffness of the RVE RVE.homogenizeStiffness(C); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Compute thermal stresses first Parent::computeStress(el_type, ghost_type); Array::const_scalar_iterator sigma_th_it = this->sigma_th(el_type, ghost_type).begin(); // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); // create vectors to store stress and strain in Voigt notation // for efficient computation of stress Vector voigt_strain(voigt_h::size); Vector voigt_stress(voigt_h::size); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); const Matrix & C_mat = *C_it; const Real & sigma_th = *sigma_th_it; /// copy strains in Voigt notation for (UInt I = 0; I < voigt_h::size; ++I) { /// copy stress in Real voigt_factor = voigt_h::factors[I]; UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; } // compute stresses in Voigt notation voigt_stress.mul(C_mat, voigt_strain); /// copy stresses back in full vectorised notation for (UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; sigma(i, j) = sigma(j, i) = voigt_stress(I) + (i == j) * sigma_th; } ++C_it; ++sigma_th_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); tangent.copy(*C_it); ++C_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::advanceASR( const Matrix & prestrain) { AKANTU_DEBUG_IN(); 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))) { + make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), + this->delta_T(this->el_type))) { auto & RVE = *(std::get<0>(data)); /// apply boundary conditions based on the current macroscopic displ. /// gradient RVE.applyBoundaryConditions(std::get<1>(data)); + /// apply homogeneous temperature field to each RVE to obtain thermoelastic + /// effect + RVE.applyHomogeneousTemperature(std::get<4>(data)); + + /// advance the ASR in every RVE RVE.advanceASR(prestrain); /// compute the average eigen_grad_u RVE.homogenizeEigenGradU(std::get<2>(data)); /// compute the new effective stiffness of the RVE RVE.homogenizeStiffness(std::get<3>(data)); } AKANTU_DEBUG_OUT(); } 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 64d6c4fa4..0d24d5168 100644 --- a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh +++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh @@ -1,122 +1,112 @@ /** * @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 "material.hh" #include "material_thermal.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!!!!! /* -------------------------------------------------------------------------- */ /** * MaterialFE2 * * parameters in the material files : * - mesh_file */ -// Emil - 27.01.2018 - re-inheriting from MaterialThermal and PlaneStressToolBox - -// template -// class MaterialFE2 : public virtual Material { -// /* ------------------------------------------------------------------------ -// */ -// /* Constructors/Destructors */ -// /* ------------------------------------------------------------------------ -// */ template class MaterialFE2 : public MaterialThermal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: typedef MaterialThermal Parent; - // Emil - 27.01.2018 public: MaterialFE2(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialFE2(); typedef VoigtHelper voigt_h; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); /// advance alkali-silica reaction void advanceASR(const Matrix & prestrain); private: void initialize(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying RVE at each integration point std::vector> RVEs; /// Meshes for all RVEs std::vector> meshes; /// the element type of the associated mesh (this material handles only one /// type!!) ElementType el_type; /// the name of RVE mesh file ID mesh_file; /// Elastic stiffness tensor at each Gauss point (in voigt notation) InternalField C; /// number of gel pockets in each underlying RVE UInt nb_gel_pockets; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_FE2_inline_impl.cc" } // 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 27306a2ed..45e1eedcf 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,553 +1,604 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "parser.hh" #include "sparse_matrix.hh" +#include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// find the four corner nodes of the RVE findCornerNodes(); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); // /// enforce periodicity on the displacement fluctuations // auto surface_pair_1 = std::make_pair("top", "bottom"); // auto surface_pair_2 = std::make_pair("right", "left"); // SurfacePairList surface_pairs_list; // surface_pairs_list.push_back(surface_pair_1); // surface_pairs_list.push_back(surface_pair_2); // TODO: To Nicolas correct the PBCs // this->setPBC(surface_pairs_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); - this->initMaterials(); - + //this->initMaterials(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); /// compute the volume of the RVE GhostType gt = _not_ghost; - for (auto element_type : this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { + for (auto element_type : + this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { Array Volume(this->mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem.integrate(Volume, element_type); } std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id; // << this->memory_id - 1; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); + this->addDumpField("delta_T"); this->dump(); this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will /// be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModelRVE::applyHomogeneousTemperature( + const Real & temperature) { + + for (UInt m = 0; m < this->getNbMaterials(); ++m) { + Material & mat = this->getMaterial(m); + + const ElementTypeMapArray & filter_map = mat.getElementFilter(); + + Mesh::type_iterator type_it = filter_map.firstType(spatial_dimension); + Mesh::type_iterator type_end = filter_map.lastType(spatial_dimension); + // Loop over all element types + for (; type_it != type_end; ++type_it) { + const Array & filter = filter_map(*type_it); + if (filter.size() == 0) + continue; + + Array & delta_T = mat.getArray("delta_T", *type_it); + Array::scalar_iterator delta_T_it = delta_T.begin(); + Array::scalar_iterator delta_T_end = delta_T.end(); + + for (; delta_T_it != delta_T_end; ++delta_T_it) { + *delta_T_it = temperature; + } + } + } +} + /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { AKANTU_DEBUG_IN(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, spatial_dimension))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } for (UInt i = 0; i < corner_nodes.size(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_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 - for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { + for (auto element_type : + mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } /// advance the damage MaterialDamageIterative<2> & mat_paste = dynamic_cast &>(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast &>(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", _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()); std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); if (this->nb_dumps % 10 == 0) { this->dump(); } this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); Real average = 0; - for (auto element_type : mesh.elementTypes(_element_kind = _ek_not_defined)) { + GhostType gt = _not_ghost; + for (auto element_type : + mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & stress_vec = this->materials[m]->getStress(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) - average += int_stress_vec( - k, row_index * spatial_dimension + col_index); // 3 is the value - // for the yy (in - // 3D, the value is - // 4) + average += int_stress_vec(k, row_index * spatial_dimension + + col_index); // 3 is the value + // for the yy (in + // 3D, the value is + // 4) } } else if (field_type == "strain") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & gradu_vec = this->materials[m]->getGradU(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else { AKANTU_DEBUG_ERROR("Averaging not implemented for this field!!!"); } } return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// clear the eigenstrain Matrix zero_eigengradu(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = // 0); // this->fillCracks(saved_damage); /// virtual test 1: H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", _scc_solution); this->solveStep(); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU( Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top - bottom); Real eps = box_size * 1e-6; auto tmp = std::make_shared(*this, box_size, "gel", this->nb_gel_pockets, eps); tmp->setFallback(material_selector); material_selector = tmp; } - SolidMechanicsModel::initMaterials(); + this->assignMaterialToElements(); + // synchronize the element material arrays + this->synchronize(_gst_material_id); + + for (auto & material : materials) { + /// init internals properties + const auto type = material->getID(); + if (type.find("material_FE2") != std::string::npos) + continue; + material->initMaterial(); + } + + this->synchronize(_gst_smm_init_mat); + + if (this->non_local_manager) { + this->non_local_manager->initialize(); + } + //SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = this->getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; Real E = mat->get("E"); auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks( const ElementTypeMapReal & saved_damage) { for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); - auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); + 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 c38055a1a..b509066c6 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,230 +1,234 @@ /** * @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 "aka_grid_dynamic.hh" #include "solid_mechanics_model.hh" #include /* -------------------------------------------------------------------------- */ 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); virtual ~SolidMechanicsModelRVE(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: void initFullImpl(const ModelOptions & option) override; /// initialize the materials void initMaterials() override; public: /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditions(const Matrix & displacement_gradient); - /// advance the reactions -> grow gel and apply homogenized properties + /// apply homogeneous temperature field from the macroscale level to the RVEs + virtual void + applyHomogeneousTemperature(const Real & temperature); + +/// advance the reactions -> grow gel and apply homogenized properties void advanceASR(const Matrix & prestrain); /// compute average stress or strain in the model Real averageTensorField(UInt row_index, UInt col_index, const ID & field_type); /// compute effective stiffness of the RVE void homogenizeStiffness(Matrix & C_macro); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ inline void unpackData(CommunicationBuffer & buffer, const Array & index, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(Volume, volume, Real); private: /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); void fillCracks(ElementTypeMapReal & saved_damage); void drainCracks(const ElementTypeMapReal & saved_damage); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ /// volume of the RVE Real volume; /// corner nodes 1, 2, 3, 4 (see Leonardo's thesis, page 98) Array corner_nodes; /// bottom nodes std::unordered_set bottom_nodes; /// left nodes std::unordered_set left_nodes; /// standard mat selector or user one bool use_RVE_mat_selector; /// the number of gel pockets inside the RVE UInt nb_gel_pockets; /// dump counter UInt nb_dumps; }; inline void SolidMechanicsModelRVE::unpackData(CommunicationBuffer & buffer, const Array & index, const SynchronizationTag & tag) { SolidMechanicsModel::unpackData(buffer, index, tag); - if (tag == _gst_smm_uv) { - 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)]); - } - } - } +// if (tag == _gst_smm_uv) { +// 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 { public: GelMaterialSelector(SolidMechanicsModel & model, const Real box_size, const std::string & gel_material, const UInt nb_gel_pockets, Real /*tolerance*/ = 0.) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), box_size(box_size) { Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = model.getSpatialDimension(); 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; checked_baries.insert(bary_id); el.element = bary_id; if (MeshDataMaterialSelector::operator()(el) == 1) continue; /// element belongs to paste gel_pockets.push_back(el); placed_gel_pockets += 1; } } UInt operator()(const Element & elem) { UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index == 1) return temp_index; std::vector::const_iterator iit = gel_pockets.begin(); std::vector::const_iterator eit = gel_pockets.end(); if (std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material); ; } return 0; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; }; } // namespace akantu ///#include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ */