diff --git a/src/material_FE2/solid_mechanics_model_RVE.cc b/src/material_FE2/solid_mechanics_model_RVE.cc index 9db29c93b..1f8b05deb 100644 --- a/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,487 +1,488 @@ /** * @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 /* -------------------------------------------------------------------------- */ __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) { /// 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); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() { delete static_communicator_dummy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFull(const ModelOptions & options) { 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 << "RVE_" << 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->dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions(const Matrix & displacement_gradient) { /// 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); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { 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"); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { 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-10, error, 2, false, *static_communicator_dummy); + 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(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { 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; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { 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); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { this->applyBoundaryConditions(H); /// solve system this->assembleStiffnessMatrix(); Real error = 0; - bool converged= this->solveStep<_scm_newton_raphson_tangent_not_computed, _scc_increment>(1e-12, error, 2, false, *static_communicator_dummy); + 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"); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU(Matrix & eigen_gradu_macro) { 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"); } /* -------------------------------------------------------------------------- */ 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) { #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, *static_communicator_dummy); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS SolverMumpsOptions opt(SolverMumpsOptions::_not_parallel); if(solver) solver->initialize(opt); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ 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__