diff --git a/extra_packages/asr-tools/package.cmake b/extra_packages/asr-tools/package.cmake index 44b20b5f9..3be8b4320 100644 --- a/extra_packages/asr-tools/package.cmake +++ b/extra_packages/asr-tools/package.cmake @@ -1,27 +1,36 @@ #=============================================================================== # @file package.cmake # # @author Nicolas Richart # # @brief package description for asr stuff # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # #=============================================================================== package_declare(asr_tools DESCRIPTION "ASR stuffs materials, model FE2, toolboxes" DEPENDS extra_materials) package_declare_sources(asr_tools asr_tools.cc asr_tools.hh material_FE2/material_FE2.hh material_FE2/material_FE2.cc material_FE2/material_FE2_inline_impl.cc solid_mechanics_model_RVE.hh solid_mechanics_model_RVE.cc + mesh_utils/nodes_pos_updater.hh + mesh_utils/nodes_pos_updater.cc + mesh_utils/nodes_pos_updater_inline_impl.cc + mesh_utils/nodes_flag_updater.hh + mesh_utils/nodes_flag_updater.cc + mesh_utils/nodes_flag_updater_inline_impl.cc + mesh_utils/crack_numbers_updater.hh + mesh_utils/crack_numbers_updater.cc + mesh_utils/crack_numbers_updater_inline_impl.cc ) diff --git a/extra_packages/asr-tools/src/asr_tools.cc b/extra_packages/asr-tools/src/asr_tools.cc index 80f2de069..ec5aaf6fb 100644 --- a/extra_packages/asr-tools/src/asr_tools.cc +++ b/extra_packages/asr-tools/src/asr_tools.cc @@ -1,2752 +1,3084 @@ /** * @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 +#include #include /* -------------------------------------------------------------------------- */ #include "aka_voigthelper.hh" #include "asr_tools.hh" #include "communicator.hh" +#include "crack_numbers_updater.hh" #include "material_FE2.hh" #include "material_damage_iterative_orthotropic.hh" #include "material_iterative_stiffness_reduction.hh" +#include "node_synchronizer.hh" +#include "nodes_flag_updater.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" #include "solid_mechanics_model_RVE.hh" +#include "solid_mechanics_model_cohesive.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ASRTools::ASRTools(SolidMechanicsModel & model) - : model(model), volume(0.), nb_dumps(0), - cohesive_insertion(false), + : model(model), volume(0.), nb_dumps(0), cohesive_insertion(false), doubled_facets_ready(false), doubled_nodes_ready(false), disp_stored(0, model.getSpatialDimension()), ext_force_stored(0, model.getSpatialDimension()), boun_stored(0, model.getSpatialDimension()), - tensile_homogenization(false) { + tensile_homogenization(false), modified_pos(false), border_nodes(false), + ASR_nodes(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); /// find four corner nodes of the RVE findCornerNodes(); + + /// initialize the modified_pos array with the initial number of nodes + auto nb_nodes = mesh.getNbNodes(); + modified_pos.resize(nb_nodes); + modified_pos.set(false); + border_nodes.resize(nb_nodes); + border_nodes.set(false); + ASR_nodes.resize(nb_nodes); + ASR_nodes.set(false); } /* -------------------------------------------------------------------------- */ void ASRTools::computePhaseVolumes() { /// compute volume of each phase and save it into a map for (auto && mat : model.getMaterials()) { this->phase_volumes[mat.getName()] = computePhaseVolume(mat.getName()); auto it = this->phase_volumes.find(mat.getName()); if (it == this->phase_volumes.end()) this->phase_volumes.erase(mat.getName()); } } /* -------------------------------------------------------------------------- */ 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); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(this->volume, SynchronizerOperation::_sum); } } /* ------------------------------------------------------------------------- */ 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!!!"); /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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_regular)) { 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; } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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"); /// save nodal values before test storeNodalFields(); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); 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::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(); /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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"); 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, damage_paste, crack_agg, crack_paste; computeDamageRatioPerMaterial(damage_agg, "aggregate"); computeDamageRatioPerMaterial(damage_paste, "paste"); computeCrackVolumePerMaterial(crack_agg, "aggregate"); computeCrackVolumePerMaterial(crack_paste, "paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << damage_agg << "," << damage_paste << "," << crack_agg << "," << crack_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 << "," << crack_agg << "," << crack_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"); 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, damage_paste, crack_agg, crack_paste; computeDamageRatioPerMaterial(damage_agg, "aggregate"); computeDamageRatioPerMaterial(damage_paste, "paste"); computeCrackVolumePerMaterial(crack_agg, "aggregate"); computeCrackVolumePerMaterial(crack_paste, "paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); 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 << "," << crack_agg << "," << crack_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 << "," << crack_agg << "," << crack_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAveragePropertiesCohesiveModel( 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"); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << 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 << 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"); 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 crack_agg = averageScalarField("crack_volume_ratio_agg"); Real crack_paste = averageScalarField("crack_volume_ratio_paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << crack_agg << "," << crack_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 << "," << crack_agg << "," << crack_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); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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(); // } /* -------------------------------------------------------------------------- */ Real ASRTools::computeSmallestElementSize() { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // constexpr auto dim = model.getSpatialDimension(); /// compute smallest element size const Array & pos = mesh.getNodes(); Real el_h_min = std::numeric_limits::max(); 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 /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { 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; // } // } // } // } /* -------------------------------------------------------------------------- */ void ASRTools::computeASRStrainLarive( const Real & delta_time_day, const Real & T, Real & ASRStrain, const Real & eps_inf, const Real & time_ch_ref, const Real & time_lat_ref, const Real & U_C, const Real & U_L, const Real & T_ref) { AKANTU_DEBUG_IN(); Real time_ch, time_lat, lambda, ksi, exp_ref; ksi = ASRStrain / eps_inf; if (T == 0) { ksi += 0; } else { time_ch = time_ch_ref * std::exp(U_C * (1. / T - 1. / T_ref)); time_lat = time_lat_ref * std::exp(U_L * (1. / T - 1. / T_ref)); exp_ref = std::exp(-time_lat / time_ch); lambda = (1 + exp_ref) / (ksi + exp_ref); ksi += delta_time_day / time_ch * (1 - ksi) / lambda; } ASRStrain = ksi * eps_inf; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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)) * delta_time_day; amount_reactive_particles -= std::exp(-activ_energy / (R * T)) * 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 + // node 1 - left bottom corner if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } - // node 2 + // node 2 - right bottom corner 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 + // node 3 - right top 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 + // node 4 - left top 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) /// store and clear the eigenstrain ///(to exclude stresses due to internal pressure) Array stored_eig(0, dim * dim, 0); storeASREigenStrain(stored_eig); clearASREigenStrain(); /// 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); } } /// 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 and the ASR eigenstrain restoreNodalFields(); restoreASREigenStrain(stored_eig); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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); 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); } } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(damage_ratio, SynchronizerOperation::_sum); } /// 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.; // 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); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(damage_ratio, SynchronizerOperation::_sum); } if (not this->phase_volumes.size()) computePhaseVolumes(); damage_ratio /= this->phase_volumes[material_name]; } /* -------------------------------------------------------------------------- */ void ASRTools::computeCrackVolume(Real & crack_volume_ratio) { crack_volume_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; auto extra_volume_copy = mat.getInternal("extra_volume")(element_type); crack_volume_ratio += Math::reduce(extra_volume_copy); } } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(crack_volume_ratio, SynchronizerOperation::_sum); } /// compute total model volume if (!this->volume) computeModelVolume(); crack_volume_ratio /= this->volume; } /* -------------------------------------------------------------------------- */ void ASRTools::computeCrackVolumePerMaterial(Real & crack_volume, 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(); crack_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 extra_volume_copy = mat.getInternal("extra_volume")(element_type); crack_volume += Math::reduce(extra_volume_copy); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(crack_volume, SynchronizerOperation::_sum); } if (not this->phase_volumes.size()) computePhaseVolumes(); crack_volume /= this->phase_volumes[material_name]; } /* ------------------------------------------------------------------------ */ void ASRTools::dumpRve() { // if (this->nb_dumps % 10 == 0) { model.dump(); // } // this->nb_dumps += 1; } /* ------------------------------------------------------------------------ */ // void ASRTools::applyBodyForce(const Real gravity = 9.81) { // auto spatial_dimension = model.getSpatialDimension(); // model.assembleMassLumped(); // auto & mass = model.getMass(); // auto & force = model.getExternalForce(); // Vector gravity(spatial_dimension); // gravity(1) = -gravity; // for (auto && data : zip(make_view(mass, spatial_dimension), // make_view(force, spatial_dimension))) { // const auto & mass_vec = (std::get<0>(data)); // AKANTU_DEBUG_ASSERT(mass_vec.norm(), "Mass vector components are zero"); // auto & force_vec = (std::get<1>(data)); // force_vec += gravity * mass_vec; // } // } /* ------------------------------------------------------------------------ */ void ASRTools::insertCohElemByCoords(const Matrix & positions) { 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); - auto & doubled_facets = - mesh_facets.createElementGroup("doubled_facets", dim - 1); - auto & doubled_nodes = mesh.createNodeGroup("doubled_nodes"); for (auto & position : positions) { 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); // Standard mersenne_twister_engine seeded with 0 UInt nb_element = matrix_elements.getElements(type_facets).size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_element - 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(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertFacetsRandomly(const UInt & nb_insertions, std::string facet_mat_name, - bool add_neighbors, - bool only_double_facets) { + Real gap_ratio) { AKANTU_DEBUG_IN(); + // fill in the border_nodes array + communicateFlagsOnNodes(); + /// pick central facets (and neighbors if needed) - pickFacetsRandomly(nb_insertions, facet_mat_name, add_neighbors); + pickFacetsRandomly(nb_insertions, facet_mat_name); + + insertOppositeFacets(); + + assignCrackNumbers(); - insertOppositeFacets(only_double_facets); + insertGap(gap_ratio); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertFacetsByCoords(const Matrix & positions, - bool add_neighbors, - bool only_double_facets) { + Real gap_ratio) { AKANTU_DEBUG_IN(); + // fill in the border_nodes array + communicateFlagsOnNodes(); + /// pick the facets to duplicate - pickFacetsByCoord(positions, add_neighbors); + pickFacetsByCoord(positions); + + insertOppositeFacets(); - insertOppositeFacets(only_double_facets); + assignCrackNumbers(); + + insertGap(gap_ratio); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } +/* ------------------------------------------------------------------------- + */ +void ASRTools::communicateFlagsOnNodes() { + auto & mesh = model.getMesh(); + auto & synch = mesh.getElementSynchronizer(); + NodesFlagUpdater nodes_flag_updater(mesh, synch, this->border_nodes); + nodes_flag_updater.fillPreventInsertion(); +} +/* ------------------------------------------------------------------------- + */ +void ASRTools::communicateCrackNumbers() { + AKANTU_DEBUG_IN(); + + if (model.getMesh().getCommunicator().getNbProc() == 1) + return; + + auto & coh_model = dynamic_cast(model); + const auto & coh_synch = coh_model.getCohesiveSynchronizer(); + CrackNumbersUpdater crack_numbers_updater(model, coh_synch); + crack_numbers_updater.communicateCrackNumbers(); + AKANTU_DEBUG_OUT(); +} /* ----------------------------------------------------------------------- */ -void ASRTools::pickFacetsByCoord(const Matrix & positions, - bool add_neighbors) { +void ASRTools::pickFacetsByCoord(const Matrix & positions) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(type); auto & doubled_facets = mesh_facets.createElementGroup("doubled_facets", dim - 1); auto & doubled_nodes = mesh.createNodeGroup("doubled_nodes"); auto & facet_conn = mesh_facets.getConnectivity(facet_type, gt); + auto & check_facets = inserter.getCheckFacets(facet_type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto facet_conn_it = make_view(facet_conn, nb_nodes_facet).begin(); + auto && comm = akantu::Communicator::getWorldCommunicator(); + auto prank = comm.whoAmI(); Vector bary_facet(dim); for (auto & position : positions) { Real min_dist = std::numeric_limits::max(); Element cent_facet; cent_facet.ghost_type = gt; 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; + if (check_facets(facet.element)) { + 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, _ghost_type = _not_ghost); - /// eliminate posibility of a ghost neighbor - bool ghost_neighbors{false}; + /// communicate between processors for the element closest to the position + auto local_min_dist = min_dist; + auto && comm = akantu::Communicator::getWorldCommunicator(); + comm.allReduce(min_dist, SynchronizerOperation::_min); + + /// let other processor insert the cohesive at this position + if (local_min_dist != min_dist) + continue; + + /// eliminate possibility of inserting on a partition border + bool border_facets{false}; Vector facet_nodes = facet_conn_it[cent_facet.element]; for (auto node : facet_nodes) { - if (not mesh.isLocalNode(node)) - ghost_neighbors = true; + if (this->border_nodes(node)) + border_facets = true; } - if (ghost_neighbors) { + if (border_facets) { std::cout << "Facet at the position " << position(0) << "," << position(1) - << " is connected to a ghost element. Skipping this position" + << " is close to the partition border. Skipping this position" << std::endl; continue; } - doubled_facets.add(cent_facet); - - // add all facet nodes to the group - for (auto node : arange(nb_nodes_facet)) - doubled_nodes.add(facet_conn(cent_facet.element, node)); - - if (add_neighbors) - pickFacetNeighbors(cent_facet); + /// eliminate possibility of intersecting existing ASR element + bool intersection{false}; + for (auto node : facet_nodes) { + if (this->ASR_nodes(node)) + intersection = true; + } + if (intersection) { + std::cout + << "Facet at the position " << position(0) << "," << position(1) + << " is intersecting another ASR element. Skipping this position" + << std::endl; + continue; + } else { + auto success_with_neighbors = pickFacetNeighbors(cent_facet); + + if (success_with_neighbors) { + doubled_facets.add(cent_facet); + /// add all facet nodes to the group + for (auto node : arange(nb_nodes_facet)) { + doubled_nodes.add(facet_conn(cent_facet.element, node)); + this->ASR_nodes(facet_conn(cent_facet.element, node)) = true; + } + std::cout << "Proc " << prank << " placed 1 ASR site" << std::endl; + continue; + } else + continue; + } } } /* ----------------------------------------------------------------------- */ void ASRTools::pickFacetsRandomly(UInt nb_insertions, - std::string facet_mat_name, - bool add_neighbors) { + std::string facet_mat_name) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(type); auto & doubled_facets = mesh_facets.createElementGroup("doubled_facets", dim - 1); auto & doubled_nodes = mesh.createNodeGroup("doubled_nodes"); auto & matrix_elements = mesh_facets.createElementGroup("matrix_facets", dim - 1); auto facet_material_id = model.getMaterialIndex(facet_mat_name); auto & facet_conn = mesh_facets.getConnectivity(facet_type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto facet_conn_it = make_view(facet_conn, nb_nodes_facet).begin(); + auto & check_facets = inserter.getCheckFacets(facet_type, gt); + auto && comm = akantu::Communicator::getWorldCommunicator(); + auto prank = comm.whoAmI(); if (not nb_insertions) return; Vector bary_facet(dim); 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 == facet_material_id) - matrix_elements.add(facet); + if (check_facets(facet.element)) { + mesh_facets.getBarycenter(facet, bary_facet); + auto & facet_material = coh_model.getFacetMaterial( + facet.type, facet.ghost_type)(facet.element); + if (facet_material == facet_material_id) + matrix_elements.add(facet); + } }, _spatial_dimension = dim - 1, _ghost_type = _not_ghost); UInt nb_element = matrix_elements.getElements(facet_type).size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_element - 1); if (not nb_element) { - auto && comm = akantu::Communicator::getWorldCommunicator(); - auto prank = comm.whoAmI(); std::cout << "Proc " << prank << " couldn't place " << nb_insertions << " ASR sites" << std::endl; return; } - for (UInt i = 0; i < nb_insertions; i++) { + UInt already_inserted = 0; + while (already_inserted < nb_insertions) { auto id = dis(random_generator); Element cent_facet; cent_facet.type = facet_type; cent_facet.element = matrix_elements.getElements(facet_type)(id); cent_facet.ghost_type = gt; - /// eliminate posibility of a ghost neighbor - bool ghost_neighbors{false}; + /// eliminate possibility of inserting on a partition border or intersecting + bool border_facets{false}; Vector facet_nodes = facet_conn_it[cent_facet.element]; for (auto node : facet_nodes) { - if (not mesh.isLocalNode(node)) - ghost_neighbors = true; + if (this->border_nodes(node) or this->ASR_nodes(node)) + border_facets = true; } - if (ghost_neighbors) { - i--; + if (border_facets) continue; + else { + auto success_with_neighbors = pickFacetNeighbors(cent_facet); + + if (success_with_neighbors) { + already_inserted++; + doubled_facets.add(cent_facet); + /// add all facet nodes to the group + for (auto node : arange(nb_nodes_facet)) { + doubled_nodes.add(facet_conn(cent_facet.element, node)); + this->ASR_nodes(facet_conn(cent_facet.element, node)) = true; + } + std::cout << "Proc " << prank << " placed 1 ASR site" << std::endl; + continue; + } else + continue; } - - doubled_facets.add(cent_facet); - - /// add all facet nodes to the group - for (auto node : arange(nb_nodes_facet)) { - doubled_nodes.add(facet_conn(cent_facet.element, node)); - } - - if (add_neighbors) - pickFacetNeighbors(cent_facet); } } /* ----------------------------------------------------------------------- */ -void ASRTools::pickFacetNeighbors(Element & cent_facet) { +bool ASRTools::pickFacetNeighbors(Element & cent_facet) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const auto & pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); auto & doubled_facets = mesh_facets.getElementGroup("doubled_facets"); auto & facet_conn = mesh_facets.getConnectivity(cent_facet.type, cent_facet.ghost_type); + auto & cent_facet_material = coh_model.getFacetMaterial( + cent_facet.type, cent_facet.ghost_type)(cent_facet.element); CSR nodes_to_segments; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_segments, dim - 1); + Array two_neighbors(2); + Element dummy{_not_defined, UInt(-1), _casper}; + two_neighbors.set(dummy); for (auto node : arange(2)) { /// 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; + // dot product less than -0.5 discards any < 120 deg + Real min_dot = -0.5; + // Real min_dot = std::numeric_limits::max(); Vector neighbor_facet_dir(dim); for (auto & elem : nodes_to_segments.getRow(facet_conn(cent_facet.element, node))) { if (elem.element == cent_facet.element) continue; if (elem.type != cent_facet.type) continue; if (elem.ghost_type != cent_facet.ghost_type) continue; + if (not inserter.getCheckFacets(elem.type, elem.ghost_type)(elem.element)) + continue; + // discard neighbors from different materials + auto & candidate_facet_material = + coh_model.getFacetMaterial(elem.type, elem.ghost_type)(elem.element); + if (candidate_facet_material != cent_facet_material) + continue; /// decide which node of the neighbor is the second UInt first_node{facet_conn(cent_facet.element, node)}; UInt second_node(-1); if (facet_conn(elem.element, 0) == first_node) { second_node = facet_conn(elem.element, 1); } else if (facet_conn(elem.element, 1) == first_node) { second_node = facet_conn(elem.element, 0); } else AKANTU_EXCEPTION( "Neighboring facet" << elem << " with nodes " << facet_conn(elem.element, 0) << " and " << facet_conn(elem.element, 1) << " doesn't have node in common with the central facet " << cent_facet << " with nodes " << facet_conn(cent_facet.element, 0) << " and " << facet_conn(cent_facet.element, 1)); - // /// discard facets laying on the border of the partition - // if (not mesh.isLocalNode(second_node)) { - // std::cout << "Neighbor's node is on the border of the partition" - // << std::endl; - // continue; - // } + /// discard facets intersecting other ASR elements + if (this->ASR_nodes(second_node)) + continue; neighbor_facet_dir = pos_it[second_node]; neighbor_facet_dir -= Vector(pos_it[first_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; + two_neighbors(node) = elem; } } - doubled_facets.add(neighbor); } + + // insert neighbors only if two of them were identified + if (two_neighbors.find(dummy) == UInt(-1)) { + for (auto & neighbor : two_neighbors) { + doubled_facets.add(neighbor); + for (UInt node : arange(2)) { + this->ASR_nodes(facet_conn(neighbor.element, node)) = true; + } + } + return true; + } else + return false; } /* ----------------------------------------------------------------------- */ void ASRTools::preventCohesiveInsertionInNeighbors() { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto el_type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(el_type); auto & doubled_nodes = mesh.getNodeGroup("doubled_nodes"); CSR nodes_to_elements; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_elements, dim - 1); for (auto node : doubled_nodes.getNodes()) { for (auto & elem : nodes_to_elements.getRow(node)) { if (elem.type != facet_type) continue; inserter.getCheckFacets(elem.type, gt)(elem.element) = false; } } } /* ----------------------------------------------------------------------- */ -void ASRTools::insertOppositeFacets(bool only_double_facets) { +void ASRTools::insertOppositeFacets() { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto & crack_facets = mesh.createElementGroup("crack_facets"); const auto & pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); auto & doubled_facets = mesh_facets.getElementGroup("doubled_facets"); /// sort facet numbers in doubled_facets_group to comply with new_elements /// event passed by the inserter doubled_facets.optimize(); /// instruct the inserter which facets to duplicate for (auto & type : doubled_facets.elementTypes(dim - 1)) { const auto element_ids = doubled_facets.getElements(type, gt); /// iterate over facets to duplicate for (auto && el : element_ids) { inserter.getCheckFacets(type, gt)(el) = false; Element new_facet{type, el, gt}; insertion(new_facet) = true; } } - /// duplicate facets (and insert coh els if needed) - inserter.insertElements(only_double_facets); + /// duplicate facets and insert coh els + inserter.insertElements(false); /// add facets connectivity to the mesh and the element group NewElementsEvent new_facets_event; for (auto & type : doubled_facets.elementTypes(dim - 1)) { const auto element_ids = doubled_facets.getElements(type, gt); if (not mesh.getConnectivities().exists(type, gt)) mesh.addConnectivityType(type, gt); auto & facet_conn = mesh.getConnectivity(type, gt); auto & mesh_facet_conn = mesh_facets.getConnectivity(type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto mesh_facet_conn_it = make_view(mesh_facet_conn, nb_nodes_facet).begin(); /// iterate over duplicated facets for (auto && el : element_ids) { Vector facet_nodes = mesh_facet_conn_it[el]; facet_conn.push_back(facet_nodes); Element new_facet{type, facet_conn.size() - 1, gt}; crack_facets.add(new_facet); for (auto & facet_node : facet_nodes) { crack_facets.addNode(facet_node, true); } new_facets_event.getList().push_back(new_facet); } } MeshUtils::fillElementToSubElementsData(mesh); mesh.sendEvent(new_facets_event); /// create an element group with nodes to apply Dirichlet model.getMesh().createElementGroupFromNodeGroup("doubled_nodes", "doubled_nodes", dim - 1); /// update FEEngineBoundary with new elements model.getFEEngineBoundary().initShapeFunctions(_not_ghost); model.getFEEngineBoundary().initShapeFunctions(_ghost); model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_not_ghost); model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_ghost); } /* ----------------------------------------------------------------------- */ +void ASRTools::assignCrackNumbers() { + auto & mesh = model.getMesh(); + auto & mesh_facets = mesh.getMeshFacets(); + auto dim = mesh.getSpatialDimension(); + const GhostType gt = _not_ghost; + + for (auto && type_facet : mesh_facets.elementTypes(dim - 1)) { + ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); + UInt nb_quad_cohesive = model.getFEEngine("CohesiveFEEngine") + .getNbIntegrationPoints(type_cohesive); + auto & coh_conn = mesh.getConnectivity(type_cohesive, gt); + auto nb_coh_nodes = coh_conn.getNbComponent(); + auto coh_conn_it = coh_conn.begin(coh_conn.getNbComponent()); + CSR nodes_to_elements; + MeshUtils::buildNode2Elements(mesh, nodes_to_elements, dim, _ek_cohesive); + + // initialize a graph + typedef boost::adjacency_list + Graph; + Graph graph; + + // build the graph based on the connectivity of cohesive elements + for (UInt i = 0; i < coh_conn.size(); i++) { + for (UInt j = i; j < coh_conn.size(); j++) { + auto nodes_1 = coh_conn_it[i]; + auto nodes_2 = coh_conn_it[j]; + std::vector nod_1(nb_coh_nodes); + std::vector nod_2(nb_coh_nodes); + for (UInt k : arange(nb_coh_nodes)) { + nod_1[k] = nodes_1(k); + nod_2[k] = nodes_2(k); + } + + std::sort(nod_1.begin(), nod_1.end()); + std::sort(nod_2.begin(), nod_2.end()); + std::vector common_nodes(nb_coh_nodes); + std::vector::iterator it; + + it = std::set_intersection(nod_1.begin(), nod_1.end(), nod_2.begin(), + nod_2.end(), common_nodes.begin()); + common_nodes.resize(it - common_nodes.begin()); + + /// if any common nodes are found, add corresponding edge to the graph + if (common_nodes.size()) { + boost::add_edge(i, j, graph); + } + } + } + + /// connectivity and number of components of a graph + std::vector component(boost::num_vertices(graph)); + int num = boost::connected_components(graph, &component[0]); + + /// shift the crack numbers by the number of cracks on processors with the + /// lower rank + auto && comm = akantu::Communicator::getWorldCommunicator(); + comm.exclusiveScan(num); + for (auto & component_nb : component) { + component_nb += num; + } + + /// assign a corresponding crack flag + for (UInt coh_el : arange(component.size())) { + const Array & material_index_vec = + model.getMaterialByElement(type_cohesive, gt); + const Array & material_local_numbering_vec = + model.getMaterialLocalNumbering(type_cohesive, gt); + Material & material = model.getMaterial(material_index_vec(coh_el)); + auto & crack_nb = + material.getInternal("crack_number")(type_cohesive, gt); + UInt material_local_num = material_local_numbering_vec(coh_el); + auto crack_nb_it = crack_nb.begin(); + for (UInt i : arange(nb_quad_cohesive)) { + crack_nb_it[material_local_num * nb_quad_cohesive + i] = + Real(component[coh_el]); + } + } + } +} +/* ------------------------------------------------------------------------ */ +void ASRTools::insertGap(const Real gap_ratio) { + AKANTU_DEBUG_IN(); + if (gap_ratio == 0) + return; + + auto & mesh = model.getMesh(); + MeshAccessor mesh_accessor(mesh); + auto & pos2modify = mesh_accessor.getNodes(); + + auto dim = mesh.getSpatialDimension(); + const GhostType gt = _not_ghost; + auto && group = mesh.getElementGroup("doubled_nodes"); + auto & pos = mesh.getNodes(); + auto pos_it = make_view(pos, dim).begin(); + auto pos2modify_it = make_view(pos2modify, dim).begin(); + + for (auto & type : group.elementTypes(dim - 1)) { + auto & facet_conn = mesh.getConnectivity(type, gt); + const UInt nb_nodes_facet = facet_conn.getNbComponent(); + auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); + const auto element_ids = group.getElements(type, gt); + auto && fe_engine = model.getFEEngineBoundary(); + auto nb_qpoints_per_facet = fe_engine.getNbIntegrationPoints(type, gt); + const auto & normals_on_quad = + fe_engine.getNormalsOnIntegrationPoints(type, gt); + auto normals_it = make_view(normals_on_quad, dim).begin(); + for (UInt i : arange(element_ids.size())) { + auto el_id = element_ids(i); + Vector normal = normals_it[el_id * nb_qpoints_per_facet]; + if (i < element_ids.size() / 2) + normal *= -1; + /// compute segment length + auto facet_nodes = facet_nodes_it[el_id]; + UInt A, B; + A = facet_nodes(0); + B = facet_nodes(1); + Vector AB = Vector(pos_it[B]) - Vector(pos_it[A]); + Real correction = AB.norm() * gap_ratio / 2; + Vector half_opening_vector = normal * correction; + for (auto j : arange(nb_nodes_facet)) { + Vector node_pos(pos2modify_it[facet_nodes(j)]); + node_pos += half_opening_vector; + this->modified_pos(facet_nodes(j)) = true; + } + } + } + + /// update FEEngine & FEEngineBoundary with new elements + model.getFEEngine().initShapeFunctions(_not_ghost); + model.getFEEngine().initShapeFunctions(_ghost); + model.getFEEngineBoundary().initShapeFunctions(_not_ghost); + model.getFEEngineBoundary().initShapeFunctions(_ghost); + model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_not_ghost); + model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_ghost); + + AKANTU_DEBUG_OUT(); +} +/* ----------------------------------------------------------------------- */ void ASRTools::onElementsAdded(const Array & elements, const NewElementsEvent &) { - /// function is activated only when expanding cohesive elements is on + // function is activated only when expanding cohesive elements is on if (not this->cohesive_insertion) return; if (this->doubled_facets_ready) { return; } + auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto & doubled_facets = mesh_facets.getElementGroup("doubled_facets"); 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; /// add new facets into the doubled_facets group auto & element_ids = elements_range.getElements(); for (auto && el : element_ids) { Element new_facet{type, el, ghost_type}; doubled_facets.add(new_facet); } this->doubled_facets_ready = true; } } /* -------------------------------------------------------------------------- */ void ASRTools::onNodesAdded(const Array & new_nodes, const NewNodesEvent &) { AKANTU_DEBUG_IN(); if (new_nodes.size() == 0) return; + /// increase the internal arrays by the number of new_nodes + UInt new_nb_nodes = this->modified_pos.size() + new_nodes.size(); + this->modified_pos.resize(new_nb_nodes); + this->border_nodes.resize(new_nb_nodes); + this->ASR_nodes.resize(new_nb_nodes); + /// function is activated only when expanding cohesive elements is on if (not this->cohesive_insertion) return; if (this->doubled_nodes_ready) return; auto & mesh = model.getMesh(); auto & node_group = mesh.getNodeGroup("doubled_nodes"); auto & central_nodes = node_group.getNodes(); auto pos_it = make_view(mesh.getNodes(), mesh.getSpatialDimension()).begin(); for (auto & new_node : new_nodes) { const Vector & new_node_coord = pos_it[new_node]; for (auto & central_node : central_nodes) { const Vector & central_node_coord = pos_it[central_node]; if (new_node_coord == central_node_coord) node_group.add(new_node); } } this->doubled_nodes_ready = true; AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ void ASRTools::updateElementGroup(const std::string group_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); AKANTU_DEBUG_ASSERT(mesh.elementGroupExists(group_name), "Element group is not registered in the mesh"); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto && group = mesh.getElementGroup(group_name); auto && pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); for (auto & type : group.elementTypes(dim - 1)) { auto & facet_conn = mesh.getConnectivity(type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); AKANTU_DEBUG_ASSERT( type == _segment_2, "Currently update group works only for el type _segment_2"); const auto element_ids = group.getElements(type, gt); for (auto && el_id : element_ids) { const auto connected_els = mesh.getElementToSubelement(type, gt)(el_id); for (auto && connected_el : connected_els) { auto type_solid = connected_el.type; auto & solid_conn = mesh.getConnectivity(type_solid, gt); const UInt nb_nodes_solid_el = solid_conn.getNbComponent(); const auto solid_nodes_it = make_view(solid_conn, nb_nodes_solid_el).begin(); Vector facet_nodes = facet_nodes_it[el_id]; Vector solid_nodes = solid_nodes_it[connected_el.element]; // /// to which of connected elements facet belongs // Array solid_nodes(nb_nodes_solid_el); // for (UInt node : arange(nb_nodes_solid_el)) { // solid_nodes(node) = solid_nodes_it[connected_el.element](node); // } // /// check for the central facet node (id doesn't change nb) // auto id = solid_nodes.find(facet_nodes(2)); // if (id == UInt(-1)) // continue; /// check only the corner nodes of facets - central will not change for (auto f : arange(2)) { auto facet_node = facet_nodes(f); const Vector & facet_node_coords = pos_it[facet_node]; for (auto s : arange(nb_nodes_solid_el)) { auto solid_node = solid_nodes(s); const Vector & solid_node_coords = pos_it[solid_node]; if (solid_node_coords == facet_node_coords) { if (solid_node != facet_node) { // group.removeNode(facet_node); facet_conn(el_id, f) = solid_node; // group.addNode(solid_node, true); } break; } } } } } } group.optimize(); 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::clearASREigenStrain() { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Matrix zero_eigenstrain(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_eigenstrain; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::storeASREigenStrain(Array & stored_eig) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; UInt nb_quads = 0; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const UInt nb_gp = model.getFEEngine().getNbIntegrationPoints(element_type); nb_quads += model.getMaterial("gel").getElementFilter(element_type).size() * nb_gp; } stored_eig.resize(nb_quads); auto stored_eig_it = stored_eig.begin(dim, dim); 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, ++stored_eig_it) (*stored_eig_it) = (*prestrain_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreASREigenStrain(Array & stored_eig) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; auto stored_eig_it = stored_eig.begin(dim, dim); 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, ++stored_eig_it) (*prestrain_it) = (*stored_eig_it); } AKANTU_DEBUG_OUT(); } +/* ------------------------------------------------------------------------ */ +template UInt ASRTools::insertSingleCohesiveElementPerModel() { + AKANTU_DEBUG_IN(); + + auto & coh_model = dynamic_cast(model); + CohesiveElementInserter & inserter = coh_model.getElementInserter(); + auto && comm = akantu::Communicator::getWorldCommunicator(); + auto prank = comm.whoAmI(); + + if (not coh_model.getIsExtrinsic()) { + AKANTU_EXCEPTION( + "This function can only be used for extrinsic cohesive elements"); + } + + coh_model.interpolateStress(); + + Mesh & mesh_facets = coh_model.getMeshFacets(); + MeshUtils::fillElementToSubElementsData(mesh_facets); + UInt nb_new_elements(0); + + for (auto && type_facet : mesh_facets.elementTypes(dim - 1)) { + + Real max_stress = std::numeric_limits::min(); + std::string max_stress_mat_name; + int max_stress_prank(-1); + UInt max_stress_facet(-1); + Real crack_nb(-1); + std::tuple max_stress_data(max_stress_facet, max_stress, + crack_nb); + + for (auto && mat : model.getMaterials()) { + + auto * mat_coh = + dynamic_cast *>(&mat); + + if (mat_coh == nullptr) + continue; + + // check which material has the non ghost facet with the highest stress + // which complies with the topological criteria + auto max_stress_data = mat_coh->findCriticalFacet(type_facet); + auto max_stress_mat = std::get<1>(max_stress_data); + auto max_stress_facet_mat = std::get<0>(max_stress_data); + auto crack_nb_mat = std::get<2>(max_stress_data); + + // communicate between processors the highest stress + Real local_max_stress = max_stress_mat; + comm.allReduce(max_stress_mat, SynchronizerOperation::_max); + + // write down the maximum stress, material_ID and proc # + if (max_stress_mat > max_stress) { + max_stress = max_stress_mat; + max_stress_mat_name = mat.getName(); + max_stress_facet = max_stress_facet_mat; + crack_nb = crack_nb_mat; + if (local_max_stress == max_stress_mat) + max_stress_prank = prank; + } + } + + // if max stress is low or another processor has the most stressed el ->skip + if ((max_stress < 1) or (max_stress_prank != prank)) + continue; + + // otherwise insert cohesive element in this specific material + auto & mat = coh_model.getMaterial(max_stress_mat_name); + auto * mat_coh = + dynamic_cast *>(&mat); + + mat_coh->insertSingleCohesiveElement(type_facet, max_stress_facet, crack_nb, + false); + } + + // insert cohesive elements + nb_new_elements = inserter.insertElements(); + AKANTU_DEBUG_OUT(); + return nb_new_elements; +} + +template UInt ASRTools::insertSingleCohesiveElementPerModel<2>(); +template UInt ASRTools::insertSingleCohesiveElementPerModel<3>(); } // namespace akantu diff --git a/extra_packages/asr-tools/src/asr_tools.hh b/extra_packages/asr-tools/src/asr_tools.hh index 63d84c5ba..502c4408c 100644 --- a/extra_packages/asr-tools/src/asr_tools.hh +++ b/extra_packages/asr-tools/src/asr_tools.hh @@ -1,944 +1,950 @@ /** * @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 /* -------------------------------------------------------------------------- */ +#include "material_cohesive_linear_sequential.hh" #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 volumes of each phase void computePhaseVolumes(); /// 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); /// computes and writes only displacements and strains void computeAveragePropertiesCohesiveModel(std::ofstream & file_output, Real time); /// Averages strains & displacements + damage ratios in FE2 material void computeAveragePropertiesFe2Material(std::ofstream & file_output, Real time); /// 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); Real computeSmallestElementSize(); // /// apply homogeneous temperature on Solid Mechanics Model // void applyTemperatureFieldToSolidmechanicsModel(const Real & temperature); /// compute ASR strain by a sigmoidal rule (Larive, 1998) void computeASRStrainLarive(const Real & delta_time_day, const Real & T, Real & ASRStrain, const Real & eps_inf, const Real & time_ch_ref, const Real & time_lat_ref, const Real & U_C, const Real & U_L, const Real & T_ref); /// 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 coordinates of their center void insertCohElemByCoords(const Matrix & positions); /// insert multiple cohesive elements by the limiting box void insertCohElemByLimits(const Matrix & insertion_limits, std::string coh_mat_name); /// insert multiple cohesive elements randomly void insertCohElemRandomly(const UInt & nb_coh_elem, std::string coh_mat_name, std::string matrix_mat_name); /// insert multiple facets (and cohesive elements if needed) void insertFacetsRandomly(const UInt & nb_coh_elem, - std::string matrix_mat_name, - bool add_neighbors = true, - bool only_double_facets = true); + std::string matrix_mat_name, Real gap_ratio); /// insert up to 3 facets pair based on the coord of the central one - void insertFacetsByCoords(const Matrix & positions, bool add_neighbors, - bool only_double_facets); + void insertFacetsByCoords(const Matrix & positions, Real gap_ratio = 0); + + /// communicates crack numbers from not ghosts to ghosts cohesive elements + void communicateCrackNumbers(); protected: + /// put flags on all nodes who have a ghost counterpart + void communicateFlagsOnNodes(); + /// pick facets by passed coordinates - void pickFacetsByCoord(const Matrix & positions, bool add_neighbors); + void pickFacetsByCoord(const Matrix & positions); /// pick facets by passed coordinates - void pickFacetsRandomly(UInt nb_insertions, std::string facet_mat_name, - bool add_neighbors); + void pickFacetsRandomly(UInt nb_insertions, std::string facet_mat_name); - /// pick two neighbors of a central facet - void pickFacetNeighbors(Element & cent_facet); + /// pick two neighbors of a central facet: returns true if success + bool pickFacetNeighbors(Element & cent_facet); - /// optimise doubled facets group, insert facets, update connectivities - void insertOppositeFacets(bool only_double_facets); + /// optimise doubled facets group, insert facets, and cohesive elements, + /// update connectivities + void insertOppositeFacets(); /// no cohesive elements in-between neighboring solid elements void preventCohesiveInsertionInNeighbors(); + /// assign crack tags to the clusters + void assignCrackNumbers(); + + /// change coordinates of central crack nodes to create an artificial gap + void insertGap(const Real gap_ratio); + /// 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); public: /// update the element group in case connectivity changed after cohesive /// elements insertion void updateElementGroup(const std::string group_name); + /// works only for the MaterialCohesiveLinearSequential + template UInt insertSingleCohesiveElementPerModel(); + // /// 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); /// compute ratio between total crack and RVE volumes void computeCrackVolume(Real & crack_volume_ratio); /// compute crack volume to material volume ratio void computeCrackVolumePerMaterial(Real & crack_volume, 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 clearASREigenStrain(); /// store elemental eigenstrain in an array void storeASREigenStrain(Array & stored_eig); /// restore eigenstrain in ASR sites from previously stored values void restoreASREigenStrain(Array & stored_eig); /* ------------------------------------------------------------------------ */ public: // Accessors bool isTensileHomogen() { return this->tensile_homogenization; }; /// phase volumes Real getPhaseVolume(const std::string & material_name) { if (not this->phase_volumes.size()) computePhaseVolumes(); return this->phase_volumes.find(material_name)->second; }; /// set the value of the insertion flag AKANTU_SET_MACRO(CohesiveInsertion, cohesive_insertion, bool); + /// get the corner nodes + AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); + /* --------------------------------------------------------------------- */ /* 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; /// dump counter UInt nb_dumps; /// flag to activate ASR expansion through cohesive elements bool cohesive_insertion; /// booleans for applying delta u bool doubled_facets_ready; bool doubled_nodes_ready; // 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}; /// phase volumes std::map phase_volumes; + + /// array to store flags on nodes position modification + Array modified_pos; + + /// array to store flags on nodes that are synchronized between processors + Array border_nodes; + + /// array to store flags on nodes where ASR elements are inserted + Array ASR_nodes; }; /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ class GelMaterialSelector : public MeshDataMaterialSelector { public: GelMaterialSelector(SolidMechanicsModel & model, std::string gel_material, const UInt nb_gel_pockets, std::string aggregate_material = "aggregate", bool gel_pairs = false) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), aggregate_material(aggregate_material), gel_pairs(gel_pairs) {} 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); std::set checked_baries; while (nb_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); nb_placed_gel_pockets += 1; } } } is_gel_initialized = true; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; } void initGelPocketPairs() { aggregate_material_id = model.getMaterialIndex(aggregate_material); Mesh & mesh = this->model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); 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); std::set checked_baries; while (nb_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) { auto & sub_to_element = mesh_facets.getSubelementToElement(el.type, el.ghost_type); auto sub_to_el_it = sub_to_element.begin(sub_to_element.getNbComponent()); const Vector & subelements_to_element = sub_to_el_it[el.element]; bool successful_placement{false}; for (auto & subelement : subelements_to_element) { auto && connected_elements = mesh_facets.getElementToSubelement( subelement.type, subelement.ghost_type)(subelement.element); for (auto & connected_element : connected_elements) { if (connected_element.element == el.element) continue; if (MeshDataMaterialSelector::operator()( connected_element) == aggregate_material_id) { gel_pockets.push_back(el); gel_pockets.push_back(connected_element); nb_placed_gel_pockets += 1; successful_placement = true; break; } } if (successful_placement) break; } } } } is_gel_initialized = true; std::cout << nb_placed_gel_pockets << " ASR-pocket pairs placed" << std::endl; } UInt operator()(const Element & elem) { if (not is_gel_initialized) { if (this->gel_pairs) { initGelPocketPairs(); } else { 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) { 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; std::string aggregate_material{"aggregate"}; UInt aggregate_material_id{1}; bool is_gel_initialized{false}; bool gel_pairs{false}; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ using MaterialCohesiveRules = std::map, ID>; class GelMaterialCohesiveRulesSelector : public MaterialSelector { public: GelMaterialCohesiveRulesSelector(SolidMechanicsModelCohesive & model, const MaterialCohesiveRules & rules, std::string gel_material, const UInt nb_gel_pockets, std::string aggregate_material = "aggregate", bool gel_pairs = false) : model(model), mesh_data_id("physical_names"), mesh(model.getMesh()), mesh_facets(model.getMeshFacets()), dim(model.getSpatialDimension()), rules(rules), gel_selector(model, gel_material, nb_gel_pockets, aggregate_material, gel_pairs), default_cohesive(model) {} UInt operator()(const Element & element) { if (mesh_facets.getSpatialDimension(element.type) == (dim - 1)) { const std::vector & element_to_subelement = mesh_facets.getElementToSubelement(element.type, element.ghost_type)( element.element); // Array & facets_check = model.getFacetsCheck(); const Element & el1 = element_to_subelement[0]; const Element & el2 = element_to_subelement[1]; ID id1 = model.getMaterial(gel_selector(el1)).getName(); ID id2 = id1; if (el2 != ElementNull) { id2 = model.getMaterial(gel_selector(el2)).getName(); } auto rit = rules.find(std::make_pair(id1, id2)); if (rit == rules.end()) { rit = rules.find(std::make_pair(id2, id1)); } if (rit != rules.end()) { return model.getMaterialIndex(rit->second); } } if (Mesh::getKind(element.type) == _ek_cohesive) { return default_cohesive(element); } return gel_selector(element); } private: SolidMechanicsModelCohesive & model; ID mesh_data_id; const Mesh & mesh; const Mesh & mesh_facets; UInt dim; MaterialCohesiveRules rules; GelMaterialSelector gel_selector; DefaultMaterialCohesiveSelector default_cohesive; }; /* ------------------------------------------------------------------------ */ /* 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 PressureSimple : public BC::Neumann::NeumannFunctor { public: PressureSimple(SolidMechanicsModel & model, const Real pressure, const std::string group_name) : model(model), pressure(pressure), group_name(group_name) {} inline void operator()(const IntegrationPoint & quad_point, Vector & dual, const Vector & /*coord*/, const Vector & normal) const { // get element types auto && mesh = model.getMesh(); AKANTU_DEBUG_ASSERT(mesh.elementGroupExists(group_name), "Element group is not registered in the mesh"); const GhostType gt = akantu::_not_ghost; const UInt facet_nb = quad_point.element; const ElementType type_facet = quad_point.type; auto && facet_conn = mesh.getConnectivity(type_facet, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); auto & group = mesh.getElementGroup(group_name); Array element_ids = group.getElements(type_facet); AKANTU_DEBUG_ASSERT(element_ids.size(), "Provided group doesn't contain this element type"); auto id = element_ids.find(facet_nb); AKANTU_DEBUG_ASSERT(id != UInt(-1), "Quad point doesn't belong to this element group"); auto normal_corrected = normal; if (id < element_ids.size() / 2) normal_corrected *= -1; else if (id >= element_ids.size()) AKANTU_EXCEPTION("Error in defining side of the cohesive element"); else normal_corrected *= 1; dual = pressure * normal_corrected; } protected: SolidMechanicsModel & model; const Real pressure; const std::string group_name; }; /* ------------------------------------------------------------------------ */ class PressureVolumeDependent : public BC::Neumann::NeumannFunctor { public: - PressureVolumeDependent(SolidMechanicsModel & model, const Real ASR_volume, + PressureVolumeDependent(SolidMechanicsModel & model, + const Real ASR_volume_ratio, const std::string group_name, const Real compressibility) - : model(model), ASR_volume(ASR_volume), group_name(group_name), - compressibility(compressibility) {} + : model(model), ASR_volume_ratio(ASR_volume_ratio), + group_name(group_name), compressibility(compressibility) {} inline void operator()(const IntegrationPoint & quad_point, Vector & dual, const Vector & /*coord*/, const Vector & normal) const { // get element types auto && mesh = model.getMesh(); AKANTU_DEBUG_ASSERT(mesh.elementGroupExists(group_name), "Element group is not registered in the mesh"); auto dim = mesh.getSpatialDimension(); const GhostType gt = akantu::_not_ghost; const UInt facet_nb = quad_point.element; const ElementType type_facet = quad_point.type; auto && facet_conn = mesh.getConnectivity(type_facet, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); auto & group = mesh.getElementGroup(group_name); Array element_ids = group.getElements(type_facet); auto && pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); auto && disp = model.getDisplacement(); const auto disp_it = make_view(disp, dim).begin(); AKANTU_DEBUG_ASSERT(element_ids.size(), "Provided group doesn't contain this element type"); auto id = element_ids.find(facet_nb); AKANTU_DEBUG_ASSERT(id != UInt(-1), "Quad point doesn't belong to this element group"); auto normal_corrected = normal; UInt opposite_facet_nb(-1); if (id < element_ids.size() / 2) { normal_corrected *= -1; opposite_facet_nb = element_ids(id + element_ids.size() / 2); } else if (id >= element_ids.size()) AKANTU_EXCEPTION("Error in defining side of the cohesive element"); else { normal_corrected *= 1; opposite_facet_nb = element_ids(id - element_ids.size() / 2); } /// compute current area of a gap - auto facet_nodes = facet_nodes_it[facet_nb]; Vector first_facet_nodes = facet_nodes_it[facet_nb]; Vector second_facet_nodes = facet_nodes_it[opposite_facet_nb]; /// corners of a quadrangle consequently - /// A-------B + /// A---M---B /// \ | - /// D-----C + /// D--N--C UInt A, B, C, D; A = first_facet_nodes(0); B = first_facet_nodes(1); - for (auto & first_facet_node : first_facet_nodes) { - bool pair_found{false}; - const Vector & first_facet_node_coord = pos_it[first_facet_node]; - for (auto & second_facet_node : second_facet_nodes) { - const Vector & second_facet_node_coord = - pos_it[second_facet_node]; - if (first_facet_node_coord == second_facet_node_coord) { - if (first_facet_node == A) - D = second_facet_node; - else - C = second_facet_node; - pair_found = true; - } - } - AKANTU_DEBUG_ASSERT(pair_found, "Node pair was not found"); - } + C = second_facet_nodes(0); + D = second_facet_nodes(1); /// quadrangle's area through diagonals Vector AC, BD; Vector A_pos = Vector(pos_it[A]) + Vector(disp_it[A]); Vector B_pos = Vector(pos_it[B]) + Vector(disp_it[B]); Vector C_pos = Vector(pos_it[C]) + Vector(disp_it[C]); Vector D_pos = Vector(pos_it[D]) + Vector(disp_it[D]); - - AC = C_pos - A_pos; - BD = D_pos - B_pos; - Real cos_alpha = AC.dot(BD) / AC.norm() / BD.norm(); - cos_alpha = Math::are_float_equal(cos_alpha, 1.) ? 1. : cos_alpha; - cos_alpha = Math::are_float_equal(cos_alpha, -1.) ? -1. : cos_alpha; - - /// check if the element is warped by evaluating the angle between the - /// normal to AB and the vector AC. - Real normal_dot_AC = normal_corrected.dot(AC); - Real volume_multiplicator = (normal_dot_AC <= 0) * 2 - 1; - - /// angle between two diagonals - Real alpha = acos(cos_alpha); - - Real current_volume = - 0.5 * AC.norm() * BD.norm() * sin(alpha) * volume_multiplicator; + Vector M_pos = (A_pos + B_pos) * 0.5; + Vector N_pos = (C_pos + D_pos) * 0.5; + Vector MN = M_pos - N_pos; + Vector AB = A_pos - B_pos; + Vector AB_0 = Vector(pos_it[A]) - Vector(pos_it[B]); + + // ASR volume computed as AB * thickness (AB * ratio) + Real ASR_volume = AB_0.norm() * AB_0.norm() * ASR_volume_ratio; + Real current_volume = AB.norm() * MN.norm(); current_volume = Math::are_float_equal(current_volume, 0.) ? 0 : current_volume; Real volume_change = current_volume - ASR_volume; - Real pressure_change = -volume_change / ASR_volume / this->compressibility; + Real pressure_change{0}; + if (volume_change < 0) { + pressure_change = -volume_change / ASR_volume / this->compressibility; + } dual = pressure_change * normal_corrected; } protected: SolidMechanicsModel & model; - const Real ASR_volume; + const Real ASR_volume_ratio; const std::string group_name; const Real compressibility; }; /* ------------------------------------------------------------------------ */ 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/mesh_utils/crack_numbers_updater.cc b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater.cc new file mode 100644 index 000000000..766816091 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater.cc @@ -0,0 +1,18 @@ +/* -------------------------------------------------------------------------- */ +#include "crack_numbers_updater.hh" +#include "element_synchronizer.hh" +#include "mesh_accessor.hh" +#include "mesh_utils.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +void CrackNumbersUpdater::communicateCrackNumbers() { + if (model.getMesh().getCommunicator().getNbProc() == 1) + return; + this->synchronizer.synchronizeOnce(*this, SynchronizationTag::_asr); +} + +} // namespace akantu diff --git a/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater.hh b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater.hh new file mode 100644 index 000000000..842254f54 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater.hh @@ -0,0 +1,54 @@ +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_CRACK_NUMBERS_UPDATER_HH__ +#define __AKANTU_CRACK_NUMBERS_UPDATER_HH__ + +/* -------------------------------------------------------------------------- */ +#include "data_accessor.hh" +#include "solid_mechanics_model.hh" + +/* -------------------------------------------------------------------------- */ +namespace akantu { +class ElementSynchronizer; +} // namespace akantu + +namespace akantu { + +class CrackNumbersUpdater : public DataAccessor { +public: + CrackNumbersUpdater(SolidMechanicsModel & model, + const ElementSynchronizer & synchronizer) + : model(model), synchronizer(synchronizer) {} + + void communicateCrackNumbers(); + + /* ------------------------------------------------------------------------ */ + /* Data Accessor inherited members */ + /* ------------------------------------------------------------------------ */ +public: + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; + + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; + + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; + + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +private: + /// Reference to the model + SolidMechanicsModel & model; + + /// distributed synchronizer to communicate nodes positions + const ElementSynchronizer & synchronizer; +}; + +} // namespace akantu + +#include "crack_numbers_updater_inline_impl.cc" + +#endif /* __AKANTU_CRACK_NUMBERS_UPDATER_HH__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.cc b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.cc new file mode 100644 index 000000000..495385ec7 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.cc @@ -0,0 +1,106 @@ +/* -------------------------------------------------------------------------- */ +#include "communicator.hh" +#include "crack_numbers_updater.hh" +#include "mesh.hh" +#include "mesh_accessor.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_CRACK_NUMBERS_UPDATER_INLINE_IMPL_CC__ +#define __AKANTU_CRACK_NUMBERS_UPDATER_INLINE_IMPL_CC__ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +inline UInt +CrackNumbersUpdater::getNbData(const Array & elements, + const SynchronizationTag & tag) const { + UInt size = 0; + if (tag == SynchronizationTag::_asr) { + size += sizeof(int); + // size += 2 * sizeof(Real); + auto & mesh = model.getMesh(); + for (auto elements_range : MeshElementsByTypes(elements)) { + auto type = elements_range.getType(); + if (mesh.getKind(type) == _ek_cohesive) { + UInt nb_elements = elements_range.getElements().size(); + UInt nb_quad_cohesive = + model.getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type); + size += nb_elements * nb_quad_cohesive * sizeof(Real); + } + } + } + return size; +} + +/* -------------------------------------------------------------------------- */ +inline void +CrackNumbersUpdater::packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const { + if (tag != SynchronizationTag::_asr) + return; + + int prank = model.getMesh().getCommunicator().whoAmI(); + buffer << prank; + + auto & mesh = model.getMesh(); + for (auto elements_range : MeshElementsByTypes(elements)) { + auto type = elements_range.getType(); + auto gt = elements_range.getGhostType(); + if ((mesh.getKind(type) == _ek_cohesive) and (gt == _not_ghost)) { + UInt nb_quad_cohesive = + model.getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type); + auto & coh_el_ids = elements_range.getElements(); + + for (UInt coh_el : coh_el_ids) { + const Array & material_index_vec = + model.getMaterialByElement(type, gt); + Material & material = model.getMaterial(material_index_vec(coh_el)); + auto & crack_nbs = material.getInternal("crack_number")(type, gt); + auto crack_nbs_it = crack_nbs.begin(); + for (auto i : arange(nb_quad_cohesive)) { + buffer << crack_nbs_it[coh_el * nb_quad_cohesive + i]; + } + } + } + } +} +/* -------------------------------------------------------------------------- + */ +inline void CrackNumbersUpdater::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { + if (tag != SynchronizationTag::_asr) + return; + + int proc; + buffer >> proc; + + auto & mesh = model.getMesh(); + for (auto elements_range : MeshElementsByTypes(elements)) { + auto type = elements_range.getType(); + auto gt = elements_range.getGhostType(); + if ((mesh.getKind(type) == _ek_cohesive) and (gt == _ghost)) { + UInt nb_quad_cohesive = + model.getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type); + auto & ghost_coh_el_ids = elements_range.getElements(); + for (UInt ghost_coh_el : ghost_coh_el_ids) { + const Array & material_index_vec = + model.getMaterialByElement(type, gt); + Material & material = + model.getMaterial(material_index_vec(ghost_coh_el)); + auto & crack_nbs = material.getInternal("crack_number")(type, gt); + auto crack_nbs_it = crack_nbs.begin(); + for (auto i : arange(nb_quad_cohesive)) { + Real unpacked_crack_nb; + buffer >> unpacked_crack_nb; + crack_nbs_it[ghost_coh_el * nb_quad_cohesive + i] = unpacked_crack_nb; + } + } + } + } +} + +} // namespace akantu + +#endif /* __AKANTU_CRACK_NUMBERS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.cc b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.cc new file mode 100644 index 000000000..f3df2ec54 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.cc @@ -0,0 +1,18 @@ +/* -------------------------------------------------------------------------- */ +#include "nodes_flag_updater.hh" +#include "element_synchronizer.hh" +#include "mesh_accessor.hh" +#include "mesh_utils.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +void NodesFlagUpdater::fillPreventInsertion() { + if (mesh.getCommunicator().getNbProc() == 1) + return; + this->synchronizer.slaveReductionOnce(*this, SynchronizationTag::_asr); +} + +} // namespace akantu diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.hh b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.hh new file mode 100644 index 000000000..710976066 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater.hh @@ -0,0 +1,66 @@ +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_NODES_FLAG_UPDATER_HH__ +#define __AKANTU_NODES_FLAG_UPDATER_HH__ + +/* -------------------------------------------------------------------------- */ +#include "data_accessor.hh" +#include "mesh.hh" + +/* -------------------------------------------------------------------------- */ +namespace akantu { +class ElementSynchronizer; +} // namespace akantu + +namespace akantu { + +class NodesFlagUpdater : public DataAccessor { +public: + NodesFlagUpdater(Mesh & mesh, ElementSynchronizer & synchronizer, + Array & prevent_insertion) + : mesh(mesh), synchronizer(synchronizer), + prevent_insertion(prevent_insertion) { + AKANTU_DEBUG_ASSERT( + mesh.getNbNodes() == prevent_insertion.size(), + "Array prevent_insertion does not have the same size as " + "the number of nodes in the mesh"); + } + + void fillPreventInsertion(); + + /* ------------------------------------------------------------------------ */ + /* Data Accessor inherited members */ + /* ------------------------------------------------------------------------ */ +public: + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; + + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; + + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; + + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +private: + /// Reference to the mesh to update + Mesh & mesh; + + /// distributed synchronizer to communicate nodes positions + ElementSynchronizer & synchronizer; + + /// Tells if a reduction is taking place or not + bool reduce{false}; + + /// tells at which nodes ASR segments shouldn't be inserted + Array & prevent_insertion; +}; + +} // namespace akantu + +#include "nodes_flag_updater_inline_impl.cc" + +#endif /* __AKANTU_NODES_FLAG_UPDATER_HH__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater_inline_impl.cc b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater_inline_impl.cc new file mode 100644 index 000000000..ba265eb10 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_flag_updater_inline_impl.cc @@ -0,0 +1,77 @@ +/* -------------------------------------------------------------------------- */ +#include "communicator.hh" +#include "mesh.hh" +#include "mesh_accessor.hh" +#include "nodes_flag_updater.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_NODES_FLAG_UPDATER_INLINE_IMPL_CC__ +#define __AKANTU_NODES_FLAG_UPDATER_INLINE_IMPL_CC__ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +inline UInt NodesFlagUpdater::getNbData(const Array & elements, + const SynchronizationTag & tag) const { + UInt size = 0; + if (tag == SynchronizationTag::_asr) { + size += + Mesh::getNbNodesPerElementList(elements) * sizeof(bool) + sizeof(int); + } + return size; +} + +/* -------------------------------------------------------------------------- */ +inline void NodesFlagUpdater::packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const { + if (tag != SynchronizationTag::_asr) + return; + + int prank = mesh.getCommunicator().whoAmI(); + buffer << prank; + + for (auto & element : elements) { + /// get element connectivity + const Vector current_conn = + const_cast(mesh).getConnectivity(element); + + /// loop on all connectivity nodes + for (auto node : current_conn) { + buffer << prevent_insertion(node); + } + } +} + +/* -------------------------------------------------------------------------- */ +inline void NodesFlagUpdater::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { + if (tag != SynchronizationTag::_asr) + return; + + MeshAccessor mesh_accessor(mesh); + auto dim = mesh.getSpatialDimension(); + auto & nodes_pos = mesh_accessor.getNodes(); + auto pos_it = nodes_pos.begin(dim); + + int proc; + buffer >> proc; + + for (auto & element : elements) { + /// get element connectivity + Vector current_conn = + const_cast(mesh).getConnectivity(element); + + /// loop on all connectivity nodes + for (auto node : current_conn) { + bool unpacked_node_flag; + buffer >> unpacked_node_flag; + prevent_insertion(node) = true; + } + } +} + +} // namespace akantu + +#endif /* __AKANTU_NODES_FLAG_UPDATER_INLINE_IMPL_CC__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.cc b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.cc new file mode 100644 index 000000000..b5047d7f0 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.cc @@ -0,0 +1,50 @@ +/** + * @file global_ids_updater.cc + * + * @author Marco Vocialta + * + * @date creation: Fri Apr 13 2012 + * @date last modification: Fri Dec 08 2017 + * + * @brief Functions of the GlobalIdsUpdater + * + * @section LICENSE + * + * Copyright (©) 2010-2018 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 "nodes_pos_updater.hh" +#include "element_synchronizer.hh" +#include "mesh_accessor.hh" +#include "mesh_utils.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +void NodesPosUpdater::updateNodesPos() { + this->reduce = true; + this->synchronizer.slaveReductionOnce(*this, SynchronizationTag::_asr); + + this->reduce = false; + this->synchronizer.synchronizeOnce(*this, SynchronizationTag::_asr); +} + +} // namespace akantu diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.hh b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.hh new file mode 100644 index 000000000..e3ed77957 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater.hh @@ -0,0 +1,97 @@ +/** + * @file global_ids_updater.hh + * + * @author Nicolas Richart + * @author Marco Vocialta + * + * @date creation: Fri Oct 02 2015 + * @date last modification: Fri Dec 08 2017 + * + * @brief Class that updates the global ids of new nodes that are + * inserted in the mesh. The functions in this class must be called + * after updating the node types + * + * @section LICENSE + * + * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_NODES_POS_UPDATER_HH__ +#define __AKANTU_NODES_POS_UPDATER_HH__ + +/* -------------------------------------------------------------------------- */ +#include "data_accessor.hh" +#include "mesh.hh" + +/* -------------------------------------------------------------------------- */ +namespace akantu { +class ElementSynchronizer; +} // namespace akantu + +namespace akantu { + +class NodesPosUpdater : public DataAccessor { +public: + NodesPosUpdater(Mesh & mesh, ElementSynchronizer & synchronizer, + const Array modified_pos) + : mesh(mesh), synchronizer(synchronizer), modified_pos(modified_pos) { + AKANTU_DEBUG_ASSERT(mesh.getNbNodes() == modified_pos.size(), + "Array modified_pos does not have the same size as " + "the number of nodes in the mesh"); + } + + void updateNodesPos(); + + /* ------------------------------------------------------------------------ */ + /* Data Accessor inherited members */ + /* ------------------------------------------------------------------------ */ +public: + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; + + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; + + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; + + /* ------------------------------------------------------------------------ */ + /* Members */ + /* ------------------------------------------------------------------------ */ +private: + /// Reference to the mesh to update + Mesh & mesh; + + /// distributed synchronizer to communicate nodes positions + ElementSynchronizer & synchronizer; + + /// Tells if a reduction is taking place or not + bool reduce{false}; + + /// which tells which nodes were modified + Array modified_pos; +}; + +} // namespace akantu + +#include "nodes_pos_updater_inline_impl.cc" + +#endif /* __AKANTU_NODES_POS_UPDATER_HH__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater_inline_impl.cc b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater_inline_impl.cc new file mode 100644 index 000000000..9423421a0 --- /dev/null +++ b/extra_packages/asr-tools/src/mesh_utils/nodes_pos_updater_inline_impl.cc @@ -0,0 +1,131 @@ +/** + * @file global_ids_updater_inline_impl.cc + * + * @author Marco Vocialta + * + * @date creation: Fri Oct 02 2015 + * @date last modification: Sun Aug 13 2017 + * + * @brief Implementation of the inline functions of GlobalIdsUpdater + * + * @section LICENSE + * + * Copyright (©) 2015-2018 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 "communicator.hh" +#include "mesh.hh" +#include "mesh_accessor.hh" +#include "nodes_pos_updater.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_NODES_POS_UPDATER_INLINE_IMPL_CC__ +#define __AKANTU_NODES_POS_UPDATER_INLINE_IMPL_CC__ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +inline UInt NodesPosUpdater::getNbData(const Array & elements, + const SynchronizationTag & tag) const { + UInt size = 0; + if (tag == SynchronizationTag::_asr) { + auto dim = mesh.getSpatialDimension(); + size += Mesh::getNbNodesPerElementList(elements) * dim * sizeof(Real) + + Mesh::getNbNodesPerElementList(elements) * sizeof(bool) + + sizeof(int); + } + return size; +} + +/* -------------------------------------------------------------------------- */ +inline void NodesPosUpdater::packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const { + if (tag != SynchronizationTag::_asr) + return; + + int prank = mesh.getCommunicator().whoAmI(); + buffer << prank; + + auto dim = mesh.getSpatialDimension(); + auto & positions = mesh.getNodes(); + auto pos_it = positions.begin(dim); + + for (auto & element : elements) { + /// get element connectivity + const Vector current_conn = + const_cast(mesh).getConnectivity(element); + + /// loop on all connectivity nodes + for (auto node : current_conn) { + Vector coord(pos_it[node]); + buffer << coord; + buffer << modified_pos(node); + } + } +} + +/* -------------------------------------------------------------------------- */ +inline void NodesPosUpdater::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { + if (tag != SynchronizationTag::_asr) + return; + + MeshAccessor mesh_accessor(mesh); + auto dim = mesh.getSpatialDimension(); + auto & nodes_pos = mesh_accessor.getNodes(); + auto pos_it = nodes_pos.begin(dim); + + int proc; + buffer >> proc; + + for (auto & element : elements) { + /// get element connectivity + Vector current_conn = + const_cast(mesh).getConnectivity(element); + + /// loop on all connectivity nodes + for (auto node : current_conn) { + Vector current_pos(pos_it[node]); + auto current_modif_flag = modified_pos(node); + Vector unpacked_pos(dim); + bool unpacked_mod_flag; + buffer >> unpacked_pos; + buffer >> unpacked_mod_flag; + if (this->reduce) { + /// if slave is modified -> take position of slave; if both slave and + /// master are modified take the average; if master is modified or + /// noneof them are modified -> do nothing + if (unpacked_mod_flag and not current_modif_flag) { + current_pos = unpacked_pos; + } else if (unpacked_mod_flag and current_modif_flag) { + current_pos += unpacked_pos; + current_pos /= 2; + } + } else { + current_pos = unpacked_pos; + } + } + } +} + +} // namespace akantu + +#endif /* __AKANTU_NODES_POS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.cc b/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.cc index f9162f52d..321fdb1d8 100644 --- a/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.cc +++ b/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.cc @@ -1,297 +1,410 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Wed Feb 21 2018 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2018 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 "material_cohesive_linear_sequential.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearSequential:: MaterialCohesiveLinearSequential(SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model, id), normal_stresses("normal_stresses", *this), normal_tractions("normal_tractions", *this), - effective_stresses("effective_stresses", *this) { + effective_stresses("effective_stresses", *this), + crack_number("crack_number", *this) { AKANTU_DEBUG_IN(); this->registerParam( "insertion_threshold", insertion_threshold, Real(0.0), _pat_parsable | _pat_readable, "Portion of elements below the highest stress to be inserted"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearSequential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::initMaterial(); normal_stresses.initialize(1); normal_tractions.initialize(spatial_dimension); effective_stresses.initialize(1); + crack_number.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearSequential::checkInsertion( bool check_only) { AKANTU_DEBUG_IN(); - const Mesh & mesh_facets = this->model->getMeshFacets(); - CohesiveElementInserter & inserter = this->model->getElementInserter(); + Mesh & mesh_facets = this->model->getMeshFacets(); + MeshUtils::fillElementToSubElementsData(mesh_facets); for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { - ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); - const auto & facets_check = inserter.getCheckFacets(type_facet); - auto & f_insertion = inserter.getInsertionFacets(type_facet); - auto & f_filter = this->facet_filter(type_facet); - auto & sig_c_eff = this->sigma_c_eff(type_cohesive); - auto & del_c = this->delta_c_eff(type_cohesive); - auto & ins_stress = this->insertion_stress(type_cohesive); - auto & trac_old = this->tractions.previous(type_cohesive); - auto & norm_stresses = normal_stresses(type_facet); - auto & norm_tractions = normal_tractions(type_facet); - const auto & f_stress = this->model->getStressOnFacets(type_facet); - const auto & sigma_limits = this->sigma_c(type_facet); - auto & eff_stresses = effective_stresses(type_facet); - auto & facet_conn = mesh_facets.getConnectivity(type_facet); - - UInt nb_quad_facet = this->model->getFEEngine("FacetsFEEngine") - .getNbIntegrationPoints(type_facet); + + // find the facet with the highest tensile stress + auto max_stress_data = findCriticalFacet(type_facet); + auto max_stress_facet = std::get<0>(max_stress_data); + auto max_stress = std::get<1>(max_stress_data); + auto crack_number = std::get<2>(max_stress_data); + + // communicate between processors the highest stress + Real local_max_stress = max_stress; + auto && comm = akantu::Communicator::getWorldCommunicator(); + comm.allReduce(max_stress, SynchronizerOperation::_max); + + // if the effective max stress is 0 or < 1 ->skip + if (max_stress < 1) + continue; + + // if the max stress at this proc != global max stress -> skip + if (local_max_stress != max_stress) + continue; + + // duplicate single facet and insert cohesive element + insertSingleCohesiveElement(type_facet, max_stress_facet, crack_number, + check_only); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------*/ +template +std::tuple +MaterialCohesiveLinearSequential::findCriticalFacet( + const ElementType & type_facet) { + AKANTU_DEBUG_IN(); + + Mesh & mesh = this->model->getMesh(); + Mesh & mesh_facets = this->model->getMeshFacets(); + MeshUtils::fillElementToSubElementsData(mesh_facets); + const auto & pos = mesh.getNodes(); + const auto pos_it = make_view(pos, spatial_dimension).begin(); + CohesiveElementInserter & inserter = this->model->getElementInserter(); + CSR nodes_to_elements; + MeshUtils::buildNode2Elements(mesh, nodes_to_elements, spatial_dimension, + _ek_cohesive); + + UInt max_stress_facet(-1); + Real max_stress = std::numeric_limits::min(); + Real potential_crack_nb(-1); + + ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); + const auto & facets_check = inserter.getCheckFacets(type_facet); + auto & f_filter = this->facet_filter(type_facet); + auto & norm_stresses = normal_stresses(type_facet); + auto & norm_tractions = normal_tractions(type_facet); + const auto & f_stress = this->model->getStressOnFacets(type_facet); + const auto & sigma_limits = this->sigma_c(type_facet); + auto & eff_stresses = effective_stresses(type_facet); + auto & facet_conn = mesh_facets.getConnectivity(type_facet); + + UInt nb_quad_facet = this->model->getFEEngine("FacetsFEEngine") + .getNbIntegrationPoints(type_facet); #ifndef AKANTU_NDEBUG - UInt nb_quad_cohesive = this->model->getFEEngine("CohesiveFEEngine") - .getNbIntegrationPoints(type_cohesive); + UInt nb_quad_cohesive = this->model->getFEEngine("CohesiveFEEngine") + .getNbIntegrationPoints(type_cohesive); - AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, - "The cohesive element and the corresponding facet do " - "not have the same numbers of integration points"); + AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, + "The cohesive element and the corresponding facet do " + "not have the same numbers of integration points"); #endif + // skip if no facets of this type are present + UInt nb_facet = f_filter.size(); + if (nb_facet == 0) + return std::make_tuple(max_stress_facet, max_stress, potential_crack_nb); + + Matrix stress_tmp(spatial_dimension, spatial_dimension); + UInt sp2 = spatial_dimension * spatial_dimension; + + const auto & tangents = this->model->getTangents(type_facet); + const auto & normals = this->model->getFEEngine("FacetsFEEngine") + .getNormalsOnIntegrationPoints(type_facet); + auto normal_begin = normals.begin(spatial_dimension); + auto tangent_begin = tangents.begin(tangents.getNbComponent()); + auto facet_stress_begin = + f_stress.begin(spatial_dimension, spatial_dimension * 2); + + // loop over each facet belonging to this material + for (auto && data : + zip(f_filter, sigma_limits, eff_stresses, + make_view(norm_stresses, nb_quad_facet), + make_view(norm_tractions, spatial_dimension, nb_quad_facet))) { + auto facet = std::get<0>(data); + auto & sigma_limit = std::get<1>(data); + auto & eff_stress = std::get<2>(data); + auto & stress_check = std::get<3>(data); + auto & normal_traction = std::get<4>(data); + + // skip facets where check shouldn't be inserted (or already inserted) + if (!facets_check(facet)) + continue; - Matrix stress_tmp(spatial_dimension, spatial_dimension); - UInt sp2 = spatial_dimension * spatial_dimension; - - const auto & tangents = this->model->getTangents(type_facet); - const auto & normals = this->model->getFEEngine("FacetsFEEngine") - .getNormalsOnIntegrationPoints(type_facet); - auto normal_begin = normals.begin(spatial_dimension); - auto tangent_begin = tangents.begin(tangents.getNbComponent()); - auto facet_stress_begin = - f_stress.begin(spatial_dimension, spatial_dimension * 2); - - std::vector new_sigmas; - std::vector> new_normal_traction; - std::vector new_delta_c; - std::vector> facet_stress; - - // loop over each facet belonging to this material - for (auto && data : - zip(f_filter, sigma_limits, eff_stresses, - make_view(norm_stresses, nb_quad_facet), - make_view(norm_tractions, spatial_dimension, nb_quad_facet))) { - auto facet = std::get<0>(data); - auto & sigma_limit = std::get<1>(data); - auto & eff_stress = std::get<2>(data); - auto & stress_check = std::get<3>(data); - auto & normal_traction = std::get<4>(data); - - // skip facets where check shouldn't be realized - if (!facets_check(facet)) - continue; - - // compute the effective norm on each quadrature point of the facet - for (UInt q : arange(nb_quad_facet)) { - UInt current_quad = facet * nb_quad_facet + q; - const Vector & normal = normal_begin[current_quad]; - const Vector & tangent = tangent_begin[current_quad]; - const Matrix & facet_stress_it = facet_stress_begin[current_quad]; - - // compute average stress on the current quadrature point - Matrix stress_1(facet_stress_it.storage(), spatial_dimension, - spatial_dimension); - - Matrix stress_2(facet_stress_it.storage() + sp2, - spatial_dimension, spatial_dimension); - - stress_tmp.copy(stress_1); - stress_tmp += stress_2; - stress_tmp /= 2.; - - Vector normal_traction_vec(normal_traction(q)); - - // compute normal and effective stress - stress_check(q) = this->computeEffectiveNorm( - stress_tmp, normal, tangent, normal_traction_vec); + // check to how many cohesive elements this facet is connected + std::vector> coh_neighbors(2); + Real coh_crack_nb(-1); + Array facet_nodes(2); + for (UInt i : arange(2)) { + facet_nodes(i) = facet_conn(facet, i); + for (auto & elem : nodes_to_elements.getRow(facet_nodes(i))) { + coh_neighbors[i].emplace(elem); } - - // verify if the effective stress overcomes the threshold - Real final_stress = stress_check.mean(); - if (this->max_quad_stress_insertion) - final_stress = *std::max_element( - stress_check.storage(), stress_check.storage() + nb_quad_facet); - - // normalize by the limit stress - eff_stress = final_stress / sigma_limit; - facet_stress.push_back(std::make_pair(facet, eff_stress)); } - // sort facets by the value of effective stress - std::sort(facet_stress.begin(), facet_stress.end(), - [](const std::pair & lhs, - const std::pair & rhs) { - return lhs.second > rhs.second; - }); + // see if a single tip crack is present + bool single_tip{false}; + Array single_tip_node; + for (UInt i : arange(2)) { + if (coh_neighbors[i].size() == 1) { + single_tip = true; + single_tip_node.push_back(i); + } + } - // continue to the next el type if no elements to insert - if (facet_stress[0].second <= 1.) + // if no coh els are connected or no single tip present - skip facet + if (not single_tip) continue; - // second loop to activate certain portion of elements - const Mesh & mesh = this->model->getMesh(); - CSR nodes_to_elements; - MeshUtils::buildNode2ElementsElementTypeMap(mesh, nodes_to_elements, - _cohesive_2d_4, _not_ghost); - - // for (auto && data : - // zip(f_filter, sigma_limits, eff_stresses, - // make_view(norm_stresses, nb_quad_facet), - // make_view(norm_tractions, spatial_dimension, nb_quad_facet))) { - UInt nb_coh_inserted{0}; - for (auto && pair : facet_stress) { - // allow insertion of only 1 cohesive element - if (nb_coh_inserted == 1) - break; - - auto & facet = std::get<0>(pair); - // get facet's local id - auto local_id = f_filter.find(facet); - AKANTU_DEBUG_ASSERT(local_id != UInt(-1), - "mismatch between global and local facet numbering"); - - // skip facets where check shouldn't be realized - if (!facets_check(facet)) - continue; - - auto & eff_stress = std::get<1>(pair); - // break when descending eff_stresses drop below 1 - if (eff_stress < 1.) - break; - - auto & sigma_limit = sigma_limits(local_id); - Vector stress_check = - make_view(norm_stresses, nb_quad_facet).begin()[local_id]; - Matrix normal_traction = - make_view(norm_tractions, spatial_dimension, nb_quad_facet) - .begin()[local_id]; - - // Real reduced_threshold = - // max_eff_stress - (max_eff_stress - 1) * this->insertion_threshold; - // reduced_threshold = std::max(reduced_threshold, 1.); - - // if (eff_stress >= reduced_threshold) { - // check if this facet is connected to a single cohesive else - std::set coh_neighbors; - for (UInt i : arange(2)) { - const UInt facet_node = facet_conn(facet, i); - for (auto & elem : nodes_to_elements.getRow(facet_node)) { - // if (mesh.getKind(elem.type) == _ek_cohesive) - coh_neighbors.emplace(elem); + // if single tip is present, compute angle with the candidate facet; if + // two single tips - compute two angles; if angle(s) are sharp - skip + bool sharp_angle{false}; + for (auto tip_node : single_tip_node) { + auto coh_el = *coh_neighbors[tip_node].begin(); + auto & crack_nbs = this->crack_number(coh_el.type, coh_el.ghost_type); + auto crack_nb_it = crack_nbs.begin(); + auto crack_nb = crack_nb_it[coh_el.element * nb_quad_cohesive]; + coh_crack_nb = crack_nb; + + auto & sub_to_element = + mesh_facets.getSubelementToElement(coh_el.type, coh_el.ghost_type); + auto sub_to_el_it = sub_to_element.begin(sub_to_element.getNbComponent()); + const Vector & subelements_to_element = + sub_to_el_it[coh_el.element]; + auto & coh_facet_conn = + mesh_facets.getConnectivity(type_facet, coh_el.ghost_type); + UInt common_node(-1); + UInt candidate_facet_node(-1); + UInt coh_facet_node(-1); + bool nodes_found{false}; + + // iterate through 2 facets of cohesive element + for (auto & coh_facet : subelements_to_element) { + Array coh_facet_nodes(2); + for (UInt node : arange(2)) { + coh_facet_nodes(node) = coh_facet_conn(coh_facet.element, node); + } + // find which node is common and which are unique + for (UInt i : arange(2)) { + auto id = coh_facet_nodes.find(facet_nodes(i)); + if (id == UInt(-1)) + continue; + nodes_found = true; + common_node = facet_nodes(i); + coh_facet_node = coh_facet_nodes(-1 * id + 1); + candidate_facet_node = facet_nodes(-1 * i + 1); } + if (nodes_found) + break; } - // if no or more than one coh neighbors - continue - if (coh_neighbors.size() != 1) - continue; - - f_insertion(facet) = true; - nb_coh_inserted++; + AKANTU_DEBUG_ASSERT(nodes_found, + "Common and unique nodes were not found between " + "candidate facet and a cohesive one"); + + // determine the direction of each facet (works only in 2D) + Vector coh_facet_dir(spatial_dimension); + Vector candidate_facet_dir(spatial_dimension); + coh_facet_dir = pos_it[coh_facet_node]; + coh_facet_dir -= Vector(pos_it[common_node]); + candidate_facet_dir = pos_it[candidate_facet_node]; + candidate_facet_dir -= Vector(pos_it[common_node]); + + // compute dot product between 2 vectors and discard sharp angles + Real dot = coh_facet_dir.dot(candidate_facet_dir); + if (dot >= 0.) + sharp_angle = true; + } - if (check_only) - continue; + if (sharp_angle) + continue; - // store the new cohesive material parameters for each quadrature - // point - for (UInt q = 0; q < nb_quad_facet; ++q) { - Real new_sigma = stress_check(q); - Vector normal_traction_vec(normal_traction(q)); + // compute the effective norm on each quadrature point of the facet + for (UInt q : arange(nb_quad_facet)) { + UInt current_quad = facet * nb_quad_facet + q; + const Vector & normal = normal_begin[current_quad]; + const Vector & tangent = tangent_begin[current_quad]; + const Matrix & facet_stress_it = facet_stress_begin[current_quad]; - if (spatial_dimension != 3) - normal_traction_vec *= -1.; + // compute average stress on the current quadrature point + Matrix stress_1(facet_stress_it.storage(), spatial_dimension, + spatial_dimension); - new_sigmas.push_back(new_sigma); - new_normal_traction.push_back(normal_traction_vec); + Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, + spatial_dimension); - Real new_delta; + stress_tmp.copy(stress_1); + stress_tmp += stress_2; + stress_tmp /= 2.; - // set delta_c in function of G_c or a given delta_c value - if (Math::are_float_equal(this->delta_c, 0.)) - new_delta = 2 * this->G_c / new_sigma; - else - new_delta = sigma_limit / new_sigma * this->delta_c; + Vector normal_traction_vec(normal_traction(q)); - new_delta_c.push_back(new_delta); - } - // } + // compute normal and effective stress + stress_check(q) = this->computeEffectiveNorm(stress_tmp, normal, tangent, + normal_traction_vec); } - // update material data for the new elements - UInt old_nb_quad_points = sig_c_eff.size(); - UInt new_nb_quad_points = new_sigmas.size(); - sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); - ins_stress.resize(old_nb_quad_points + new_nb_quad_points); - trac_old.resize(old_nb_quad_points + new_nb_quad_points); - del_c.resize(old_nb_quad_points + new_nb_quad_points); - - for (UInt q = 0; q < new_nb_quad_points; ++q) { - sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; - del_c(old_nb_quad_points + q) = new_delta_c[q]; - for (UInt dim = 0; dim < spatial_dimension; ++dim) { - ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); - trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); - } + // verify if the effective stress overcomes the threshold + Real final_stress = stress_check.mean(); + if (this->max_quad_stress_insertion) + final_stress = *std::max_element(stress_check.storage(), + stress_check.storage() + nb_quad_facet); + + // normalize by the limit stress + eff_stress = final_stress / sigma_limit; + + // update the max stress and corresponding facet + if (eff_stress > max_stress) { + max_stress_facet = facet; + max_stress = eff_stress; + potential_crack_nb = coh_crack_nb; } } + // if (local_max_stress == max_stress) + return std::make_tuple(max_stress_facet, max_stress, potential_crack_nb); + // else + // return std::make_tuple(-1, -1, -1); +} + +/* -------------------------------------------------------------------------*/ +template +void MaterialCohesiveLinearSequential:: + insertSingleCohesiveElement(const ElementType & type_facet, UInt facet_nb, + Real crack_nb, bool check_only) { + + CohesiveElementInserter & inserter = this->model->getElementInserter(); + + ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); + auto & f_filter = this->facet_filter(type_facet); + auto & f_insertion = inserter.getInsertionFacets(type_facet); + auto & norm_stresses = normal_stresses(type_facet); + auto & norm_tractions = normal_tractions(type_facet); + const auto & sigma_limits = this->sigma_c(type_facet); + auto & sig_c_eff = this->sigma_c_eff(type_cohesive); + auto & del_c = this->delta_c_eff(type_cohesive); + auto & ins_stress = this->insertion_stress(type_cohesive); + auto & trac_old = this->tractions.previous(type_cohesive); + auto & crack_nbs_current_facet = this->crack_number(type_cohesive); + UInt nb_quad_facet = this->model->getFEEngine("FacetsFEEngine") + .getNbIntegrationPoints(type_facet); + + Vector new_sigmas(nb_quad_facet); + Array new_normal_traction(nb_quad_facet, spatial_dimension); + Vector new_delta_c(nb_quad_facet); + Vector new_crack_nb(nb_quad_facet); + + // get facet's local id + auto local_id = f_filter.find(facet_nb); + AKANTU_DEBUG_ASSERT(local_id != UInt(-1), + "mismatch between global and local facet numbering"); + + // mark the insertion of the cohesive element + f_insertion(facet_nb) = true; + + if (check_only) + return; + + auto & sigma_limit = sigma_limits(local_id); + Vector stress_check = + make_view(norm_stresses, nb_quad_facet).begin()[local_id]; + Matrix normal_traction = + make_view(norm_tractions, spatial_dimension, nb_quad_facet) + .begin()[local_id]; + + // store the new cohesive material parameters for each quad point + for (UInt q = 0; q < nb_quad_facet; ++q) { + Real new_sigma = stress_check(q); + Vector normal_traction_vec(normal_traction(q)); + + if (spatial_dimension != 3) + normal_traction_vec *= -1.; + + new_sigmas(q) = new_sigma; + new_crack_nb(q) = crack_nb; + for (UInt i : arange(spatial_dimension)) + new_normal_traction(q, i) = normal_traction_vec(i); + + Real new_delta; + + // set delta_c in function of G_c or a given delta_c value + if (Math::are_float_equal(this->delta_c, 0.)) + new_delta = 2 * this->G_c / new_sigma; + else + new_delta = sigma_limit / new_sigma * this->delta_c; + + new_delta_c(q) = new_delta; + } + + // update material data for the new elements + UInt old_nb_quad_points = sig_c_eff.size(); + sig_c_eff.resize(old_nb_quad_points + nb_quad_facet); + ins_stress.resize(old_nb_quad_points + nb_quad_facet); + trac_old.resize(old_nb_quad_points + nb_quad_facet); + del_c.resize(old_nb_quad_points + nb_quad_facet); + crack_nbs_current_facet.resize(old_nb_quad_points + nb_quad_facet); + + for (UInt q = 0; q < nb_quad_facet; ++q) { + sig_c_eff(old_nb_quad_points + q) = new_sigmas(q); + del_c(old_nb_quad_points + q) = new_delta_c(q); + crack_nbs_current_facet(old_nb_quad_points + q) = new_crack_nb(q); + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + ins_stress(old_nb_quad_points + q, dim) = new_normal_traction(q, dim); + trac_old(old_nb_quad_points + q, dim) = new_normal_traction(q, dim); + } + } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear_sequential, MaterialCohesiveLinearSequential); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.hh b/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.hh index 8a1714251..ef557f3cd 100644 --- a/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.hh +++ b/extra_packages/extra-materials/src/material_cohesive/material_cohesive_linear_sequential.hh @@ -1,102 +1,115 @@ /** * @file material_cohesive_linear.hh * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Feb 21 2018 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2018 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 "material_cohesive_linear.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_SEQUENTIAL_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_SEQUENTIAL_HH__ namespace akantu { /** * Cohesive material linear damage for extrinsic case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: * 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) * - penalty : stiffness in compression to prevent penetration */ template class MaterialCohesiveLinearSequential : public MaterialCohesiveLinear { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ + /* -----------------------------------------------------------------------*/ + /* Constructors/Destructors */ + /* -----------------------------------------------------------------------*/ public: MaterialCohesiveLinearSequential(SolidMechanicsModel & model, const ID & id = ""); - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ + /* -----------------------------------------------------------------------*/ + /* Methods */ + /* -----------------------------------------------------------------------*/ public: /// initialise material void initMaterial() override; /// check stress for cohesive elements' insertion void checkInsertion(bool check_only = false) override; + /// identify a facet with the highest tensile stress complying with + /// topological criterion + std::tuple + findCriticalFacet(const ElementType & type_facet); + + /// duplicate single facet and insert cohesive if check_only!=false + void insertSingleCohesiveElement(const ElementType & type_facet, + UInt facet_nb, Real crack_nb, + bool check_only); + protected: - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ + /* -----------------------------------------------------------------------*/ + /* Accessors */ + /* -----------------------------------------------------------------------*/ - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ + /* -----------------------------------------------------------------------*/ + /* Class Members */ + /* -----------------------------------------------------------------------*/ protected: /// portion of elements with the stress above their threshold to be inserted Real insertion_threshold; /// normal stresses on facets FacetInternalField normal_stresses; /// normal stresses on facets FacetInternalField normal_tractions; /// normal stresses normalized by the stress limit FacetInternalField effective_stresses; + + /// crack number + CohesiveInternalField crack_number; }; -/* -------------------------------------------------------------------------- */ -/* inline functions */ -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------*/ +/* inline functions */ +/* -------------------------------------------------------------------------*/ } // namespace akantu #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_SEQUENTIAL_HH__ */ diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 895a6dec1..6497f05d6 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,633 +1,637 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2018 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 . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { __attribute__((unused)) constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG __attribute__((unused)) constexpr Real REAL_INIT_VALUE{0.}; #else __attribute__((unused)) constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; using MemoryID = UInt; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (fluid_diffusion_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (ce_insertion_order) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ (htm_pressure) \ (htm_gradient_pressure) \ (htm_phi) \ (htm_gradient_phi) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ - (solver_solution) + (solver_solution) \ + (asr) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups _ce_insertion_order, ///< synchronization of the order of insertion of - ///cohesive elements + /// cohesive elements // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- FluidDiffusion tags --- _htm_pressure, ///< synchronization of the nodal pressure _htm_gradient_pressure, ///< synchronization of the element gradient /// pressure // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- - _solver_solution ///< synchronization of the solution obained with the - /// PETSc solver + _solver_solution, ///< synchronization of the solution obained with the + /// PETSc solver + + // --- ASRTools tags --- + _asr ///< synchronization for asr (nodes, node flags) }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT ' ' #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic * number is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */