diff --git a/extra_packages/asr-tools/src/asr_tools.cc b/extra_packages/asr-tools/src/asr_tools.cc index 3913a941d..5e9b88a97 100644 --- a/extra_packages/asr-tools/src/asr_tools.cc +++ b/extra_packages/asr-tools/src/asr_tools.cc @@ -1,2315 +1,2273 @@ /** * @file ASR_tools.cc * @author Aurelia Cuba Ramos * @author Emil Gallyamov * * @brief implementation tools for the analysis of ASR samples * * @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 /* -------------------------------------------------------------------------- */ #include "aka_voigthelper.hh" #include "asr_tools.hh" #include "communicator.hh" #include "material_FE2.hh" #include "material_damage_iterative_orthotropic.hh" #include "material_iterative_stiffness_reduction.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ASRTools::ASRTools(SolidMechanicsModel & model) : model(model), volume(0.), stress_limit(0), nb_dumps(0), doubled_facets_ready(false), doubled_nodes_ready(false), node_pairs(0), disp_stored(0, model.getSpatialDimension()), ext_force_stored(0, model.getSpatialDimension()), boun_stored(0, model.getSpatialDimension()), tensile_homogenization(false) { // register event handler for asr tools auto & mesh = model.getMesh(); // auto const dim = model.getSpatialDimension(); mesh.registerEventHandler(*this, _ehp_lowest); if (mesh.hasMeshFacets()) mesh.getMeshFacets().registerEventHandler(*this, _ehp_lowest); // resize the array of tuples node_pairs.resize(2); // TODO make it dependent on number of gel sites /// find four corner nodes of the RVE findCornerNodes(); // /// resize stress limit according to the dimension size // switch (dim) { // case 2: { // stress_limit.resize(VoigtHelper<2>::size * VoigtHelper<2>::size); // break; // } // case 3: { // stress_limit.resize(VoigtHelper<3>::size * VoigtHelper<3>::size); // break; // } // } } /* -------------------------------------------------------------------------- */ void ASRTools::computeModelVolume() { auto const dim = model.getSpatialDimension(); auto & mesh = model.getMesh(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); GhostType gt = _not_ghost; this->volume = 0; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { Array Volume(mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume += fem.integrate(Volume, element_type); } // auto && comm = akantu::Communicator::getWorldCommunicator(); // std::cout << "starting reduction"; // comm.allReduce(this->volume, SynchronizerOperation::_sum); // std::cout << " ending reduction" << std::endl; } /* ------------------------------------------------------------------------- */ void ASRTools::applyFreeExpansionBC() { /// boundary conditions const auto & mesh = model.getMesh(); const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); const UInt dim = mesh.getSpatialDimension(); const Array & pos = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & boun = model.getBlockedDOFs(); /// accessing bounds Real bottom = lowerBounds(1); Real left = lowerBounds(0); Real top = upperBounds(1); Real eps = std::abs((top - bottom) * 1e-6); switch (dim) { case 2: { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; if (std::abs(pos(i, 0) - left) < eps) { boun(i, 0) = true; disp(i, 0) = 0.0; } } } break; } case 3: { /// accessing bounds Real back = lowerBounds(2); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; if ((std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 2) = 0.0; } } } break; } } } /* -------------------------------------------------------------------------- */ void ASRTools::applyLoadedBC(const Vector & traction, const ID & element_group, bool multi_axial) { /// boundary conditions const auto & mesh = model.getMesh(); const UInt dim = mesh.getSpatialDimension(); const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real left = lowerBounds(0); Real top = upperBounds(1); Real eps = std::abs((top - bottom) * 1e-6); const Array & pos = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & boun = model.getBlockedDOFs(); // disp.clear(); // boun.clear(); switch (dim) { case 2: { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; } if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } if (multi_axial && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } } break; } case 3: { Real back = lowerBounds(2); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 1) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; disp(i, 2) = 0.0; } if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; } if (multi_axial && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 2) = true; disp(i, 2) = 0.0; } if (multi_axial && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } } } } // try { model.applyBC(BC::Neumann::FromTraction(traction), element_group); // } catch (...) { // } } /* -------------------------------------------------------------------------- */ void ASRTools::fillNodeGroup(NodeGroup & node_group, bool multi_axial) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); const Vector & upperBounds = mesh.getUpperBounds(); const Vector & lowerBounds = mesh.getLowerBounds(); Real top = upperBounds(1); Real right = upperBounds(0); Real bottom = lowerBounds(1); Real front; if (dim == 3) front = upperBounds(2); Real eps = std::abs((top - bottom) * 1e-6); const Array & pos = mesh.getNodes(); if (multi_axial) { /// fill the NodeGroup with the nodes on the left, bottom and back /// surface for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 0) - right) < eps) node_group.add(i); if (std::abs(pos(i, 1) - top) < eps) node_group.add(i); if (dim == 3 && std::abs(pos(i, 2) - front) < eps) node_group.add(i); } } /// fill the NodeGroup with the nodes on the bottom surface else { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - top) < eps) node_group.add(i); } } } /* -------------------------------------------------------------------------- */ Real ASRTools::computeAverageDisplacement(SpatialDirection direction) { Real av_displ = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); const Vector & upperBounds = mesh.getUpperBounds(); const Vector & lowerBounds = mesh.getLowerBounds(); Real right = upperBounds(0); Real left = lowerBounds(0); Real top = upperBounds(1); Real front; if (dim == 3) front = upperBounds(2); Real eps = std::abs((right - left) * 1e-6); const Array & pos = mesh.getNodes(); const Array & disp = model.getDisplacement(); UInt nb_nodes = 0; if (direction == _x) { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 0) - right) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 0); ++nb_nodes; } } } else if (direction == _y) { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - top) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 1); ++nb_nodes; } } } else if ((direction == _z) && (model.getSpatialDimension() == 3)) { AKANTU_DEBUG_ASSERT(model.getSpatialDimension() == 3, "no z-direction in 2D problem"); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 2) - front) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 2); ++nb_nodes; } } } else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(av_displ, SynchronizerOperation::_sum); comm.allReduce(nb_nodes, SynchronizerOperation::_sum); return av_displ / nb_nodes; } /* -------------------------------------------------------------------------- */ Real ASRTools::computeVolumetricExpansion(SpatialDirection direction) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); UInt tensor_component = 0; if (dim == 2) { if (direction == _x) tensor_component = 0; else if (direction == _y) tensor_component = 3; else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); } else if (dim == 3) { if (direction == _x) tensor_component = 0; else if (direction == _y) tensor_component = 4; else if (direction == _z) tensor_component = 8; else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); } else AKANTU_EXCEPTION("This example does not work for 1D!!!"); Real gradu_tot = 0; Real tot_volume = 0; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const FEEngine & fe_engine = model.getFEEngine(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { const ElementTypeMapUInt & element_filter_map = model.getMaterial(m).getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = model.getMaterial(m).getElementFilter(element_type); if (!elem_filter.size()) continue; // const Array & gradu_vec = // model.getMaterial(m).getGradU(element_type); Array gradu_vec(elem_filter.size() * fe_engine.getNbIntegrationPoints(element_type), dim * dim, 0.); fe_engine.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vec, dim, element_type, gt, elem_filter); Array int_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fe_engine.integrate(gradu_vec, int_gradu_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) { gradu_tot += int_gradu_vec(k, tensor_component); } } Array Volume(mesh.getNbElement(element_type) * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); Real int_volume = fe_engine.integrate(Volume, element_type); tot_volume += int_volume; } auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(gradu_tot, SynchronizerOperation::_sum); comm.allReduce(tot_volume, SynchronizerOperation::_sum); return gradu_tot / tot_volume; } /* -------------------------------------------------------------------------- */ Real ASRTools::computeDamagedVolume(const ID & mat_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Real total_damage = 0; GhostType gt = _not_ghost; const Material & mat = model.getMaterial(mat_name); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const FEEngine & fe_engine = model.getFEEngine(); const ElementTypeMapUInt & element_filter_map = mat.getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = mat.getElementFilter(element_type); if (!elem_filter.size()) continue; const Array & damage = mat.getInternal("damage")(element_type); total_damage += fe_engine.integrate(damage, element_type, gt, elem_filter); } auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(total_damage, SynchronizerOperation::_sum); return total_damage; } /* -------------------------------------------------------------------------- */ void ASRTools::computeStiffnessReduction(std::ofstream & file_output, Real time, bool tension) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); /// save nodal values before test storeNodalFields(); if (dim == 2) { Real int_residual_x = 0; Real int_residual_y = 0; // try { int_residual_x = performLoadingTest(_x, tension); int_residual_y = performLoadingTest(_y, tension); // } catch (...) { // } if (prank == 0) file_output << time << "," << int_residual_x << "," << int_residual_y << std::endl; } else { Real int_residual_x = performLoadingTest(_x, tension); Real int_residual_y = performLoadingTest(_y, tension); Real int_residual_z = performLoadingTest(_z, tension); if (prank == 0) file_output << time << "," << int_residual_x << "," << int_residual_y << "," << int_residual_z << std::endl; } /// return the nodal values restoreNodalFields(); } -/* -------------------------------------------------------------------------- - */ -void ASRTools::computeStiffnessReductionFe2Material(std::ofstream & file_output, - Real time, bool tension) { - - /// variables for parallel execution - auto && comm = akantu::Communicator::getWorldCommunicator(); - auto prank = comm.whoAmI(); - - /// access the FE2 material and enable usage of homogenized stiffness - const UInt dim = 2; - MaterialFE2 & mat = - dynamic_cast &>(model.getMaterial("FE2_mat")); - mat.enableHomogenStiffness(); - - /// save nodal values before test - storeNodalFields(); - - if (dim == 2) { - Real int_residual_x = performLoadingTest(_x, tension); - Real int_residual_y = performLoadingTest(_y, tension); - if (prank == 0) - file_output << time << "," << int_residual_x << "," << int_residual_y - << std::endl; - } - - else { - Real int_residual_x = performLoadingTest(_x, tension); - Real int_residual_y = performLoadingTest(_y, tension); - Real int_residual_z = performLoadingTest(_z, tension); - if (prank == 0) - file_output << time << "," << int_residual_x << "," << int_residual_y - << "," << int_residual_z << std::endl; - } - - /// return the nodal values - restoreNodalFields(); - - /// disable usage of homogenized stiffness - mat.disableHomogenStiffness(); -} - /* -------------------------------------------------------------------------- */ void ASRTools::storeNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); this->disp_stored.copy(disp); this->boun_stored.copy(boun); this->ext_force_stored.copy(ext_force); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); disp.copy(this->disp_stored); boun.copy(this->boun_stored); ext_force.copy(this->ext_force_stored); /// update grad_u // model.assembleInternalForces(); } /* ---------------------------------------------------------------------- */ void ASRTools::resetNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); disp.clear(); boun.clear(); ext_force.clear(); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreInternalFields() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); mat.restorePreviousState(); } } /* -------------------------------------------------------------------------- */ void ASRTools::resetInternalFields() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); mat.resetInternalsWithHistory(); } } /* -------------------------------------------------------------------------- */ void ASRTools::storeDamageField() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.isInternal("damage_stored", _ek_regular) && mat.isInternal("reduction_step_stored", _ek_regular)) { for (auto & el_type : mat.getElementFilter().elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto & red_stored = mat.getInternal("reduction_step_stored")(el_type, _not_ghost); auto & red_current = mat.getInternal("reduction_step")(el_type, _not_ghost); red_stored.copy(red_current); auto & dam_stored = mat.getInternal("damage_stored")(el_type, _not_ghost); auto & dam_current = mat.getInternal("damage")(el_type, _not_ghost); dam_stored.copy(dam_current); } } } } /* -------------------------------------------------------------------------- */ void ASRTools::restoreDamageField() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.isInternal("damage_stored", _ek_regular) && mat.isInternal("reduction_step_stored", _ek_regular)) { for (auto & el_type : mat.getElementFilter().elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto & red_stored = mat.getInternal("reduction_step_stored")(el_type, _not_ghost); auto & red_current = mat.getInternal("reduction_step")(el_type, _not_ghost); red_current.copy(red_stored); auto & dam_stored = mat.getInternal("damage_stored")(el_type, _not_ghost); auto & dam_current = mat.getInternal("damage")(el_type, _not_ghost); dam_current.copy(dam_stored); } } } } /* -------------------------------------------------------------------------- */ Real ASRTools::performLoadingTest(SpatialDirection direction, bool tension) { UInt dir; if (direction == _x) dir = 0; if (direction == _y) dir = 1; if (direction == _z) { AKANTU_DEBUG_ASSERT(model.getSpatialDimension() == 3, "Error in problem dimension!!!"); dir = 2; } const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); /// indicate to material that its in the loading test UInt nb_materials = model.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", true); } } /// boundary conditions const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real left = lowerBounds(0); Real back; if (dim == 3) back = lowerBounds(2); Real eps = std::abs((top - bottom) * 1e-6); Real imposed_displacement = std::abs(lowerBounds(dir) - upperBounds(dir)); imposed_displacement *= 0.001; const auto & pos = mesh.getNodes(); auto & disp = model.getDisplacement(); auto & boun = model.getBlockedDOFs(); auto & ext_force = model.getExternalForce(); UInt nb_nodes = mesh.getNbNodes(); disp.clear(); boun.clear(); ext_force.clear(); if (dim == 3) { for (UInt i = 0; i < nb_nodes; ++i) { /// fix one corner node to avoid sliding if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 1) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; disp(i, 2) = 0.0; } if ((std::abs(pos(i, dir) - lowerBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = 0.0; } if ((std::abs(pos(i, dir) - upperBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = (2 * tension - 1) * imposed_displacement; } } } else { for (UInt i = 0; i < nb_nodes; ++i) { /// fix one corner node to avoid sliding if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; boun(i, 1) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; } if ((std::abs(pos(i, dir) - lowerBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = 0.0; } if ((std::abs(pos(i, dir) - upperBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = (2 * tension - 1) * imposed_displacement; } } } try { model.solveStep(); } catch (debug::Exception & e) { auto & solver = model.getNonLinearSolver("static"); int nb_iter = solver.get("nb_iterations"); std::cout << "Loading test did not converge in " << nb_iter << " iterations." << std::endl; throw e; } /// compute the force (residual in this case) along the edge of the /// imposed displacement Real int_residual = 0.; const Array & residual = model.getInternalForce(); for (UInt n = 0; n < nb_nodes; ++n) { if (std::abs(pos(n, dir) - upperBounds(dir)) < eps && mesh.isLocalOrMasterNode(n)) { int_residual += -residual(n, dir); } } /// indicate to material that its out of the loading test for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", false); } } /// restore historical internal fields restoreInternalFields(); auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(int_residual, SynchronizerOperation::_sum); return int_residual; } /* -------------------------------------------------------------------------- */ void ASRTools::computeAverageProperties(std::ofstream & file_output) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real damage_agg = computeDamagedVolume("aggregate"); Real damage_paste = computeDamagedVolume("paste"); if (dim == 2) { if (prank == 0) file_output << av_strain_x << " " << av_strain_y << " " << av_displ_x << " " << av_displ_y << " " << damage_agg << " " << damage_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << av_strain_x << " " << av_strain_y << " " << av_strain_z << " " << av_displ_x << " " << av_displ_y << " " << av_displ_z << " " << damage_agg << " " << damage_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAverageProperties(std::ofstream & file_output, Real time) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real damage_agg = computeDamagedVolume("aggregate"); Real damage_paste = computeDamagedVolume("paste"); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << damage_agg << "," << damage_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << "," << damage_agg << "," << damage_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAveragePropertiesFe2Material(std::ofstream & file_output, Real time) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real damage_paste = averageScalarField("damage_ratio_paste"); Real damage_agg = averageScalarField("damage_ratio_agg"); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << damage_agg << "," << damage_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << "," << damage_agg << "," << damage_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::saveState(UInt current_step) { /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); GhostType gt = _not_ghost; /// save the reduction step on each quad UInt nb_materials = model.getNbMaterials(); /// open the output file std::stringstream file_stream; file_stream << "./restart/reduction_steps_file_" << prank << "_" << current_step << ".txt"; std::string reduction_steps_file = file_stream.str(); std::ofstream file_output; file_output.open(reduction_steps_file.c_str()); if (!file_output.is_open()) AKANTU_EXCEPTION("Could not create the file " + reduction_steps_file + ", does its folder exist?"); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = model.getMaterial(m); if (mat.getName() == "gel") continue; /// get the reduction steps internal field const InternalField & reduction_steps = mat.getInternal("damage_step"); /// loop over all types in that material for (auto element_type : reduction_steps.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; const Array & reduction_step_array = reduction_steps(element_type, gt); for (UInt i = 0; i < reduction_step_array.size(); ++i) { file_output << reduction_step_array(i) << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + reduction_steps_file); } } } /// close the file file_output.close(); comm.barrier(); /// write the number of the last successfully saved step in a file if (prank == 0) { std::string current_step_file = "./restart/current_step.txt"; std::ofstream file_output; file_output.open(current_step_file.c_str()); file_output << current_step << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + current_step_file); file_output.close(); } } /* -------------------------------------------------------------------------- */ bool ASRTools::loadState(UInt & current_step) { current_step = 0; bool restart = false; /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); GhostType gt = _not_ghost; /// proc 0 has to save the current step if (prank == 0) { std::string line; std::string current_step_file = "./restart/current_step.txt"; std::ifstream file_input; file_input.open(current_step_file.c_str()); if (file_input.is_open()) { std::getline(file_input, line); std::stringstream sstr(line); sstr >> current_step; file_input.close(); } } comm.broadcast(current_step); if (!current_step) return restart; if (prank == 0) std::cout << "....Restarting simulation" << std::endl; restart = true; /// save the reduction step on each quad UInt nb_materials = model.getNbMaterials(); /// open the output file std::stringstream file_stream; file_stream << "./restart/reduction_steps_file_" << prank << "_" << current_step << ".txt"; std::string reduction_steps_file = file_stream.str(); std::ifstream file_input; file_input.open(reduction_steps_file.c_str()); if (!file_input.is_open()) AKANTU_EXCEPTION("Could not open file " + reduction_steps_file); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = model.getMaterial(m); if (mat.getName() == "gel") continue; /// get the material parameters Real E = mat.getParam("E"); Real max_damage = mat.getParam("max_damage"); Real max_reductions = mat.getParam("max_reductions"); /// get the internal field that need to be set InternalField & reduction_steps = const_cast &>(mat.getInternal("damage_step")); InternalField & Sc = const_cast &>(mat.getInternal("Sc")); InternalField & damage = const_cast &>(mat.getInternal("damage")); /// loop over all types in that material Real reduction_constant = mat.getParam("reduction_constant"); for (auto element_type : reduction_steps.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; Array & reduction_step_array = reduction_steps(element_type, gt); Array & Sc_array = Sc(element_type, gt); Array & damage_array = damage(element_type, gt); const Array & D_tangent = mat.getInternal("tangent")(element_type, gt); const Array & eps_u = mat.getInternal("ultimate_strain")(element_type, gt); std::string line; for (UInt i = 0; i < reduction_step_array.size(); ++i) { std::getline(file_input, line); if (file_input.fail()) AKANTU_EXCEPTION("Could not read data from file " + reduction_steps_file); std::stringstream sstr(line); sstr >> reduction_step_array(i); if (reduction_step_array(i) == 0) continue; if (reduction_step_array(i) == max_reductions) { damage_array(i) = max_damage; Real previous_damage = 1. - (1. / std::pow(reduction_constant, reduction_step_array(i) - 1)); Sc_array(i) = eps_u(i) * (1. - previous_damage) * E * D_tangent(i) / ((1. - previous_damage) * E + D_tangent(i)); } else { damage_array(i) = 1. - (1. / std::pow(reduction_constant, reduction_step_array(i))); Sc_array(i) = eps_u(i) * (1. - damage_array(i)) * E * D_tangent(i) / ((1. - damage_array(i)) * E + D_tangent(i)); } } } } /// close the file file_input.close(); return restart; } /* -------------------------------------------------------------------------- */ Real ASRTools::computePhaseVolume(const ID & mat_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Real total_volume = 0; GhostType gt = _not_ghost; const Material & mat = model.getMaterial(mat_name); for (auto element_type : mesh.elementTypes(dim, gt, _ek_regular)) { const FEEngine & fe_engine = model.getFEEngine(); const ElementTypeMapUInt & element_filter_map = mat.getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = mat.getElementFilter(element_type); if (!elem_filter.size()) continue; Array volume(elem_filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); total_volume += fe_engine.integrate(volume, element_type, gt, elem_filter); } auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(total_volume, SynchronizerOperation::_sum); return total_volume; } /* -------------------------------------------------------------------------- */ void ASRTools::applyEigenGradUinCracks( const Matrix prescribed_eigen_grad_u, const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); const Real max_damage = mat.getParam("max_damage"); const ElementTypeMapUInt & element_filter = mat.getElementFilter(); for (auto element_type : element_filter.elementTypes(dim, gt, _ek_not_defined)) { if (!element_filter(element_type, gt).size()) continue; const Array & damage = mat.getInternal("damage")(element_type, gt); const Real * dam_it = damage.storage(); Array & eigen_gradu = mat.getInternal("eigen_grad_u")(element_type, gt); Array::matrix_iterator eigen_it = eigen_gradu.begin(dim, dim); Array::matrix_iterator eigen_end = eigen_gradu.end(dim, dim); for (; eigen_it != eigen_end; ++eigen_it, ++dam_it) { if (Math::are_float_equal(max_damage, *dam_it)) { Matrix & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::applyEigenGradUinCracks( const Matrix prescribed_eigen_grad_u, const ElementTypeMapUInt & critical_elements, bool to_add) { GhostType gt = _not_ghost; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const Array & critical_elements_vec = critical_elements(element_type, gt); const Array & material_index_vec = model.getMaterialByElement(element_type, gt); const Array & material_local_numbering_vec = model.getMaterialLocalNumbering(element_type, gt); for (UInt e = 0; e < critical_elements_vec.size(); ++e) { UInt element = critical_elements_vec(e); Material & material = model.getMaterial(material_index_vec(element)); Array & eigen_gradu = material.getInternal("eigen_grad_u")(element_type, gt); UInt material_local_num = material_local_numbering_vec(element); Array::matrix_iterator eigen_it = eigen_gradu.begin(dim, dim); eigen_it += material_local_num; Matrix & current_eigengradu = *eigen_it; if (to_add) current_eigengradu += prescribed_eigen_grad_u; else current_eigengradu = prescribed_eigen_grad_u; } } } // /* // -------------------------------------------------------------------------- // */ void ASRTools::fillCracks(ElementTypeMapReal & saved_damage) { // const UInt dim = model.getSpatialDimension(); // const Material & mat_gel = model.getMaterial("gel"); // const Real E_gel = mat_gel.getParam("E"); // Real E_homogenized = 0.; // GhostType gt = _not_ghost; // const auto & mesh = model.getMesh(); // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // if (mat.getName() == "gel") // continue; // const Real E = mat.getParam("E"); // InternalField & damage = mat.getInternal("damage"); // for (auto element_type : mesh.elementTypes(dim, gt, _ek_regular)) { // const Array & elem_filter = // mat.getElementFilter(element_type, gt); if (!elem_filter.size()) // continue; // Array & damage_vec = damage(element_type, gt); // Array & saved_damage_vec = saved_damage(element_type, gt); // for (UInt i = 0; i < damage_vec.size(); ++i) { // saved_damage_vec(elem_filter(i)) = damage_vec(i); // E_homogenized = (E_gel - E) * damage_vec(i) + E; // damage_vec(i) = 1. - (E_homogenized / E); // } // } // } // } // /* // -------------------------------------------------------------------------- // */ void ASRTools::drainCracks(const ElementTypeMapReal & // saved_damage) { // // model.dump(); // const UInt dim = model.getSpatialDimension(); // const auto & mesh = model.getMesh(); // GhostType gt = _not_ghost; // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // if (mat.getName() == "gel") // continue; // else { // InternalField & damage = mat.getInternal("damage"); // for (auto element_type : mesh.elementTypes(dim, gt, // _ek_not_defined)) { // const Array & elem_filter = // mat.getElementFilter(element_type, gt); // if (!elem_filter.size()) // continue; // Array & damage_vec = damage(element_type, gt); // const Array & saved_damage_vec = // saved_damage(element_type, gt); for (UInt i = 0; i < // damage_vec.size(); ++i) { // damage_vec(i) = saved_damage_vec(elem_filter(i)); // } // } // } // } // // model.dump(); // } /* -------------------------------------------------------------------------- */ template Real ASRTools::computeSmallestElementSize() { const auto & mesh = model.getMesh(); // constexpr auto dim = model.getSpatialDimension(); /// compute smallest element size const Array & pos = mesh.getNodes(); Real el_h_min = 10000.; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_element = mesh.getNbElement(element_type); Array X(0, nb_nodes_per_element * dim); model.getFEEngine().extractNodalToElementField(mesh, pos, X, element_type, _not_ghost); Array::matrix_iterator X_el = X.begin(dim, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el) { Real el_h = model.getFEEngine().getElementInradius(*X_el, element_type); if (el_h < el_h_min) el_h_min = el_h; } } /// find global Gauss point with highest stress auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(el_h_min, SynchronizerOperation::_min); return el_h_min; } // /* // -------------------------------------------------------------------------- // */ // /// apply homogeneous temperature on the whole specimen // void ASRTools::applyTemperatureFieldToSolidmechanicsModel( // const Real & temperature) { // UInt dim = model.getMesh().getSpatialDimension(); // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // for (auto el_type : mat.getElementFilter().elementTypes(dim)) { // const auto & filter = mat.getElementFilter()(el_type); // if (filter.size() == 0) // continue; // auto & delta_T = mat.getArray("delta_T", el_type); // auto dt = delta_T.begin(); // auto dt_end = delta_T.end(); // for (; dt != dt_end; ++dt) { // *dt = temperature; // } // } // } // } /* -------------------------------------------------------------------------- */ Real ASRTools::computeDeltaGelStrainThermal(const Real delta_time_day, const Real k, const Real activ_energy, const Real R, const Real T, Real & amount_reactive_particles, const Real saturation_const) { /// compute increase in gel strain value for interval of time delta_time /// as temperatures are stored in C, conversion to K is done Real delta_strain = amount_reactive_particles * k * std::exp(-activ_energy / (R * (T + 273.15))) * delta_time_day; amount_reactive_particles -= std::exp(-activ_energy / (R * (T + 273.15))) * delta_time_day / saturation_const; if (amount_reactive_particles < 0.) amount_reactive_particles = 0.; return delta_strain; } /* -----------------------------------------------------------------------*/ Real ASRTools::computeDeltaGelStrainLinear(const Real delta_time, const Real k) { /// compute increase in gel strain value for dt simply by deps = k * /// delta_time Real delta_strain = k * delta_time; return delta_strain; } /* ---------------------------------------------------------------------- */ void ASRTools::applyBoundaryConditionsRve( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // /// transform into Green strain to exclude rotations // auto F = displacement_gradient + Matrix::eye(dim); // auto E = 0.5 * (F.transpose() * F - Matrix::eye(dim)); /// 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(dim); Vector appl_disp(dim); const auto & lower_bounds = mesh.getLowerBounds(); for (auto node : this->corner_nodes) { x(0) = pos(node, 0); x(1) = pos(node, 1); x -= lower_bounds; appl_disp.mul(displacement_gradient, x); (model.getBlockedDOFs())(node, 0) = true; (model.getDisplacement())(node, 0) = appl_disp(0); (model.getBlockedDOFs())(node, 1) = true; (model.getDisplacement())(node, 1) = appl_disp(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::applyHomogeneousTemperature(const Real & temperature) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); for (auto el_type : mat.getElementFilter().elementTypes(dim)) { const auto & filter = mat.getElementFilter()(el_type); if (filter.size() == 0) continue; auto & deltas_T = mat.getArray("delta_T", el_type); for (auto && delta_T : deltas_T) { delta_T = temperature; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::removeTemperature() { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); for (auto el_type : mat.getElementFilter().elementTypes(dim)) { const auto & filter = mat.getElementFilter()(el_type); if (filter.size() == 0) continue; auto & deltas_T = mat.getArray("delta_T", el_type); for (auto && delta_T : deltas_T) { delta_T = 0.; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::findCornerNodes() { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); switch (dim) { case 2: { corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, dim))) { 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; } } break; } case 3: { corner_nodes.resize(8); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, dim))) { 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)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { 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)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { 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)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { corner_nodes(3) = node; } // node 5 if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(4) = node; } // node 6 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(5) = node; } // node 7 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(6) = node; } // node 8 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(7) = node; } } break; } } // AKANTU_DEBUG_ASSERT(dim == 2, "This is 2D only!"); // for (UInt i = 0; i < corner_nodes.size(); ++i) { // if (corner_nodes(i) == UInt(-1)) // AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real ASRTools::averageScalarField(const ID & field_name) { AKANTU_DEBUG_IN(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); Real average = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); if (!elem_filter.size()) continue; const auto & scalar_field = model.getMaterial(m).getInternal(field_name)(element_type); Array int_scalar_vec(elem_filter.size(), 1, "int_of_scalar"); fem.integrate(scalar_field, int_scalar_vec, 1, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_scalar_vec(k); } } /// compute total model volume if (!this->volume) computeModelVolume(); return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real ASRTools::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); Real average = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & stress_vec = model.getMaterial(m).getStress(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), dim * dim, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec(k, row_index * dim + 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 < model.getNbMaterials(); ++m) { const auto & gradu_vec = model.getMaterial(m).getGradU(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, dim * dim, 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 * dim + col_index) + int_gradu_vec(k, col_index * dim + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & eigen_gradu_vec = model.getMaterial(m).getInternal( "eigen_grad_u")(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, dim * dim, 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 * dim + col_index); } } else { AKANTU_ERROR("Averaging not implemented for this field!!!"); } } /// compute total model volume if (!this->volume) computeModelVolume(); return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::homogenizeStressField(Matrix & stress) { AKANTU_DEBUG_IN(); stress(0, 0) = averageTensorField(0, 0, "stress"); stress(1, 1) = averageTensorField(1, 1, "stress"); stress(0, 1) = averageTensorField(0, 1, "stress"); stress(1, 0) = averageTensorField(1, 0, "stress"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::setStiffHomogenDir(Matrix & stress) { AKANTU_DEBUG_IN(); auto dim = stress.rows(); Vector eigenvalues(dim); stress.eig(eigenvalues); Real hydrostatic_stress = 0; UInt denominator = 0; for (UInt i = 0; i < dim; ++i, ++denominator) hydrostatic_stress += eigenvalues(i); hydrostatic_stress /= denominator; AKANTU_DEBUG_OUT(); this->tensile_homogenization = (hydrostatic_stress > 0); } /* -------------------------------------------------------------------------- */ void ASRTools::homogenizeStiffness(Matrix & C_macro, bool tensile_test) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim == 2, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (0.001;0;0) 2. eps_el = (0,0.001,0) 3. eps_el = (0,0,0.001) /// clear the eigenstrain (to exclude stresses due to internal pressure) clearGelEigenStrain(); /// save nodal values before tests storeNodalFields(); /// storage for results of 3 different loading states UInt voigt_size = 1; switch (dim) { case 2: { voigt_size = VoigtHelper<2>::size; break; } case 3: { voigt_size = VoigtHelper<3>::size; break; } } 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); /// indicate to material that its in the loading test UInt nb_materials = model.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", true); } } /// virtual test 1: H(0, 0) = 0.001 * (2 * tensile_test - 1); performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); H(1, 1) = 0.001 * (2 * tensile_test - 1); performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); H(0, 1) = 0.001; H(1, 0) = 0.001; performVirtualTesting(H, stresses, strains, 2); /// indicate to material that its out of the loading test for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", false); } } // /// set up the stress limit at 10% of stresses in undamaged state // if (first_time) { // auto && str_lim_mat_it = // make_view(this->stress_limit, voigt_size, voigt_size).begin(); // *str_lim_mat_it = stresses * 0.1; // } // /// compare stresses with the lower limit and update if needed // auto && str_lim_vec_it = make_view(this->stress_limit, // voigt_size).begin(); for (UInt i = 0; i != voigt_size; ++i, // ++str_lim_vec_it) { // Vector stress = stresses(i); // // Vector stress_lim = stress_limit(i); // Real stress_norm = stress.norm(); // Real stress_limit_norm = (*str_lim_vec_it).norm(); // if (stress_norm < stress_limit_norm) // // stresses(i) = *str_lim_vec_it; // return; // } /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); /// Make C matrix symmetric Matrix C_direct(voigt_size, voigt_size); C_direct.mul(stresses, eps_inverse); for (UInt i = 0; i != voigt_size; ++i) { for (UInt j = 0; j != voigt_size; ++j) { C_macro(i, j) = 0.5 * (C_direct(i, j) + C_direct(j, i)); C_macro(j, i) = C_macro(i, j); } } /// return the nodal values restoreNodalFields(); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void ASRTools::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); auto & disp = model.getDisplacement(); auto & boun = model.getBlockedDOFs(); auto & ext_force = model.getExternalForce(); disp.clear(); boun.clear(); ext_force.clear(); applyBoundaryConditionsRve(H); // auto & solver = model.getNonLinearSolver(); // solver.set("max_iterations", 50); // solver.set("threshold", 1e-5); // solver.set("convergence_type", SolveConvergenceCriteria::_solution); model.solveStep(); /// get average stress and strain eff_stresses(0, test_no) = averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * averageTensorField(1, 0, "strain"); /// restore historical internal fields restoreInternalFields(); AKANTU_DEBUG_OUT(); } /* --------------------------------------------------- */ void ASRTools::homogenizeEigenGradU(Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = model.getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; Real E = mat.get("E"); auto & damage = mat.getInternal("damage"); GhostType gt = _not_ghost; for (auto && type : damage.elementTypes(gt)) { const auto & elem_filter = mat.getElementFilter(type); auto nb_integration_point = model.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 ASRTools::drainCracks(const ElementTypeMapReal & saved_damage) { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; auto & damage = mat.getInternal("damage"); GhostType gt = _not_ghost; for (auto && type : damage.elementTypes(gt)) { const auto & elem_filter = mat.getElementFilter(type); auto nb_integration_point = model.getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::computeDamageRatio(Real & damage_ratio) { damage_ratio = 0.; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; GhostType gt = _not_ghost; const ElementTypeMapArray & filter_map = mat.getElementFilter(); // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; const FEEngine & fe_engine = model.getFEEngine(); auto & damage_array = mat.getInternal("damage")(element_type); damage_ratio += fe_engine.integrate(damage_array, element_type, gt, filter); } } /// compute total model volume if (!this->volume) computeModelVolume(); damage_ratio /= this->volume; } /* -------------------------------------------------------------------------- */ void ASRTools::computeDamageRatioPerMaterial(Real & damage_ratio, const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); const ElementTypeMapArray & filter_map = mat.getElementFilter(); const FEEngine & fe_engine = model.getFEEngine(); damage_ratio = 0.; Real mat_volume = 0.; // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; auto & damage_array = mat.getInternal("damage")(element_type); damage_ratio += fe_engine.integrate(damage_array, element_type, gt, filter); Array volume( filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); mat_volume += fe_engine.integrate(volume, element_type, gt, filter); } damage_ratio /= mat_volume; } /* -------------------------------------------------------------------------*/ void ASRTools::dumpRve() { // if (this->nb_dumps % 10 == 0) { model.dump(); // } // this->nb_dumps += 1; } /* -------------------------------------------------------------------------- */ void ASRTools::applyBodyForce() { auto spatial_dimension = model.getSpatialDimension(); model.assembleMassLumped(); auto & mass = model.getMass(); auto & force = model.getExternalForce(); Vector gravity(spatial_dimension); gravity(1) = -9.81; for (auto && data : zip(make_view(mass, spatial_dimension), make_view(force, spatial_dimension))) { const auto & mass_vec = (std::get<0>(data)); auto & force_vec = (std::get<1>(data)); force_vec += gravity * mass_vec; } } /* ------------------------------------------------------------------------- */ void ASRTools::insertCohElemByCoord(const Vector & position) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); const auto gt = _not_ghost; // insert cohesive elements auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh_facets = inserter.getMeshFacets(); Vector bary_facet(dim); Real min_dist = std::numeric_limits::max(); Element source_facet; for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto dist = bary_facet.distance(position); if (dist < min_dist) { min_dist = dist; source_facet = facet; } }, _spatial_dimension = dim - 1); inserter.getCheckFacets(source_facet.type, gt)(source_facet.element) = false; insertion(source_facet) = true; inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertCohElemByLimits(const Matrix & insertion_limits, std::string coh_mat_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); auto pos_it = make_view(mesh.getNodes(), dim).begin(); // insert cohesive elements auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh_facets = inserter.getMeshFacets(); auto coh_material = model.getMaterialIndex(coh_mat_name); auto tolerance = Math::getTolerance(); Vector bary_facet(dim); for_each_element( mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < dim and bary_facet(coord_in_limit) > (insertion_limits(coord_in_limit, 0) - tolerance) and bary_facet(coord_in_limit) < (insertion_limits(coord_in_limit, 1) + tolerance)) ++coord_in_limit; if (coord_in_limit == dim) { coh_model.getFacetMaterial(facet.type, facet.ghost_type)( facet.element) = coh_material; inserter.getCheckFacets(facet.type, facet.ghost_type)(facet.element) = false; insertion(facet) = true; } }, _spatial_dimension = dim - 1); inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertCohElemRandomly(const UInt & nb_coh_elem, std::string coh_mat_name, std::string matrix_mat_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); auto & coh_model = dynamic_cast(model); // get element types const GhostType gt = _not_ghost; auto type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto type_facets = Mesh::getFacetType(type); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); Vector bary_facet(dim); auto coh_material = model.getMaterialIndex(coh_mat_name); auto & mesh_facets = inserter.getMeshFacets(); auto matrix_material_id = model.getMaterialIndex(matrix_mat_name); auto & matrix_elements = mesh_facets.createElementGroup("matrix_facets", dim - 1); for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto & facet_material = coh_model.getFacetMaterial( facet.type, facet.ghost_type)(facet.element); if (facet_material == matrix_material_id) matrix_elements.add(facet); }, _spatial_dimension = dim - 1); // Will be used to obtain a seed for the random number engine std::random_device rd; // Standard mersenne_twister_engine seeded with rd() std::mt19937 random_generator(rd()); std::uniform_int_distribution<> dis( 0, matrix_elements.getElements(type_facets).size() - 1); for (auto _[[gnu::unused]] : arange(nb_coh_elem)) { auto id = dis(random_generator); Element facet; facet.type = type_facets; facet.element = matrix_elements.getElements(type_facets)(id); facet.ghost_type = gt; coh_model.getFacetMaterial(facet.type, facet.ghost_type)(facet.element) = coh_material; inserter.getCheckFacets(facet.type, facet.ghost_type)(facet.element) = false; insertion(facet) = true; } inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ const Array ASRTools::insertCohElOrFacetsByCoord(const Vector & position, bool add_neighbors, bool only_double_facets) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); const auto gt = _not_ghost; const auto & pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); auto & node_group = mesh.createNodeGroup("doubled_nodes"); // insert cohesive elements auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh_facets = inserter.getMeshFacets(); Vector bary_facet(dim); Real min_dist = std::numeric_limits::max(); Element cent_facet; for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto dist = std::abs(bary_facet.distance(position)); if (dist < min_dist) { min_dist = dist; cent_facet = facet; } }, _spatial_dimension = dim - 1); inserter.getCheckFacets(cent_facet.type, gt)(cent_facet.element) = false; insertion(cent_facet) = true; Array facets_added; facets_added.push_back(cent_facet); if (add_neighbors) { auto & facet_conn = mesh_facets.getConnectivity(cent_facet.type, gt); CSR nodes_to_elements; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_elements, dim - 1); // Synchronizer element connected to the nodes of the facet for (auto node : arange(2)) { // add corner nodes to the group node_group.add(facet_conn(cent_facet.element, node)); // add corner node to the tuples std::get<0>(node_pairs(node)) = facet_conn(cent_facet.element, node); // vector of the central facet Vector cent_facet_dir(pos_it[facet_conn(cent_facet.element, !node)], true); cent_facet_dir -= Vector(pos_it[facet_conn(cent_facet.element, node)]); cent_facet_dir /= cent_facet_dir.norm(); Real min_dot = std::numeric_limits::max(); Element neighbor; Vector neighbor_facet_dir(dim); for (auto & elem : nodes_to_elements.getRow(facet_conn(cent_facet.element, node))) { if (elem.element == cent_facet.element) continue; if (elem.type != cent_facet.type) continue; // vector of the neighboring facet if (facet_conn(elem.element, 0) == facet_conn(cent_facet.element, node)) { neighbor_facet_dir = pos_it[facet_conn(elem.element, 1)]; } else if (facet_conn(elem.element, 1) == facet_conn(cent_facet.element, node)) { neighbor_facet_dir = pos_it[facet_conn(elem.element, 0)]; } else AKANTU_EXCEPTION("Neighboring facet doesn't have node in common with " "the central one."); neighbor_facet_dir -= Vector(pos_it[facet_conn(cent_facet.element, node)]); neighbor_facet_dir /= neighbor_facet_dir.norm(); Real dot = cent_facet_dir.dot(neighbor_facet_dir); if (dot < min_dot) { min_dot = dot; neighbor = elem; } } inserter.getCheckFacets(neighbor.type, gt)(neighbor.element) = false; insertion(neighbor) = true; facets_added.push_back(neighbor); } } inserter.insertElements(only_double_facets); // flag internal facets auto & doubled_facets = mesh_facets.getData("doubled_facets"); doubled_facets.initialize(mesh_facets, _spatial_dimension = dim - 1, _with_nb_element = true, _default_value = false); for (auto & facet : facets_added) { Array & doubled_facets_array = doubled_facets(facet.type, facet.ghost_type); doubled_facets_array(facet.element) = true; } // create empty element group with nodes to apply Dirichlet model.getMesh().createElementGroupFromNodeGroup("doubled_nodes", "doubled_nodes", dim - 1); AKANTU_DEBUG_OUT(); return facets_added; } /* -------------------------------------------------------------------------- */ void ASRTools::onElementsAdded(const Array & elements, const NewElementsEvent &) { if (this->doubled_facets_ready) return; auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); auto & mesh_facets = mesh.getMeshFacets(); auto & doubled_facets = mesh_facets.getData("doubled_facets"); doubled_facets.initialize(mesh_facets, _spatial_dimension = dim - 1, _with_nb_element = true, _default_value = false); for (auto elements_range : akantu::MeshElementsByTypes(elements)) { auto type = elements_range.getType(); auto ghost_type = elements_range.getGhostType(); if (mesh.getKind(type) != _ek_regular) continue; if (ghost_type != _not_ghost) continue; auto & doubled_facets_array = mesh_facets.getData("doubled_facets", type, ghost_type); auto & element_ids = elements_range.getElements(); for (auto && el : element_ids) { doubled_facets_array(el) = true; } } this->doubled_facets_ready = true; } /* -------------------------------------------------------------------------- */ void ASRTools::onNodesAdded(const Array & new_nodes, const NewNodesEvent &) { AKANTU_DEBUG_IN(); if (new_nodes.size() == 0) return; if (this->doubled_nodes_ready) return; auto & mesh = model.getMesh(); auto & node_group = mesh.getNodeGroup("doubled_nodes"); auto pos_it = make_view(mesh.getNodes(), mesh.getSpatialDimension()).begin(); for (auto & node : new_nodes) { node_group.add(node); // complete pairs with lower nodes Vector lower_coord = pos_it[node]; for (auto & pair : node_pairs) { Vector upper_coord = pos_it[std::get<0>(pair)]; if (upper_coord == lower_coord) std::get<1>(pair) = node; } } this->doubled_nodes_ready = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::applyDeltaU(Real delta_u) { // get element group with nodes to apply Dirichlet auto & crack_facets = model.getMesh().getElementGroup("doubled_nodes"); DeltaU delta_u_bc(model, delta_u, getNodePairs()); model.applyBC(delta_u_bc, crack_facets); } /* -------------------------------------------------------------------------- */ void ASRTools::applyGelStrain(const Matrix & prestrain) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(dim, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::clearGelEigenStrain() { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); 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 &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/asr-tools/src/asr_tools.hh b/extra_packages/asr-tools/src/asr_tools.hh index ac9659d45..0dedd4a90 100644 --- a/extra_packages/asr-tools/src/asr_tools.hh +++ b/extra_packages/asr-tools/src/asr_tools.hh @@ -1,594 +1,590 @@ /** * @file asr_tools.hh * @author Aurelia Cuba Ramos * @date Tue Jan 16 10:26:53 2014 * * @brief tools for the analysis of ASR samples * * @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 #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #ifdef AKANTU_COHESIVE_ELEMENT #include "solid_mechanics_model_cohesive.hh" #endif #ifdef AKANTU_FLUID_DIFFUSION #include "fluid_diffusion_model.hh" #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ASR_TOOLS_HH__ #define __AKANTU_ASR_TOOLS_HH__ namespace akantu { class NodeGroup; class SolidMechanicsModel; } // namespace akantu namespace akantu { class ASRTools : public MeshEventHandler { public: ASRTools(SolidMechanicsModel & model); virtual ~ASRTools() = default; /// This function is used in case of uni-axial boundary conditions: /// It detects the nodes, on which a traction needs to be applied and stores /// them in a node group void fillNodeGroup(NodeGroup & node_group, bool multi_axial = false); /// compute volume of the model passed to ASRTools (RVE if FE2_mat) void computeModelVolume(); /// Apply free expansion boundary conditions void applyFreeExpansionBC(); /// Apply loaded boundary conditions void applyLoadedBC(const Vector & traction, const ID & element_group, bool multi_axial); /// This function computes the average displacement along one side of the /// sample, /// caused by the expansion of gel pockets Real computeAverageDisplacement(SpatialDirection direction); /// This function computes a compoment of average volumetric strain tensor, /// i.e. eps_macro_xx or eps_macro_yy or eps_macro_zz Real computeVolumetricExpansion(SpatialDirection direction); /// This function computes the amount of damage in one material phase Real computeDamagedVolume(const ID & mat_name); /// This function is used to compute the average stiffness by performing a /// virtual loading test Real performLoadingTest(SpatialDirection direction, bool tension); /// perform tension tests and integrate the internal force on the upper /// surface void computeStiffnessReduction(std::ofstream & file_output, Real time, bool tension = true); /// just the average properties, NO tension test void computeAverageProperties(std::ofstream & file_output); /// Same as the last one but writes a csv with first column having time in /// days void computeAverageProperties(std::ofstream & file_output, Real time); /// Averages strains & displacements + damage ratios in FE2 material void computeAveragePropertiesFe2Material(std::ofstream & file_output, Real time); - /// perform tension tests of a specimen composed of fe2 material - void computeStiffnessReductionFe2Material(std::ofstream & file_output, - Real time, bool tension = true); - /// Save the state of the simulation for subsequent restart void saveState(UInt current_step); /// Load to previous state of the simulation for restart bool loadState(UInt & current_step); /// compute the volume occupied by a given material Real computePhaseVolume(const ID & mat_name); /// apply eigenstrain in cracks, that are filled with gel void applyEigenGradUinCracks(const Matrix prescribed_eigen_grad_u, const ID & material_name); /// apply eigenstrain in the cracks, that advanced in the last step void applyEigenGradUinCracks(const Matrix prescribed_eigen_grad_u, const ElementTypeMapUInt & critical_elements, bool to_add = false); /// fill fully cracked elements with gel void fillCracks(ElementTypeMapReal & saved_damage); /// drain the gel in fully cracked elements void drainCracks(const ElementTypeMapReal & saved_damage); template Real computeSmallestElementSize(); // /// apply homogeneous temperature on Solid Mechanics Model // void applyTemperatureFieldToSolidmechanicsModel(const Real & temperature); /// compute increase in gel strain within 1 timestep Real computeDeltaGelStrainThermal(const Real delta_time_day, const Real k, const Real activ_energy, const Real R, const Real T, Real & amount_reactive_particles, const Real saturation_const); /// compute linear increase in gel strain Real computeDeltaGelStrainLinear(const Real delta_time, const Real k); /// insert single cohesive element by the coordinate of its center void insertCohElemByCoord(const Vector & position); /// insert multiple cohesive elements by the limiting box void insertCohElemByLimits(const Matrix & insertion_limits, std::string coh_mat_name); /// insert multiple cohesive elements by the limiting box void insertCohElemRandomly(const UInt & nb_coh_elem, std::string coh_mat_name, std::string matrix_mat_name); /// insert up to 3 facets pair based on the coord of the central one const Array insertCohElOrFacetsByCoord(const Vector & position, bool add_neighbors = true, bool only_double_facets = false); /// on elements added for asr-tools void onElementsAdded(const Array & elements, const NewElementsEvent & element_event); void onNodesAdded(const Array & new_nodes, const NewNodesEvent & nodes_event); /// apply self-weight force void applyBodyForce(); /// apply delta u on nodes void applyDeltaU(Real delta_u); /// apply eigenstrain on gel material void applyGelStrain(const Matrix & prestrain); /* ------------------------------------------------------------------------ */ /// RVE part /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditionsRve(const Matrix & displacement_gradient); /// apply homogeneous temperature field from the macroscale level to the RVEs virtual void applyHomogeneousTemperature(const Real & temperature); /// remove temperature from RVE on the end of ASR advancement virtual void removeTemperature(); /// averages scalar field over the WHOLE!!! volume of a model Real averageScalarField(const ID & field_name); /// 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 (in tension by default) void homogenizeStiffness(Matrix & C_macro, bool tensile_test = true); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /// compute damage to the RVE area ratio void computeDamageRatio(Real & damage); /// compute damage to material area ratio void computeDamageRatioPerMaterial(Real & damage_ratio, const ID & material_name); /// dump the RVE void dumpRve(); /// compute average stress in the RVE void homogenizeStressField(Matrix & stress); /// compute hydrostatic part of the stress and assign homogen direction void setStiffHomogenDir(Matrix & stress); /// storing nodal fields before tests void storeNodalFields(); /// restoring nodal fields before tests void restoreNodalFields(); /// resetting nodal fields void resetNodalFields(); /// restoring internal fields after tests void restoreInternalFields(); /// resetting internal fields with history void resetInternalFields(); /// store damage in materials who have damage void storeDamageField(); /// assign stored values to the current ones void restoreDamageField(); /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); /// clear the eigenstrain (to exclude stresses due to internal pressure) void clearGelEigenStrain(); /* ------------------------------------------------------------------------ */ public: // Accessors inline Array> getNodePairs() const { return node_pairs; } bool isTensileHomogen() { return this->tensile_homogenization; }; /* --------------------------------------------------------------------- */ /* Members */ /* --------------------------------------------------------------------- */ private: SolidMechanicsModel & model; // /// 2D hardcoded - no 3D support currently // using voigt_h = VoigtHelper<2>; protected: /// 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; /// lower limit for stresses Array stress_limit; /// dump counter UInt nb_dumps; /// booleans for applying delta u bool doubled_facets_ready; bool doubled_nodes_ready; // array of tuples to store nodes pairs:1st- is the one on the upper facet Array> node_pairs; // arrays to store nodal values during virtual tests Array disp_stored; Array ext_force_stored; Array boun_stored; /// if stiffness homogenization will be done in tension bool tensile_homogenization{false}; }; /* -------------------------------------------------------------------------- */ /* 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, std::string aggregate_material = "aggregate", 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), aggregate_material(aggregate_material) {} void initGelPocket() { aggregate_material_id = model.getMaterialIndex(aggregate_material); Mesh & mesh = this->model.getMesh(); UInt dim = model.getSpatialDimension(); // Element el{_triangle_3, 0, _not_ghost}; for (auto el_type : model.getMaterial(aggregate_material) .getElementFilter() .elementTypes(dim)) { const auto & filter = model.getMaterial(aggregate_material).getElementFilter()(el_type); if (!filter.size() == 0) AKANTU_EXCEPTION("Check the element type for aggregate material"); Element el{el_type, 0, _not_ghost}; UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); Array barycenter(nb_element, dim); for (auto && data : enumerate(make_view(barycenter, dim))) { el.element = std::get<0>(data); auto & bary = std::get<1>(data); mesh.getBarycenter(el, bary); } /// generate the gel pockets srand(0.); Vector center(dim); 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) == aggregate_material_id) { gel_pockets.push_back(el); placed_gel_pockets += 1; } } } is_gel_initialized = true; } UInt operator()(const Element & elem) { /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (not is_gel_initialized) initGelPocket(); UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index != aggregate_material_id) return temp_index; auto iit = gel_pockets.begin(); auto eit = gel_pockets.end(); if (std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; if (prank == 0) std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material); } return temp_index; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; std::string aggregate_material{"aggregate"}; UInt aggregate_material_id{1}; bool is_gel_initialized{false}; }; /* -------------------------------------------------------------------------- */ /* Boundary conditions functors */ /* -------------------------------------------------------------------------- */ class Pressure : public BC::Neumann::NeumannFunctor { public: Pressure(SolidMechanicsModel & model, const Array & pressure_on_qpoint, const Array & quad_coords) : model(model), pressure_on_qpoint(pressure_on_qpoint), quad_coords(quad_coords) {} inline void operator()(const IntegrationPoint & quad_point, Vector & dual, const Vector & coord, const Vector & normal) const { // get element types auto && mesh = model.getMesh(); const UInt dim = mesh.getSpatialDimension(); const GhostType gt = akantu::_not_ghost; const UInt facet_nb = quad_point.element; const ElementKind cohesive_kind = akantu::_ek_cohesive; const ElementType type_facet = quad_point.type; const ElementType type_coh = FEEngine::getCohesiveElementType(type_facet); auto && cohesive_conn = mesh.getConnectivity(type_coh, gt); const UInt nb_nodes_coh_elem = cohesive_conn.getNbComponent(); auto && facet_conn = mesh.getConnectivity(type_facet, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); auto && fem_boundary = model.getFEEngineBoundary(); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type_facet, gt); auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); AKANTU_DEBUG_ASSERT(nb_nodes_coh_elem == 2 * nb_nodes_facet, "Different number of nodes belonging to one cohesive " "element face and facet"); // const akantu::Mesh &mesh_facets = cohesive_mesh.getMeshFacets(); const Array> & elem_to_subelem = mesh.getElementToSubelement(type_facet, gt); const auto & adjacent_elems = elem_to_subelem(facet_nb); auto normal_corrected = normal; // loop over all adjacent elements UInt coh_elem_nb; if (not adjacent_elems.empty()) { for (UInt f : arange(adjacent_elems.size())) { if (adjacent_elems[f].kind() != cohesive_kind) continue; coh_elem_nb = adjacent_elems[f].element; Array upper_nodes(nb_nodes_coh_elem / 2); Array lower_nodes(nb_nodes_coh_elem / 2); for (UInt node : arange(nb_nodes_coh_elem / 2)) { upper_nodes(node) = cohesive_conn(coh_elem_nb, node); lower_nodes(node) = cohesive_conn(coh_elem_nb, node + nb_nodes_coh_elem / 2); } bool upper_face = true; bool lower_face = true; Vector facet_nodes = facet_nodes_it[facet_nb]; for (UInt s : arange(nb_nodes_facet)) { auto idu = upper_nodes.find(facet_nodes(s)); auto idl = lower_nodes.find(facet_nodes(s)); if (idu == UInt(-1)) upper_face = false; else if (idl == UInt(-1)) lower_face = false; } if (upper_face && not lower_face) normal_corrected *= -1; else if (not upper_face && lower_face) normal_corrected *= 1; else AKANTU_EXCEPTION("Error in defining side of the cohesive element"); break; } } auto flow_qpoint_it = make_view(quad_coords, dim).begin(); bool node_found; for (auto qp : arange(nb_quad_points)) { const Vector flow_quad_coord = flow_qpoint_it[coh_elem_nb * nb_quad_points + qp]; if (flow_quad_coord != coord) continue; Real P = pressure_on_qpoint(coh_elem_nb * nb_quad_points + qp); // if (P < 0) // P = 0.; dual = P * normal_corrected; node_found = true; } if (not node_found) AKANTU_EXCEPTION("Quad point is not found in the flow mesh"); } protected: SolidMechanicsModel & model; const Array & pressure_on_qpoint; const Array & quad_coords; }; /* -------------------------------------------------------------------------- */ class DeltaU : public BC::Dirichlet::DirichletFunctor { public: DeltaU(const SolidMechanicsModel & model, const Real delta_u, const Array> & node_pairs) : model(model), delta_u(delta_u), node_pairs(node_pairs) { displacement = model.getDisplacement(); } inline void operator()(UInt node, Vector & flags, Vector & primal, const Vector &) const { // get element types auto && mesh = model.getMesh(); const UInt dim = mesh.getSpatialDimension(); auto && mesh_facets = mesh.getMeshFacets(); auto disp_it = make_view(displacement, dim).begin(); CSR nodes_to_elements; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_elements, dim - 1); // get actual distance between two nodes Vector node_disp(disp_it[node]); Vector other_node_disp(dim); bool upper_face = false; bool lower_face = false; for (auto && pair : this->node_pairs) { if (node == std::get<0>(pair)) { other_node_disp = disp_it[std::get<1>(pair)]; upper_face = true; break; } else if (node == std::get<1>(pair)) { other_node_disp = disp_it[std::get<0>(pair)]; lower_face = true; break; } } AKANTU_DEBUG_ASSERT(upper_face == true or lower_face == true, "Error in identifying the node in tuple"); Real sign = -upper_face + lower_face; // compute normal at node (averaged between two surfaces) Vector normal(dim); for (auto & elem : nodes_to_elements.getRow(node)) { if (mesh.getKind(elem.type) != _ek_regular) continue; if (elem.ghost_type != _not_ghost) continue; auto & doubled_facets_array = mesh_facets.getData( "doubled_facets", elem.type, elem.ghost_type); if (doubled_facets_array(elem.element) != true) continue; auto && fe_engine_facet = model.getFEEngine("FacetsFEEngine"); auto nb_qpoints_per_facet = fe_engine_facet.getNbIntegrationPoints(elem.type, elem.ghost_type); const auto & normals_on_quad = fe_engine_facet.getNormalsOnIntegrationPoints(elem.type, elem.ghost_type); auto normals_it = make_view(normals_on_quad, dim).begin(); normal += sign * Vector(normals_it[elem.element * nb_qpoints_per_facet]); } normal /= normal.norm(); // get distance between two nodes in normal direction Real node_disp_norm = node_disp.dot(normal); Real other_node_disp_norm = other_node_disp.dot(-1. * normal); Real dist = node_disp_norm + other_node_disp_norm; Real prop_factor = dist == 0. ? 0.5 : node_disp_norm / dist; // get correction displacement Real correction = delta_u - dist; // apply absolute value of displacement primal += normal * correction * prop_factor; flags.set(false); } protected: const SolidMechanicsModel & model; const Real delta_u; const Array> node_pairs; Array displacement; }; } // namespace akantu #endif /* __AKANTU_ASR_TOOLS_HH__ */ diff --git a/extra_packages/asr-tools/src/material_FE2/material_FE2.cc b/extra_packages/asr-tools/src/material_FE2/material_FE2.cc index b3f7cff5d..17ab99d09 100644 --- a/extra_packages/asr-tools/src/material_FE2/material_FE2.cc +++ b/extra_packages/asr-tools/src/material_FE2/material_FE2.cc @@ -1,654 +1,603 @@ /** * @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 "aka_iterators.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), gelstrain("gelstrain", *this), non_reacted_gel("non_reacted_gel", *this), damage_ratio("damage_ratio", *this), damage_ratio_paste("damage_ratio_paste", *this), damage_ratio_agg("damage_ratio_agg", *this) { AKANTU_DEBUG_IN(); this->C.initialize(voigt_h::size * voigt_h::size); this->gelstrain.initialize(spatial_dimension * spatial_dimension); this->non_reacted_gel.initialize(1); this->non_reacted_gel.setDefaultValue(1.0); this->damage_ratio.initialize(1); this->damage_ratio_paste.initialize(1); this->damage_ratio_agg.initialize(1); 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"); this->registerParam("k", k, _pat_parsable | _pat_modifiable, "pre-exponential factor of Arrhenius law"); this->registerParam("activ_energy", activ_energy, _pat_parsable | _pat_modifiable, "activation energy of ASR in Arrhenius law"); this->registerParam("R", R, _pat_parsable | _pat_modifiable, "universal gas constant R in Arrhenius law"); this->registerParam("saturation_constant", sat_const, Real(0.0), _pat_parsable | _pat_modifiable, "saturation constant"); this->registerParam("reset_damage", reset_damage, false, _pat_parsable | _pat_modifiable, "reset damage"); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// create a Mesh and SolidMechanicsModel on each integration point of the /// material auto & mesh = this->model.getMesh(); auto & g_ids = mesh.template getData("global_ids", this->el_type); // AKANTU_DEBUG_ASSERT(g_ids.size() > 0, "Global numbering array is empty"); auto const & element_filter = this->getElementFilter()(this->el_type); for (auto && data : zip(arange(element_filter.size()), element_filter, make_view(C(this->el_type), voigt_h::size, voigt_h::size))) { UInt mat_el_id = std::get<0>(data); UInt proc_el_id = std::get<1>(data); UInt gl_el_id = g_ids(proc_el_id); auto & C = std::get<2>(data); meshes.emplace_back(std::make_unique( spatial_dimension, "RVE_mesh_" + std::to_string(gl_el_id), mat_el_id + 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(gl_el_id), mat_el_id + 1)); auto & RVE = *RVEs.back(); RVE.initFull(_analysis_method = _static); /// compute intial stiffness of the RVE RVE.homogenizeStiffness(C, RVE.isTensileHomogen()); } 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::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Compute thermal stresses first Parent::computeStress(el_type, ghost_type); - /// Iterate between macro- and meso-scales for most of the cases - if (not use_homogenized_stiffness) { - for (auto && data : - zip(RVEs, - make_view(this->gradu(el_type), spatial_dimension, - spatial_dimension), - make_view(this->stress(el_type), spatial_dimension, - spatial_dimension), - this->sigma_th(el_type), - make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), - make_view(this->gelstrain(this->el_type), spatial_dimension, - spatial_dimension), - this->delta_T(this->el_type))) { - auto & RVE = *(std::get<0>(data)); - - /// reset nodal and internal fields - RVE.resetNodalFields(); - RVE.resetInternalFields(); - - /// reset the damage to the previously converged state - if (reset_damage) - RVE.restoreDamageField(); - - /// apply boundary conditions based on the current macroscopic displ. - /// gradient - RVE.applyBoundaryConditionsRve(std::get<1>(data)); - - // /// apply temperature only for the output - // RVE.applyHomogeneousTemperature(std::get<7>(data)); - - /// advance the ASR in every RVE based on the new gel strain - RVE.advanceASR(std::get<5>(data)); - - /// compute the average average rVE stress - RVE.homogenizeStressField(std::get<2>(data)); - - // /// compute the new effective stiffness of the RVE - // if (not reset_damage) { - // /// decide whether stiffness homogenization is done via tension - // RVE.setStiffHomogenDir(std::get<2>(data)); - // /// compute the new effective stiffness of the RVE - // RVE.homogenizeStiffness(std::get<4>(data), RVE.isTensileHomogen()); - // } - } - } - /// use homogen stiffness to solve macro-problem (for residual check) - else { - 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); + for (auto && data : + zip(RVEs, + make_view(this->gradu(el_type), spatial_dimension, + spatial_dimension), + make_view(this->stress(el_type), spatial_dimension, + spatial_dimension), + this->sigma_th(el_type), + make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), + make_view(this->gelstrain(this->el_type), spatial_dimension, + spatial_dimension), + this->delta_T(this->el_type))) { + auto & RVE = *(std::get<0>(data)); - const Matrix & C_mat = *C_it; - const Real & sigma_th = *sigma_th_it; + /// reset nodal and internal fields + RVE.resetNodalFields(); + RVE.resetInternalFields(); - /// 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]; + /// reset the damage to the previously converged state + if (reset_damage) + RVE.restoreDamageField(); - voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; - } + /// apply boundary conditions based on the current macroscopic displ. + /// gradient + RVE.applyBoundaryConditionsRve(std::get<1>(data)); - // compute stresses in Voigt notation - voigt_stress.mul(C_mat, voigt_strain); + // /// apply temperature only for the output + // RVE.applyHomogeneousTemperature(std::get<7>(data)); - /// 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; - } + /// advance the ASR in every RVE based on the new gel strain + RVE.advanceASR(std::get<5>(data)); - ++C_it; - ++sigma_th_it; + /// compute the average average rVE stress + RVE.homogenizeStressField(std::get<2>(data)); - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + // /// compute the new effective stiffness of the RVE + // if (not reset_damage) { + // /// decide whether stiffness homogenization is done via tension + // RVE.setStiffHomogenDir(std::get<2>(data)); + // /// compute the new effective stiffness of the RVE + // RVE.homogenizeStiffness(std::get<4>(data), RVE.isTensileHomogen()); + // } } - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::increaseGelStrain(Real & dt_day) { for (auto && data : zip(this->delta_T(this->el_type), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension), this->non_reacted_gel(this->el_type))) { /// compute new gel strain for every element if (this->sat_const) computeNewGelStrainTimeDependent(std::get<1>(data), dt_day, std::get<0>(data), std::get<2>(data)); else computeNewGelStrain(std::get<1>(data), dt_day, std::get<0>(data)); } } /* -------------------------------------------------------------------------- */ template void MaterialFE2::beforeSolveStep() { AKANTU_DEBUG_IN(); - for (const auto & type : - this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { - Array & elem_filter = this->element_filter(type, _not_ghost); + // for (const auto & type : + // this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { + // Array & elem_filter = this->element_filter(type, _not_ghost); - if (elem_filter.size() == 0) - return; - } + // if (elem_filter.size() == 0) + // return; + // } // for (auto && data : // zip(RVEs, // make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), // make_view(this->stress(el_type), spatial_dimension, // spatial_dimension))) { // auto & RVE = *(std::get<0>(data)); // if (reset_damage) // RVE.storeDamageField(); // if (reset_damage) { // /// decide whether stiffness homogenization is done via tension // RVE.setStiffHomogenDir(std::get<2>(data)); // /// compute the new effective stiffness of the RVE // RVE.homogenizeStiffness(std::get<1>(data), RVE.isTensileHomogen()); // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::afterSolveStep() { AKANTU_DEBUG_IN(); for (const auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { Array & elem_filter = this->element_filter(type, _not_ghost); if (elem_filter.size() == 0) return; } for (auto && data : zip(RVEs, make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), make_view(this->stress(el_type), spatial_dimension, spatial_dimension), this->delta_T(this->el_type), this->damage_ratio(this->el_type), this->damage_ratio_paste(this->el_type), this->damage_ratio_agg(this->el_type))) { auto & RVE = *(std::get<0>(data)); if (reset_damage) RVE.storeDamageField(); // if (reset_damage) { /// decide whether stiffness homogenization is done via tension RVE.setStiffHomogenDir(std::get<2>(data)); /// compute the new effective stiffness of the RVE RVE.homogenizeStiffness(std::get<1>(data), RVE.isTensileHomogen()); // } /// compute damage ratio in each RVE RVE.computeDamageRatio(std::get<4>(data)); /// compute damage ratio per material RVE.computeDamageRatioPerMaterial(std::get<5>(data), "paste"); RVE.computeDamageRatioPerMaterial(std::get<6>(data), "aggregate"); // /// apply temperature only for the output // RVE.applyHomogeneousTemperature(std::get<1>(data)); } 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), this->delta_T(this->el_type), this->damage_ratio(this->el_type))) { auto & RVE = *(std::get<0>(data)); /// apply boundary conditions based on the current macroscopic displ. /// gradient RVE.applyBoundaryConditionsRve(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 damage volume in each rve RVE.computeDamageRatio(std::get<5>(data)); /// 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), RVE.isTensileHomogen()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::advanceASR(const Real & delta_time) { 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), this->delta_T(this->el_type), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension), this->non_reacted_gel(this->el_type), this->damage_ratio(this->el_type))) { auto & RVE = *(std::get<0>(data)); /// apply boundary conditions based on the current macroscopic displ. /// gradient RVE.applyBoundaryConditionsRve(std::get<1>(data)); /// apply homogeneous temperature field to each RVE to obtain /// thermoelastic effect RVE.applyHomogeneousTemperature(std::get<4>(data)); /// compute new gel strain for every element if (this->sat_const) computeNewGelStrainTimeDependent(std::get<5>(data), delta_time, std::get<4>(data), std::get<6>(data)); else computeNewGelStrain(std::get<5>(data), delta_time, std::get<4>(data)); /// advance the ASR in every RVE based on the new gel strain RVE.advanceASR(std::get<5>(data)); /// compute damage volume in each rve RVE.computeDamageRatio(std::get<7>(data)); /// remove temperature field - not to mess up with the stiffness /// homogenization further RVE.removeTemperature(); /// 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), RVE.isTensileHomogen()); /// apply temperature back for the output RVE.applyHomogeneousTemperature(std::get<4>(data)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::computeNewGelStrain( Matrix & gelstrain, const Real & delta_time_day, const Real & T) { AKANTU_DEBUG_IN(); const auto & k = this->k; const auto & Ea = this->activ_energy; const auto & R = this->R; /// compute increase in gel strain value for interval of time delta_time /// as temperatures are stored in C, conversion to K is done Real delta_strain = k * std::exp(-Ea / (R * (T + 273.15))) * delta_time_day; for (UInt i = 0; i != spatial_dimension; ++i) gelstrain(i, i) += delta_strain; AKANTU_DEBUG_OUT(); } /* --------------------------------------------------------------------*/ template void MaterialFE2::computeNewGelStrainTimeDependent( Matrix & gelstrain, const Real & delta_time_day, const Real & T, Real & non_reacted_gel) { AKANTU_DEBUG_IN(); const auto & k = this->k; const auto & Ea = this->activ_energy; const auto & R = this->R; const auto & sat_const = this->sat_const; /// compute increase in gel strain value for interval of time delta_time /// as temperatures are stored in C, conversion to K is done Real delta_strain = non_reacted_gel * k * std::exp(-Ea / (R * (T + 273.15))) * delta_time_day; for (UInt i = 0; i != spatial_dimension; ++i) gelstrain(i, i) += delta_strain; non_reacted_gel -= std::exp(-Ea / (R * (T + 273.15))) * delta_time_day / sat_const; if (non_reacted_gel < 0.) non_reacted_gel = 0.; AKANTU_DEBUG_OUT(); } /* --------------------------------------------------------------------*/ template void MaterialFE2::resetGelStrain( const Real & old_time_step) { AKANTU_DEBUG_IN(); const auto & k = this->k; const auto & Ea = this->activ_energy; const auto & R = this->R; const auto & sat_const = this->sat_const; for (auto && data : zip(this->delta_T(this->el_type), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension), this->non_reacted_gel(this->el_type))) { auto & gelstrain = std::get<1>(data); auto & T = std::get<0>(data); auto & non_reacted_gel = std::get<2>(data); non_reacted_gel += std::exp(-Ea / (R * (T + 273.15))) * old_time_step / sat_const; Real prev_delta_strain = non_reacted_gel * k * std::exp(-Ea / (R * (T + 273.15))) * old_time_step; for (UInt i = 0; i != spatial_dimension; ++i) gelstrain(i, i) += -prev_delta_strain; if (non_reacted_gel > 1.) non_reacted_gel = 1.; } AKANTU_DEBUG_OUT(); } /* ----------------------------------------------------------- */ template Real MaterialFE2::computeAverageGelStrain() { AKANTU_DEBUG_IN(); Real av_gelstrain = 0; UInt nb_RVEs = 0; for (auto && data : enumerate(make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension))) { av_gelstrain += std::get<1>(data)(0, 0); nb_RVEs = std::get<0>(data) + 1; } auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(av_gelstrain, SynchronizerOperation::_sum); comm.allReduce(nb_RVEs, SynchronizerOperation::_sum); return av_gelstrain /= nb_RVEs; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::setDirectoryToRveDumper( const std::string & directory) { AKANTU_DEBUG_IN(); for (auto && data : RVEs) { /// set default directory to all RVEs data->setDirectory(directory); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialFE2::getNbRVEs() { AKANTU_DEBUG_IN(); UInt nb_RVEs = 0; for (auto && data : enumerate(RVEs)) { nb_RVEs = std::get<0>(data) + 1; } return nb_RVEs; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::dump() { AKANTU_DEBUG_IN(); for (auto && RVE : RVEs) { /// update stress field before dumping RVE->assembleInternalForces(); /// dump all the RVEs RVE->dumpRve(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::setTimeStep(Real time_step) { AKANTU_DEBUG_IN(); for (auto && RVE : RVEs) { /// set time step to all the RVEs RVE->setTimeStep(time_step); } AKANTU_DEBUG_OUT(); } INSTANTIATE_MATERIAL(material_FE2, MaterialFE2); } // namespace akantu diff --git a/extra_packages/asr-tools/src/material_FE2/material_FE2.hh b/extra_packages/asr-tools/src/material_FE2/material_FE2.hh index bd3e0563a..d74076298 100644 --- a/extra_packages/asr-tools/src/material_FE2/material_FE2.hh +++ b/extra_packages/asr-tools/src/material_FE2/material_FE2.hh @@ -1,190 +1,181 @@ /** * @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 */ template class MaterialFE2 : public MaterialThermal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: typedef MaterialThermal Parent; 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); /// gel strain inrease linear with time, exponential with temperature void computeNewGelStrain(Matrix & gelstrain, const Real & delta_time_day, const Real & temp); /// assymptotic gel strain - time curve void computeNewGelStrainTimeDependent(Matrix & gelstrain, const Real & delta_time_day, const Real & T, Real & non_reacted_gel); /// reset the gel strain value to a previous value void resetGelStrain(const Real & old_time_step); /// advance alkali-silica reaction by the user-provided gel strain void advanceASR(const Matrix & prestrain); /// advance alkali-silica reaction based on delta time and temperature- /// dependent reaction rate void advanceASR(const Real & delta_time); /// compute amount of gel strain averaged across all RVEs Real computeAverageGelStrain(); /// set default dumper directory to all rves void setDirectoryToRveDumper(const std::string & directory); /// compute number of RVEs UInt getNbRVEs(); /// dump all the rves void dump(); /// increase gel strain according to time step void increaseGelStrain(Real & dt); /// set time step to all RVEs void setTimeStep(Real dt); /// Store previously converged damage values in damage_stored virtual void beforeSolveStep(); /// update damage ratio after converged step virtual void afterSolveStep(); - /// enable usage of homogenized stiffness at the macro-scale - void enableHomogenStiffness() { this->use_homogenized_stiffness = true; }; - - /// disable usage of homogenized stiffness at the macro-scale - void disableHomogenStiffness() { this->use_homogenized_stiffness = false; }; - 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; /// pre-exponential factor of Arrhenius law Real k; /// activation energy of ASR in Arrhenius law Real activ_energy; /// universal gas constant; Real R; /// saturation constant for time dependent gel strain increase Real sat_const; /// current gelstrain due to ASR at each Gauss point InternalField gelstrain; /// percent of yet non-reacted gel (for time-dependent asr simulation) InternalField non_reacted_gel; /// ratio between area of damaged elements weighted by damage value /// and the total area of RVE InternalField damage_ratio; InternalField damage_ratio_paste; InternalField damage_ratio_agg; - /// flag for using the homogenized stiffness to solve macro-scale problem - bool use_homogenized_stiffness; /// TODO remove it and associated functions - /// flag to reset damage to previously converged values on each iteration - bool reset_damage; + bool reset_damage{false}; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_FE2_inline_impl.cc" } // namespace akantu #endif /* __AKANTU_MATERIAL_FE_2_HH__ */