diff --git a/src/material_FE2/material_FE2.cc b/src/material_FE2/material_FE2.cc index c612daf42..f494b103f 100644 --- a/src/material_FE2/material_FE2.cc +++ b/src/material_FE2/material_FE2.cc @@ -1,182 +1,188 @@ /** * @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), 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(); /// 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(); // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation - AKANTU_DEBUG_IN(); 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; /// 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); } ++C_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 new effective stiffness of the RVE (*RVE_it)->homogenizeStiffness(*C_it); /// compute the average eigen_grad_u (*RVE_it)->homogenizeEigenGradU(*eigen_gradu_it); } + AKANTU_DEBUG_OUT(); } INSTANTIATE_MATERIAL(MaterialFE2); __END_AKANTU__ diff --git a/src/material_FE2/solid_mechanics_model_RVE.cc b/src/material_FE2/solid_mechanics_model_RVE.cc index db9f1dace..38cfc27d5 100644 --- a/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,492 +1,514 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" #include "material_damage_iterative.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif +#ifdef AKANTU_USE_PETSC +#include "solver_petsc.hh" +#include "petsc_matrix.hh" +#endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), static_communicator_dummy(StaticCommunicator::getStaticCommunicatorDummy()), - nb_gel_pockets(nb_gel_pockets) { - + nb_gel_pockets(nb_gel_pockets), + nb_dumps(0) { + AKANTU_DEBUG_IN(); /// create node groups for PBCs mesh.createGroupsFromMeshData("physical_names"); /// find the four corner nodes of the RVE findCornerNodes(); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const ElementGroup & bottom = mesh.getElementGroup("bottom"); bottom_nodes.insert( bottom.node_begin(), bottom.node_end() ); const ElementGroup & left = mesh.getElementGroup("left"); left_nodes.insert( left.node_begin(), left.node_end() ); /// enforce periodicity on the displacement fluctuations SurfacePair surface_pair_1 = std::make_pair("top","bottom"); SurfacePair surface_pair_2 = std::make_pair("right","left"); SurfacePairList surface_pairs_list; surface_pairs_list.push_back(surface_pair_1); surface_pairs_list.push_back(surface_pair_2); this->setPBC(surface_pairs_list); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() { delete static_communicator_dummy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFull(const ModelOptions & options) { + AKANTU_DEBUG_IN(); SolidMechanicsModel::initFull(options); this->initMaterials(); /// compute the volume of the RVE FEEngine * fem = this->fems["SolidMechanicsFEEngine"]; GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; Array Volume(this->mesh.getNbElement(element_type) * this->fems["SolidMechanicsFEEngine"]->getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem->integrate(Volume, element_type); } std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id << this->memory_id - 1; this->setBaseName (base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField ("stress" ); this->addDumpField ("grad_u" ); this->addDumpField ("eigen_grad_u" ); this->addDumpField ("blocked_dofs" ); this->addDumpField ("material_index" ); this->addDumpField ("damage" ); this->addDumpField ("Sc"); this->addDumpField ("force"); this->addDumpField ("equivalent_stress"); this->addDumpField ("residual"); this->dump(); + this->nb_dumps +=1; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions(const Matrix & displacement_gradient) { - + AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node,0); x(1) = pos(node,1); appl_disp.mul(displacement_gradient,x); (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = appl_disp(0); (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node,0); x(1) = pos(node,1); appl_disp.mul(displacement_gradient,x); (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = appl_disp(0); (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node,0); x(1) = pos(node,1); appl_disp.mul(displacement_gradient,x); (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = appl_disp(0); (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node,0); x(1) = pos(node,1); appl_disp.mul(displacement_gradient,x); (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = appl_disp(0); (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = appl_disp(1); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { + AKANTU_DEBUG_IN(); mesh.computeBoundingBox(); // find corner nodes const Array & position = mesh.getNodes(); const Vector & lower_bounds = mesh.getLowerBounds(); const Vector & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { // node 1 if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && Math::are_float_equal(position(n, 1), lower_bounds(1))) { corner_nodes(0) = n; } // node 2 else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) && Math::are_float_equal(position(n, 1), lower_bounds(1))) { corner_nodes(1) = n; } // node 3 else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) && Math::are_float_equal(position(n, 1), upper_bounds(1))) { corner_nodes(2) = n; } // node 4 else if (Math::are_float_equal(position(n, 0), lower_bounds(0)) && Math::are_float_equal(position(n, 1), upper_bounds(1))) { corner_nodes(3) = n; } } for (UInt i = 0; i < corner_nodes.getSize(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_DEBUG_ERROR("The corner node " << i + 1 << " wasn't found"); } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { - + AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal("eigen_grad_u")(element_type, gt)); Array::iterator< Matrix > prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } /// advance the damage MaterialDamageIterative<2> & mat_paste = dynamic_cast & >(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast & >(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; Real error = 0; bool converged = false; do { converged = this->solveStep<_scm_newton_raphson_tangent, _scc_increment>(1e-6, error, 2, false, *static_communicator_dummy); AKANTU_DEBUG_ASSERT(converged, "Did not converge"); std::cout << "the error is " << error << std::endl; /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) nb_damaged_elements = mat_aggregate.updateDamage(); else nb_damaged_elements = mat_paste.updateDamage(); std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); - this->dump(); + if (this->nb_dumps % 10 == 0) { + this->dump(); + this->nb_dumps += 1; + } + + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { - + AKANTU_DEBUG_IN(); FEEngine * fem = this->fems["SolidMechanicsFEEngine"]; Real average = 0; GhostType gt = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { const ElementType element_type = *it; if (field_type == "stress") { for(UInt m = 0; m < this->materials.size(); ++m) { const Array & stress_vec = this->materials[m]->getStress(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_stress"); fem->integrate(stress_vec, int_stress_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.getSize(); ++k) average += int_stress_vec(k, row_index * spatial_dimension + col_index); //3 is the value for the yy (in 3D, the value is 4) } } else if (field_type == "strain") { for(UInt m = 0; m < this->materials.size(); ++m) { const Array & gradu_vec = this->materials[m]->getGradU(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_gradu"); fem->integrate(gradu_vec, int_gradu_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.getSize(); ++k) /// averaging is done only for normal components, so stress and strain are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for(UInt m = 0; m < this->materials.size(); ++m) { const Array & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const Array & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_gradu"); fem->integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.getSize(); ++k) /// averaging is done only for normal components, so stress and strain are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else AKANTU_DEBUG_ERROR("Averaging not implemented for this field!!!"); } return average/this->volume; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { - + AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// 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.); /// 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); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { - + AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); /// solve system this->assembleStiffnessMatrix(); Real error = 0; bool converged= this->solveStep<_scm_newton_raphson_tangent_not_computed, _scc_increment>(1e-6, error, 2, false, *static_communicator_dummy); AKANTU_DEBUG_ASSERT(converged, "Did not converge"); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0,0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0,0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1,1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1,1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1,0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1,0, "strain"); - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU(Matrix & eigen_gradu_macro) { + AKANTU_DEBUG_IN(); eigen_gradu_macro(0,0) = this->averageTensorField(0,0, "eigen_grad_u"); eigen_gradu_macro(1,1) = this->averageTensorField(1,1, "eigen_grad_u"); eigen_gradu_macro(0,1) = this->averageTensorField(0,1, "eigen_grad_u"); eigen_gradu_macro(1,0) = this->averageTensorField(1,0, "eigen_grad_u"); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if(!are_materials_instantiated) instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top-bottom); Real eps = box_size * 1e-6; if(is_default_material_selector) delete material_selector; material_selector = new GelMaterialSelector(*this, box_size, "gel", this->nb_gel_pockets, eps); is_default_material_selector = false; } SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initSolver(__attribute__((unused)) SolverOptions & options) { + AKANTU_DEBUG_IN(); #if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #else jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _unsymmetric, sstr.str(), memory_id, 1); #endif //AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id, 1); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #endif //AKANTU_USE_PETSC } delete solver; std::stringstream sstr_solv; sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str(), memory_id); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str(), memory_id); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS SolverMumpsOptions opt(SolverMumpsOptions::_serial_split); if(solver) solver->initialize(opt); #endif //AKANTU_HAS_SOLVER + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension, *static_communicator_dummy); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/material_FE2/solid_mechanics_model_RVE.hh b/src/material_FE2/solid_mechanics_model_RVE.hh index 1342b8d71..a88903c1d 100644 --- a/src/material_FE2/solid_mechanics_model_RVE.hh +++ b/src/material_FE2/solid_mechanics_model_RVE.hh @@ -1,273 +1,276 @@ /** * @file solid_mechanics_model_RVE.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 14:54:18 2016 * * @brief SMM for RVE computations in FE2 simulations * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "aka_grid_dynamic.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class SolidMechanicsModelRVE : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector = true, UInt nb_gel_pockets = 400, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModelRVE(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(const ModelOptions & options = SolidMechanicsModelOptions(_static, true)); /// initialize the materials virtual void initMaterials(); /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditions(const Matrix & displacement_gradient); /// advance the reactions -> grow gel and apply homogenized properties void advanceASR(const Matrix & prestrain); /// compute average stress or strain in the model Real averageTensorField(UInt row_index, UInt col_index, const ID & field_type); /// compute effective stiffness of the RVE void homogenizeStiffness(Matrix & C_macro); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /// initialize the solver and the jacobian_matrix (called by initImplicit) virtual void initSolver(SolverOptions & options = _solver_no_options); /// allocate all vectors virtual void initArrays(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(Volume, volume, Real); private: /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ /// volume of the RVE Real volume; /// corner nodes 1, 2, 3, 4 (see Leonardo's thesis, page 98) Array corner_nodes; /// bottom nodes std::unordered_set bottom_nodes; /// left nodes std::unordered_set left_nodes; /// standard mat selector or user one bool use_RVE_mat_selector; StaticCommunicator * static_communicator_dummy; /// the number of gel pockets inside the RVE UInt nb_gel_pockets; + + /// dump counter + UInt nb_dumps; }; inline void SolidMechanicsModelRVE::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) { SolidMechanicsModel::unpackData(buffer, index, tag); if (tag == _gst_smm_uv) { Array::vector_iterator disp_it = displacement->begin(spatial_dimension); Vector current_disp(disp_it[index]); // if node is at the bottom, u_bottom = u_top +u_2 -u_3 if ( bottom_nodes.count(index) ) { current_disp += Vector(disp_it[corner_nodes(1)]); current_disp -= Vector(disp_it[corner_nodes(2)]); } // if node is at the left, u_left = u_right +u_4 -u_3 else if ( left_nodes.count(index) ) { current_disp += Vector(disp_it[corner_nodes(3)]); current_disp -= Vector(disp_it[corner_nodes(2)]); } } } /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ class GelMaterialSelector : public MeshDataMaterialSelector { public: GelMaterialSelector(SolidMechanicsModel & model, const Real box_size, const std::string & gel_material, const UInt nb_gel_pockets, Real tolerance = 0.) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), box_size(box_size) { Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = model.getSpatialDimension(); const Vector & lower_bounds = mesh.getLowerBounds(); const Vector & upper_bounds = mesh.getUpperBounds(); Vector gcenter(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { gcenter[i] = (upper_bounds[i] + lower_bounds[i]) / 2.; } Real grid_box_size = upper_bounds[0] - lower_bounds[0]; Real grid_spacing = grid_box_size/7; Vector gspacing(spatial_dimension, grid_spacing); SpatialGrid grid(spatial_dimension, gspacing, gcenter); 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,2); 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); el.element = elem; grid.insert(el, *bary_it); } /// generate the gel pockets srand(0.); Vector center(model.getSpatialDimension()); for (UInt i = 0; i < this->nb_gel_pockets; ++i) { center.clear(); center(0) = -box_size/2. + (box_size) * ((Real) rand() / (RAND_MAX)); center(1) = -box_size/2. + (box_size) * ((Real) rand() / (RAND_MAX)); Real min_dist = box_size; el.element = 0; bary_it = barycenter.begin(spatial_dimension); /// find cell in which current bary center lies SpatialGrid::CellID cell_id = grid.getCellID(center); SpatialGrid::Cell::const_iterator first_el = grid.beginCell(cell_id); SpatialGrid::Cell::const_iterator last_el = grid.endCell(cell_id); /// loop over all the elements in that cell for (;first_el != last_el; ++first_el){ const Element & elem = *first_el; Vector bary = bary_it[elem.element]; if (center.distance(bary) <= min_dist) { min_dist = center.distance(bary); el.element = elem.element; } } gel_pockets.push_back(el); } } UInt operator()(const Element & elem) { UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index == 1) return temp_index; std::vector::const_iterator iit = gel_pockets.begin(); std::vector::const_iterator eit = gel_pockets.end(); if(std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material);; } return 0; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; }; __END_AKANTU__ ///#include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_RVE_HH__ */