diff --git a/extra_packages/asr-tools/src/asr_tools.cc b/extra_packages/asr-tools/src/asr_tools.cc index 197ec6b48..30adaa91b 100644 --- a/extra_packages/asr-tools/src/asr_tools.cc +++ b/extra_packages/asr-tools/src/asr_tools.cc @@ -1,5745 +1,5943 @@ /** * @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 #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_cohesive_linear_sequential.hh" #include "material_damage_iterative_orthotropic.hh" #include "material_iterative_stiffness_reduction.hh" #include "mesh_utils.hh" #include "node_synchronizer.hh" #include "nodes_eff_stress_updater.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 #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ASRTools::ASRTools(SolidMechanicsModel & model) : model(model), volume(0.), nb_dumps(0), cohesive_insertion(false), - doubled_facets_ready(false), doubled_nodes_ready(false), + doubled_facets_ready(false), asr_central_nodes_ready(false), disp_stored(0, model.getSpatialDimension()), ext_force_stored(0, model.getSpatialDimension()), boun_stored(0, model.getSpatialDimension()), tensile_homogenization(false), modified_pos(false), partition_border_nodes(false), nodes_eff_stress(0), ASR_nodes(false) { // register event handler for asr tools auto & mesh = model.getMesh(); - // auto const dim = model.getSpatialDimension(); - mesh.registerEventHandler(*this, _ehp_lowest); + mesh.registerEventHandler(*this, _ehp_synchronizer); + // mesh.registerEventHandler(*this, _ehp_lowest); if (mesh.hasMeshFacets()) { - mesh.getMeshFacets().registerEventHandler(*this, _ehp_lowest); + mesh.getMeshFacets().registerEventHandler(*this, _ehp_synchronizer); + // 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); partition_border_nodes.resize(nb_nodes); partition_border_nodes.set(false); nodes_eff_stress.resize(nb_nodes); nodes_eff_stress.set(0); 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_regular)) { 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 (!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.empty()) { 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.empty()) { 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 (const 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 (const 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; } } 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 - left bottom corner if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // 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 - 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 - 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.zero(); H(1, 1) = 0.001 * (2 * tensile_test - 1); performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.zero(); 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() { model.dump(); } /* ------------------------------------------------------------------------ */ // 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::insertASRCohesivesRandomly(const UInt & nb_insertions, std::string facet_mat_name, Real gap_ratio) { AKANTU_DEBUG_IN(); // fill in the partition_border_nodes array communicateFlagsOnNodes(); /// pick central facets and neighbors pickFacetsRandomly(nb_insertions, facet_mat_name); insertOppositeFacetsAndCohesives(); assignCrackNumbers(); insertGap(gap_ratio); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------- */ void ASRTools::insertASRCohesivesRandomly3D(const UInt & nb_insertions, std::string facet_mat_name, Real /*gap_ratio*/) { AKANTU_DEBUG_IN(); // fill in the partition_border_nodes array communicateFlagsOnNodes(); /// pick central facets and neighbors pickFacetsRandomly(nb_insertions, facet_mat_name); insertOppositeFacetsAndCohesives(); assignCrackNumbers(); // insertGap3D(gap_ratio); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------- */ void ASRTools::insertASRCohesiveLoops3D(const UInt & nb_insertions, std::string facet_mat_name, Real /*gap_ratio*/) { AKANTU_DEBUG_IN(); // fill in the partition_border_nodes array communicateFlagsOnNodes(); /// pick central facets and neighbors auto inserted = closedFacetsLoopAroundPoint(nb_insertions, facet_mat_name); auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(inserted, SynchronizerOperation::_sum); AKANTU_DEBUG_ASSERT(inserted, "No ASR sites were inserted across the domain"); insertOppositeFacetsAndCohesives(); assignCrackNumbers(); communicateCrackNumbers(); // insertGap3D(gap_ratio); // preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------- */ template void ASRTools::insertPreCracks(const UInt & nb_insertions, std::string facet_mat_name) { AKANTU_DEBUG_IN(); // fill in the partition_border_nodes array communicateFlagsOnNodes(); /// pick central facets and neighbors auto inserted = insertCohesiveLoops(nb_insertions, facet_mat_name); auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(inserted, SynchronizerOperation::_sum); AKANTU_DEBUG_ASSERT(inserted, "No ASR sites were inserted across the domain"); communicateCrackNumbers(); AKANTU_DEBUG_OUT(); } template void ASRTools::insertPreCracks<2>(const UInt & nb_insertions, std::string mat_name); template void ASRTools::insertPreCracks<3>(const UInt & nb_insertions, std::string mat_name); /* ------------------------------------------------------------------------- */ void ASRTools::insertASRCohesivesByCoords(const Matrix & positions, Real gap_ratio) { AKANTU_DEBUG_IN(); // fill in the border_nodes array communicateFlagsOnNodes(); /// pick the facets to duplicate pickFacetsByCoord(positions); insertOppositeFacetsAndCohesives(); 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->partition_border_nodes); nodes_flag_updater.fillPreventInsertion(); } /* ------------------------------------------------------------------- */ void ASRTools::communicateEffStressesOnNodes() { auto & mesh = model.getMesh(); auto & synch = mesh.getElementSynchronizer(); NodesEffStressUpdater nodes_eff_stress_updater(mesh, synch, this->nodes_eff_stress); nodes_eff_stress_updater.updateMaxEffStressAtNodes(); } /* ------------------------------------------------------------------- */ 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) { 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 (const 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) { 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); /// 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 (this->partition_border_nodes(node)) border_facets = true; } if (border_facets) { std::cout << "Facet at the position " << position(0) << "," << position(1) << " is close to the partition border. Skipping this position" << std::endl; continue; } /// 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) { 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) { 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) { std::cout << "Proc " << prank << " couldn't place " << nb_insertions << " ASR sites" << std::endl; return; } 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 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 (this->partition_border_nodes(node) or this->ASR_nodes(node)) border_facets = true; } 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; } } } /* ------------------------------------------------------------------- */ UInt ASRTools::closedFacetsLoopAroundPoint(UInt nb_insertions, std::string 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 element_type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(element_type); auto & doubled_facets = mesh_facets.createElementGroup("doubled_facets", dim - 1); auto & doubled_nodes = mesh.createNodeGroup("doubled_nodes"); auto material_id = model.getMaterialIndex(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 & node_conn = mesh_facets.getConnectivity(_point_1, gt); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); const Real max_dot{0.7}; if (not nb_insertions) return 0; CSR nodes_to_facets; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_facets, dim - 1); std::set matrix_nodes; for_each_element( mesh_facets, [&](auto && point_el) { // get the node number auto node = node_conn(point_el.element); // discard ghost nodes if (mesh.isPureGhostNode(node)) goto nextnode; for (auto & facet : nodes_to_facets.getRow(node)) { auto & facet_material = coh_model.getFacetMaterial( facet.type, facet.ghost_type)(facet.element); if (facet_material != material_id) { goto nextnode; } } matrix_nodes.emplace(node); nextnode:; }, _spatial_dimension = 0, _ghost_type = _not_ghost); UInt nb_nodes = matrix_nodes.size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_nodes - 1); if (not nb_nodes) { std::cout << "Proc " << prank << " couldn't place " << nb_insertions << " ASR sites" << std::endl; return 0; } UInt already_inserted = 0; std::set left_nodes = matrix_nodes; while (already_inserted < nb_insertions) { if (not left_nodes.size()) { std::cout << "Proc " << prank << " inserted " << already_inserted << " ASR sites out of " << nb_insertions << std::endl; return already_inserted; } auto id = dis(random_generator); auto it = matrix_nodes.begin(); std::advance(it, id); UInt cent_node(*it); left_nodes.erase(*it); // not touching other cracks if (this->ASR_nodes(cent_node)) continue; // поехали CSR segments_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, segments_to_nodes, dim - 2); auto && segments_to_node = segments_to_nodes.getRow(cent_node); Array segments_list; for (auto && segment : segments_to_node) { segments_list.push_back(segment); } // pick the first non-ghost segment in the list Element starting_segment; for (auto segment : segments_list) { if (segment.ghost_type == _not_ghost) { starting_segment = segment; break; } } // go to next node if starting segment could not be picked if (starting_segment == ElementNull) { std::cout << "Could not pick the starting segment, switching to the next node" << std::endl; continue; } auto & facets_to_segment = mesh_facets.getElementToSubelement(starting_segment); // pick the first good facet in the list Element starting_facet{ElementNull}; for (const auto & facet : facets_to_segment) { // check if the facet and its nodes are ok if (isFacetAndNodesGood(facet, material_id, true)) { starting_facet = facet; break; } } // go to next node if starting facet could not be picked if (starting_facet == ElementNull) { std::cout << "Could not pick the starting facet, switching to the next node" << std::endl; continue; } // call the actual function that build the bridge auto facets_in_loop = findFacetsLoopFromSegment2Segment( starting_facet, starting_segment, starting_segment, cent_node, max_dot, material_id, true); // loop broke -> try another point if (not facets_in_loop.size()) { std::cout << "Couldn't close the facet loop, taking the next node" << std::endl; continue; } // otherwise champaigne for (auto & facet : facets_in_loop) { doubled_facets.add(facet); /// add all facet nodes to the group for (auto node : arange(nb_nodes_facet)) { doubled_nodes.add(facet_conn(facet.element, node)); this->ASR_nodes(facet_conn(facet.element, node)) = true; } } this->asr_central_nodes.push_back(cent_node); // store the asr facets // this->ASR_facets_from_mesh_facets.push_back(facets_in_loop); already_inserted++; } std::cout << "Proc " << prank << " placed " << already_inserted << " ASR site" << std::endl; return already_inserted; } /* ------------------------------------------------------------------- */ template UInt ASRTools::insertCohesiveLoops(UInt nb_insertions, std::string mat_name) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); const GhostType gt = _not_ghost; auto element_type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(element_type); auto material_id = model.getMaterialIndex(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 & node_conn = mesh_facets.getConnectivity(_point_1, gt); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); const Real max_dot{0.7}; // instantiate crack_numbers array for (auto && type_facet : mesh_facets.elementTypes(dim - 1)) { ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.getDataPointer("crack_numbers", type_cohesive); } CSR nodes_to_facets; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_facets, dim - 1); std::set matrix_nodes; for_each_element( mesh_facets, [&](auto && point_el) { // get the node number auto node = node_conn(point_el.element); // discard ghost nodes if (mesh.isPureGhostNode(node)) goto nextnode; for (auto & facet : nodes_to_facets.getRow(node)) { auto & facet_material = coh_model.getFacetMaterial( facet.type, facet.ghost_type)(facet.element); if (facet_material != material_id) { goto nextnode; } } matrix_nodes.emplace(node); nextnode:; }, _spatial_dimension = 0, _ghost_type = _not_ghost); UInt nb_nodes = matrix_nodes.size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_nodes - 1); UInt already_inserted = 0; std::set left_nodes = matrix_nodes; Array> loops; while (already_inserted < nb_insertions and nb_nodes and left_nodes.size()) { auto id = dis(random_generator); auto it = matrix_nodes.begin(); std::advance(it, id); UInt cent_node(*it); left_nodes.erase(*it); // not touching other cracks if (this->ASR_nodes(cent_node)) continue; // поехали CSR segments_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, segments_to_nodes, dim - 2); auto && segments_to_node = segments_to_nodes.getRow(cent_node); Array segments_list; for (auto && segment : segments_to_node) { segments_list.push_back(segment); } // pick the first non-ghost segment in the list Element starting_segment; for (auto segment : segments_list) { if (segment.ghost_type == _not_ghost) { starting_segment = segment; break; } } // go to next node if starting segment could not be picked if (starting_segment == ElementNull) { std::cout << "Could not pick the starting segment, switching to the next node" << std::endl; continue; } auto & facets_to_segment = mesh_facets.getElementToSubelement(starting_segment); // pick the first good facet in the list Element starting_facet{ElementNull}; for (const auto & facet : facets_to_segment) { // check if the facet and its nodes are ok if (isFacetAndNodesGood(facet, material_id, true)) { starting_facet = facet; break; } } // go to next node if starting facet could not be picked if (starting_facet == ElementNull) { std::cout << "Could not pick the starting facet, switching to the next node" << std::endl; continue; } // call the actual function that build the bridge auto facets_in_loop = findFacetsLoopFromSegment2Segment( starting_facet, starting_segment, starting_segment, cent_node, max_dot, material_id, true); // loop broke -> try another point if (not facets_in_loop.size()) { std::cout << "Couldn't close the facet loop, taking the next node" << std::endl; continue; } // otherwise champaigne loops.push_back(facets_in_loop); for (auto & facet : facets_in_loop) { /// add all facet nodes to the group for (auto node : arange(nb_nodes_facet)) { this->ASR_nodes(facet_conn(facet.element, node)) = true; } } this->asr_central_nodes.push_back(cent_node); already_inserted++; } + /// shift crack numbers by the number of cracks on lower processors + auto num_cracks = loops.size(); + comm.exclusiveScan(num_cracks); + // convert loops into a map std::map facet_nbs_crack_nbs; for (auto && data : enumerate(loops)) { - auto crack_nb = std::get<0>(data); + auto crack_nb = std::get<0>(data) + num_cracks; auto && facet_loop = std::get<1>(data); for (auto && facet : facet_loop) { facet_nbs_crack_nbs[facet.element] = crack_nb; } } // insert cohesives in a corresponding material auto & mat = model.getMaterial(material_id); auto * mat_coh = dynamic_cast *>(&mat); AKANTU_DEBUG_ASSERT(mat_coh != nullptr, "Chosen material is not linear sequential"); mat_coh->insertCohesiveElements(facet_nbs_crack_nbs, facet_type, false); inserter.insertElements(); // fill up ASR facets vector this->asr_facets.resize(loops.size()); for (auto && data : zip(loops, this->asr_facets)) { auto && facet_loop = std::get<0>(data); auto && asr_facet_loop = std::get<1>(data); for (auto & facet : facet_loop) { Element cohesive{ElementNull}; auto & connected_els = mesh_facets.getElementToSubelement(facet); for (auto & connected_el : connected_els) { if (mesh.getKind(connected_el.type) == _ek_cohesive) { cohesive = connected_el; break; } } AKANTU_DEBUG_ASSERT(cohesive != ElementNull, "Connected cohesive element not identified"); auto && connected_subels = mesh_facets.getSubelementToElement(cohesive); for (auto & connected_subel : connected_subels) { if (connected_subel.type != facet.type) continue; asr_facet_loop.push_back(connected_subel); } } } + // split nodes into pairs and compute averaged normals + identifyASRCentralNodePairsAndNormals(); + std::cout << "Proc " << prank << " placed " << already_inserted << " ASR site" << std::endl; return already_inserted; } template UInt ASRTools::insertCohesiveLoops<2>(UInt nb_insertions, std::string mat_name); template UInt ASRTools::insertCohesiveLoops<3>(UInt nb_insertions, std::string mat_name); /* ------------------------------------------------------------------- */ Array ASRTools::findFacetsLoopFromSegment2Segment( Element directional_facet, Element starting_segment, Element ending_segment, UInt cent_node, Real max_dot, UInt material_id, bool check_asr_nodes) { auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT( dim == 3, "findFacetsLoopFromSegment2Segment currently supports only 3D"); AKANTU_DEBUG_ASSERT(Mesh::getSpatialDimension(directional_facet.type) == dim - 1, "Directional facet provided has wrong spatial dimension"); if (Mesh::getSpatialDimension(starting_segment.type) != dim - 2 or Mesh::getSpatialDimension(ending_segment.type) != dim - 2) { AKANTU_EXCEPTION("Starting facet provided has wrong spatial dimension"); } // get the point el corresponding to the central node CSR points_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, points_to_nodes, 0); auto && points_to_node = points_to_nodes.getRow(cent_node); auto point_el = *points_to_node.begin(); AKANTU_DEBUG_ASSERT(point_el.type == _point_1, "Provided node element type is wrong"); // check if the provided node is shared by two segments auto & segments_to_node = mesh_facets.getElementToSubelement(point_el); auto ret1 = std::find(segments_to_node.begin(), segments_to_node.end(), starting_segment); AKANTU_DEBUG_ASSERT(ret1 != segments_to_node.end(), "Starting segment doesn't contain the central node"); auto ret2 = std::find(segments_to_node.begin(), segments_to_node.end(), ending_segment); AKANTU_DEBUG_ASSERT(ret2 != segments_to_node.end(), "Ending segment doesn't contain the central node"); // check if the directional facet contain starting segment auto & facets_to_segment = mesh_facets.getElementToSubelement(starting_segment); auto ret3 = std::find(facets_to_segment.begin(), facets_to_segment.end(), directional_facet); AKANTU_DEBUG_ASSERT(ret3 != facets_to_segment.end(), "Starting facet doesn't contain the starting segment"); // loop through facets to arrive from the starting segment to the endin bool loop_closed{false}; auto current_facet(directional_facet); auto border_segment(starting_segment); Array facets_in_loop; Array segments_in_loop; while (not loop_closed and facets_in_loop.size() < 10) { // loop through all the neighboring facets and pick the prettiest Element next_facet(ElementNull); Real current_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, current_facet); bool almost_closed{false}; auto neighbor_facets = mesh_facets.getElementToSubelement(border_segment); for (auto neighbor_facet : neighbor_facets) { if (neighbor_facet == current_facet) continue; if (not isFacetAndNodesGood(neighbor_facet, material_id, check_asr_nodes)) continue; // first check if it has the ending segment -> loop closed bool ending_segment_touched{false}; if (facets_in_loop.size()) { const Vector neighbors_segments = mesh_facets.getSubelementToElement(neighbor_facet); auto ret = std::find(neighbors_segments.begin(), neighbors_segments.end(), ending_segment); if (ret != neighbors_segments.end()) { ending_segment_touched = true; almost_closed = true; } } // conditions on the angle between facets Real neighbor_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, neighbor_facet); Real hypotenuse = sqrt(neighbor_facet_indiam * neighbor_facet_indiam + current_facet_indiam * current_facet_indiam) / 2; auto dist = MeshUtils::distanceBetweenIncentersCorrected( mesh_facets, current_facet, neighbor_facet); if (dist < hypotenuse) { continue; } // get abs of dot product between two normals Real dot = MeshUtils::cosSharpAngleBetween2Facets(model, current_facet, neighbor_facet); if (dot > max_dot) { next_facet = neighbor_facet; if (ending_segment_touched) { facets_in_loop.push_back(neighbor_facet); loop_closed = true; break; } } } if (next_facet == ElementNull) { facets_in_loop.clear(); return facets_in_loop; } if (almost_closed & not loop_closed) { facets_in_loop.clear(); return facets_in_loop; } if (loop_closed) { return facets_in_loop; } // otherwise make this facet a current one current_facet = next_facet; facets_in_loop.push_back(current_facet); segments_in_loop.push_back(border_segment); // identify the next border segment auto previous_border_segment = border_segment; border_segment = ElementNull; const Vector current_segments = mesh_facets.getSubelementToElement(current_facet); for (auto & segment : current_segments) { if (segment == previous_border_segment) continue; // check if the segment includes the central node auto && nodes_to_segment = mesh_facets.getConnectivity(segment); bool includes_cent_node{false}; for (auto & node : nodes_to_segment) { if (node == cent_node) { includes_cent_node = true; break; } } if (not includes_cent_node) continue; // check if the segment was previously added if (segments_in_loop.find(segment) != UInt(-1)) continue; border_segment = segment; break; } if (border_segment == ElementNull) { facets_in_loop.clear(); return facets_in_loop; } } return facets_in_loop; } /* ------------------------------------------------------------------- */ Array ASRTools::findFacetsLoopByGraphByDist(const Array & limit_facets, const Array & limit_segments, const UInt & cent_node) { auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); Array facets_in_loop; auto && starting_facet = limit_facets(0); auto && ending_facet = limit_facets(1); auto && starting_segment = limit_segments(0); auto && ending_segment = limit_segments(1); Real max_dist = std::numeric_limits::min(); // Create a struct to hold properties for each vertex typedef struct vertex_properties { Element segment; } vertex_properties_t; // Create a struct to hold properties for each edge typedef struct edge_properties { Element facet; Real weight; } edge_properties_t; // Define the type of the graph typedef boost::adjacency_list graph_t; typedef graph_t::vertex_descriptor vertex_descriptor_t; typedef graph_t::edge_descriptor edge_descriptor_t; typedef boost::graph_traits::edge_iterator edge_iter; typedef boost::graph_traits::vertex_iterator vertex_iter; graph_t g; // prepare the list of facets sharing the central node CSR facets_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, facets_to_nodes, dim - 1); auto && facets_to_node = facets_to_nodes.getRow(cent_node); std::set facets_added; std::map segment_vertex_map; // get the point el corresponding to the central node CSR points_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, points_to_nodes, 0); auto && points_to_node = points_to_nodes.getRow(cent_node); auto point_el = *points_to_node.begin(); // get list of all segments connected to the node auto && segm_to_point = mesh_facets.getElementToSubelement(point_el); // add starting and ending segment to the graph vertex_descriptor_t v_start = boost::add_vertex(g); g[v_start].segment = starting_segment; vertex_descriptor_t v_end = boost::add_vertex(g); g[v_end].segment = ending_segment; segment_vertex_map[starting_segment] = v_start; segment_vertex_map[ending_segment] = v_end; for (auto && facet : facets_to_node) { // discard undesired facets first if (facet == starting_facet or facet == ending_facet) continue; if (not isFacetAndNodesGood(facet)) continue; if (belong2SameElement(starting_facet, facet) or belong2SameElement(ending_facet, facet)) continue; // identify 2 side facets Array side_segments; auto && facet_segments = mesh_facets.getSubelementToElement(facet); for (auto & segment : facet_segments) { auto ret = std::find(segm_to_point.begin(), segm_to_point.end(), segment); if (ret != segm_to_point.end()) { side_segments.push_back(segment); } } AKANTU_DEBUG_ASSERT(side_segments.size() == 2, "Number of side facets is wrong"); // add corresponding vertex (if not there)and edge into graph Array v(2); for (UInt i = 0; i < 2; i++) { auto it = segment_vertex_map.find(side_segments(i)); if (it != segment_vertex_map.end()) { v(i) = it->second; } else { v(i) = boost::add_vertex(g); g[v(i)].segment = side_segments(i); segment_vertex_map[side_segments(i)] = v(i); } } // add corresponding edge into the graph auto ret = facets_added.emplace(facet); if (ret.second) { // compute facet's area to use as a weight // auto facet_area = MeshUtils::getFacetArea(model, facet); auto dist = MeshUtils::distanceBetweenIncenters(mesh_facets, starting_facet, facet); dist += MeshUtils::distanceBetweenIncenters(mesh_facets, ending_facet, facet); if (dist > max_dist) max_dist = dist; // add edge std::pair e = boost::add_edge(v(0), v(1), g); g[e.first].weight = dist; // g[e.first].weight = facet_area; g[e.first].facet = facet; } } // normalize and flip weights edge_iter ei, ei_end; for (std::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) { auto weight = g[*ei].weight; weight = max_dist - weight; // weight = std::pow(weight, 4);//give less weight for smaller values auto facet_area = MeshUtils::getFacetArea(model, g[*ei].facet); g[*ei].weight = weight + sqrt(facet_area); } // Print out some useful information std::cout << "Starting segm " << starting_segment.element << " ending segm " << ending_segment.element << std::endl; std::cout << "Graph:" << std::endl; for (std::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) { std::cout << g[*ei].facet.element << " "; } std::cout << std::endl; std::cout << "num_verts: " << boost::num_vertices(g) << std::endl; std::cout << "num_edges: " << boost::num_edges(g) << std::endl; // BGL Dijkstra's Shortest Paths here... std::vector distances(boost::num_vertices(g)); std::vector predecessors(boost::num_vertices(g)); boost::dijkstra_shortest_paths( g, v_start, boost::weight_map(boost::get(&edge_properties_t::weight, g)) .distance_map(boost::make_iterator_property_map( distances.begin(), boost::get(boost::vertex_index, g))) .predecessor_map(boost::make_iterator_property_map( predecessors.begin(), boost::get(boost::vertex_index, g)))); std::cout << "distances and parents:" << std::endl; vertex_iter vi, vend; for (boost::tie(vi, vend) = boost::vertices(g); vi != vend; ++vi) { std::cout << "distance(" << g[*vi].segment.element << ") = " << distances[*vi] << ", "; std::cout << "parent = " << g[predecessors[*vi]].segment.element << std::endl; } // Extract the shortest path from v1 to v3. typedef std::vector path_t; path_t path; vertex_descriptor_t v = v_end; for (vertex_descriptor_t u = predecessors[v]; u != v; v = u, u = predecessors[v]) { std::pair edge_pair = boost::edge(u, v, g); path.push_back(edge_pair.first); } std::cout << std::endl; std::cout << "Shortest Path from segment " << starting_segment.element << " to " << ending_segment.element << ":" << std::endl; for (path_t::reverse_iterator riter = path.rbegin(); riter != path.rend(); ++riter) { vertex_descriptor_t u_tmp = boost::source(*riter, g); vertex_descriptor_t v_tmp = boost::target(*riter, g); edge_descriptor_t e_tmp = boost::edge(u_tmp, v_tmp, g).first; std::cout << " " << g[u_tmp].segment.element << " -> " << g[v_tmp].segment.element << " through " << g[e_tmp].facet.element << " (area: " << g[e_tmp].weight << ")" << std::endl; facets_in_loop.push_back(g[e_tmp].facet); } return facets_in_loop; } /* ------------------------------------------------------------------- */ Array ASRTools::findFacetsLoopByGraphByArea( const Array & limit_facets, const Array & limit_segments, const Array & preinserted_facets, const UInt & cent_node) { auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); Array facets_in_loop; // find two best neighbors Array best_neighbors(2, 1, ElementNull); Array neighbors_borders(2, 1, ElementNull); for (UInt i = 0; i < 2; i++) { Element best_neighbor{ElementNull}; // if the facet is already provided - use it if (preinserted_facets(i) != ElementNull) { best_neighbor = preinserted_facets(i); } else { // if not - search for it auto limit_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, limit_facets(i)); auto max_dot = std::numeric_limits::min(); auto && neighbor_facets = mesh_facets.getElementToSubelement(limit_segments(i)); for (auto && neighbor_facet : neighbor_facets) { if (neighbor_facet == limit_facets(0) or neighbor_facet == limit_facets(1)) continue; if (not isFacetAndNodesGood(neighbor_facet)) continue; // conditions on the angle between facets Real neighbor_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, neighbor_facet); Real hypotenuse = sqrt(neighbor_facet_indiam * neighbor_facet_indiam + limit_facet_indiam * limit_facet_indiam) / 2; auto dist = MeshUtils::distanceBetweenIncentersCorrected( mesh_facets, limit_facets(i), neighbor_facet); if (dist < hypotenuse) continue; // find neighbor with the biggest thing // get abs of dot product between two normals Real dot = MeshUtils::cosSharpAngleBetween2Facets( model, limit_facets(i), neighbor_facet); if (dot > max_dot) { max_dot = dot; best_neighbor = neighbor_facet; } } } if (best_neighbor == ElementNull) { return facets_in_loop; } else { best_neighbors(i) = best_neighbor; } // find new limiting segments auto && neighbor_segments = mesh_facets.getSubelementToElement(best_neighbor); for (auto & segment : neighbor_segments) { if (segment == limit_segments(i)) continue; // check if the segment includes the central node auto && nodes_to_segment = mesh_facets.getConnectivity(segment); bool includes_cent_node{false}; for (auto & node : nodes_to_segment) { if (node == cent_node) { includes_cent_node = true; break; } } if (not includes_cent_node) continue; neighbors_borders(i) = segment; break; } } // check if borders are not the same segment if (neighbors_borders(0) == neighbors_borders(1)) { for (auto && best_neighbor : best_neighbors) { facets_in_loop.push_back(best_neighbor); } } // Create a struct to hold properties for each vertex typedef struct vertex_properties { Element segment; } vertex_properties_t; // Create a struct to hold properties for each edge typedef struct edge_properties { Element facet; Real weight; } edge_properties_t; // Define the type of the graph typedef boost::adjacency_list graph_t; typedef graph_t::vertex_descriptor vertex_descriptor_t; typedef graph_t::edge_descriptor edge_descriptor_t; // typedef boost::graph_traits::edge_iterator edge_iter; // typedef boost::graph_traits::vertex_iterator vertex_iter; graph_t g; // prepare the list of facets sharing the central node CSR facets_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, facets_to_nodes, dim - 1); auto && facets_to_node = facets_to_nodes.getRow(cent_node); std::set facets_added; std::map segment_vertex_map; // get the point el corresponding to the central node CSR points_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, points_to_nodes, 0); auto && points_to_node = points_to_nodes.getRow(cent_node); auto point_el = *points_to_node.begin(); // get list of all segments connected to the node auto && segm_to_point = mesh_facets.getElementToSubelement(point_el); // add starting and ending segment to the graph vertex_descriptor_t v_start = boost::add_vertex(g); g[v_start].segment = neighbors_borders(0); vertex_descriptor_t v_end = boost::add_vertex(g); g[v_end].segment = neighbors_borders(1); segment_vertex_map[neighbors_borders(0)] = v_start; segment_vertex_map[neighbors_borders(1)] = v_end; // compute inscribed circles diameters Array best_neighbors_indiam; for (auto && best_neighbor : best_neighbors) { auto neighbor_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, best_neighbor); best_neighbors_indiam.push_back(neighbor_facet_indiam); } for (auto && facet : facets_to_node) { // discard undesired facets first if (facet == best_neighbors(0) or facet == best_neighbors(1) or facet == limit_facets(0) or facet == limit_facets(1)) continue; if (not isFacetAndNodesGood(facet)) continue; if (belong2SameElement(best_neighbors(0), facet) or belong2SameElement(best_neighbors(1), facet) or belong2SameElement(limit_facets(0), facet) or belong2SameElement(limit_facets(1), facet)) continue; // discard short distances to starting and ending facets bool short_dist{false}; auto facet_indiam = MeshUtils::getInscribedCircleDiameter(model, facet); for (UInt i = 0; i < 2; i++) { auto hypotenuse = sqrt(best_neighbors_indiam(i) * best_neighbors_indiam(i) + facet_indiam * facet_indiam) / 2; auto dist = MeshUtils::distanceBetweenIncenters(mesh_facets, best_neighbors(i), facet); if (dist <= hypotenuse) { short_dist = true; break; } } if (short_dist) continue; // identify 2 side facets Array side_segments; auto && facet_segments = mesh_facets.getSubelementToElement(facet); for (auto & segment : facet_segments) { auto ret = std::find(segm_to_point.begin(), segm_to_point.end(), segment); if (ret != segm_to_point.end()) { side_segments.push_back(segment); } } AKANTU_DEBUG_ASSERT(side_segments.size() == 2, "Number of side facets is wrong"); // add corresponding vertex (if not there)and edge into graph Array v(2); for (UInt i = 0; i < 2; i++) { auto it = segment_vertex_map.find(side_segments(i)); if (it != segment_vertex_map.end()) { v(i) = it->second; } else { v(i) = boost::add_vertex(g); g[v(i)].segment = side_segments(i); segment_vertex_map[side_segments(i)] = v(i); } } // add corresponding edge into the graph auto ret = facets_added.emplace(facet); if (ret.second) { // compute facet's area to use as a weight auto facet_area = MeshUtils::getFacetArea(model, facet); // add edge std::pair e = boost::add_edge(v(0), v(1), g); g[e.first].weight = facet_area; g[e.first].facet = facet; } } // // Print out some useful information // std::cout << "Starting facet " << best_neighbors(0).element // << " ending facet " << best_neighbors(0).element << std::endl; // std::cout << "Graph:" << std::endl; // // boost::print_graph(g, boost::get(&vertex_properties_t::segment, g)); // edge_iter ei, ei_end; // for (std::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) { // std::cout << g[*ei].facet.element << ","; // } // std::cout << std::endl; // for (std::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) { // vertex_descriptor_t u = boost::source(*ei, g); // vertex_descriptor_t v = boost::target(*ei, g); // std::cout << g[*ei].facet.element << " from " << g[u].segment.element // << "->" << g[v].segment.element << " area " << g[*ei].weight // << std::endl; // } // std::cout << "num_verts: " << boost::num_vertices(g) << std::endl; // std::cout << "num_edges: " << boost::num_edges(g) << std::endl; // std::cout << std::endl; // BGL Dijkstra's Shortest Paths here... std::vector distances(boost::num_vertices(g)); std::vector predecessors(boost::num_vertices(g)); boost::dijkstra_shortest_paths( g, v_end, boost::weight_map(boost::get(&edge_properties_t::weight, g)) .distance_map(boost::make_iterator_property_map( distances.begin(), boost::get(boost::vertex_index, g))) .predecessor_map(boost::make_iterator_property_map( predecessors.begin(), boost::get(boost::vertex_index, g)))); // Extract the shortest path from v1 to v3. typedef std::vector path_t; path_t path; vertex_descriptor_t v = v_start; for (vertex_descriptor_t u = predecessors[v]; u != v; v = u, u = predecessors[v]) { std::pair edge_pair = boost::edge(u, v, g); path.push_back(edge_pair.first); } if (path.size()) { // add 2 determenistic facets for (auto && best_neighbor : best_neighbors) { facets_in_loop.push_back(best_neighbor); } // add the shortest path facets for (path_t::reverse_iterator riter = path.rbegin(); riter != path.rend(); ++riter) { vertex_descriptor_t u_tmp = boost::source(*riter, g); vertex_descriptor_t v_tmp = boost::target(*riter, g); edge_descriptor_t e_tmp = boost::edge(u_tmp, v_tmp, g).first; facets_in_loop.push_back(g[e_tmp].facet); } } return facets_in_loop; } /* ------------------------------------------------------------------- */ Array ASRTools::findStressedFacetLoopAroundNode( const Array & limit_facets, const Array & limit_segments, const UInt & cent_node, Real min_dot) { auto & coh_model = dynamic_cast(model); auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(limit_facets(0) != limit_facets(1), "Limit facets are equal"); AKANTU_DEBUG_ASSERT(limit_segments(0) != limit_segments(1), "Limit segments are equal"); AKANTU_DEBUG_ASSERT( dim == 3, "findFacetsLoopFromSegment2Segment currently supports only 3D"); // get the point el corresponding to the central node CSR points_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, points_to_nodes, 0); auto && points_to_node = points_to_nodes.getRow(cent_node); auto point_el = *points_to_node.begin(); auto & segments_to_point = mesh_facets.getElementToSubelement(point_el); AKANTU_DEBUG_ASSERT(point_el.type == _point_1, "Provided node element type is wrong"); for (UInt i = 0; i < 2; i++) { // check the dimensions AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(limit_facets(i).type) == dim - 1, "Facet" << limit_facets(i) << " has wrong spatial dimension"); AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(limit_segments(i).type) == dim - 2, "Segment" << limit_segments(i) << " has wrong spatial dimension"); // check if the provided node is shared by two segments auto ret = std::find(segments_to_point.begin(), segments_to_point.end(), limit_segments(i)); AKANTU_DEBUG_ASSERT(ret != segments_to_point.end(), "Segment " << limit_segments(i) << " doesn't contain the central node"); // check if limit facets contain limit segments auto & facets_to_segment = mesh_facets.getElementToSubelement(limit_segments(i)); auto ret1 = std::find(facets_to_segment.begin(), facets_to_segment.end(), limit_facets(i)); AKANTU_DEBUG_ASSERT(ret1 != facets_to_segment.end(), "Facet " << limit_facets(i) << " doesn't contain the segment " << limit_segments(i)); } Array facets_in_loop; // Create a struct to hold properties for each vertex typedef struct vertex_properties { Element segment; } vertex_properties_t; // Create a struct to hold properties for each edge typedef struct edge_properties { Element facet; Real weight; } edge_properties_t; // Define the type of the graph typedef boost::adjacency_list graph_t; typedef typename graph_t::vertex_descriptor vertex_descriptor_t; typedef typename graph_t::edge_descriptor edge_descriptor_t; graph_t g; // prepare the list of facets sharing the central node CSR facets_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, facets_to_nodes, dim - 1); auto && facets_to_node = facets_to_nodes.getRow(cent_node); std::set facets_added; std::map segment_vertex_map; // add starting and ending segment to the graph vertex_descriptor_t v_start = boost::add_vertex(g); g[v_start].segment = limit_segments(0); vertex_descriptor_t v_end = boost::add_vertex(g); g[v_end].segment = limit_segments(1); segment_vertex_map[limit_segments(0)] = v_start; segment_vertex_map[limit_segments(1)] = v_end; for (auto && facet : facets_to_node) { // discard undesired facets first if (facet == limit_facets(0) or facet == limit_facets(1)) continue; if (not isFacetAndNodesGood(facet)) continue; if (belong2SameElement(limit_facets(0), facet) or belong2SameElement(limit_facets(1), facet)) continue; // check the effective stress on this facet auto & facet_material_index = coh_model.getFacetMaterial(facet.type, facet.ghost_type)(facet.element); auto & mat = model.getMaterial(facet_material_index); auto * mat_coh = dynamic_cast(&mat); if (mat_coh == nullptr) continue; auto && mat_facet_filter = mat_coh->getFacetFilter(facet.type, facet.ghost_type); UInt facet_local_num = mat_facet_filter.find(facet.element); AKANTU_DEBUG_ASSERT(facet_local_num != UInt(-1), "Local facet number was not identified"); if (not mat_coh->isInternal("effective_stresses", _ek_regular)) continue; auto & eff_stress_array = mat.getInternal("effective_stresses")( facet.type, facet.ghost_type); // eff stress is always non-negative auto eff_stress = eff_stress_array(facet_local_num); if (eff_stress == 0) eff_stress = 1e-3; // if (eff_stress < 1) // continue; // identify 2 sides of a facets Array side_segments; auto && facet_segments = mesh_facets.getSubelementToElement(facet); for (auto & segment : facet_segments) { auto ret = std::find(segments_to_point.begin(), segments_to_point.end(), segment); if (ret != segments_to_point.end()) { side_segments.push_back(segment); } } AKANTU_DEBUG_ASSERT(side_segments.size() == 2, "Number of side facets is wrong"); // add corresponding vertex (if not there)and edge into graph Array v(2); for (UInt i = 0; i < 2; i++) { auto it = segment_vertex_map.find(side_segments(i)); if (it != segment_vertex_map.end()) { v(i) = it->second; } else { v(i) = boost::add_vertex(g); g[v(i)].segment = side_segments(i); segment_vertex_map[side_segments(i)] = v(i); } } // add corresponding edge into the graph auto ret = facets_added.emplace(facet); if (ret.second) { // compute facet's area to use as a weight // auto facet_area = MeshUtils::getFacetArea(model, facet); // add edge std::pair e = boost::add_edge(v(0), v(1), g); // g[e.first].weight = facet_area; g[e.first].weight = 1 / eff_stress; g[e.first].facet = facet; } } // BGL Dijkstra's Shortest Paths here... std::vector distances(boost::num_vertices(g)); std::vector predecessors(boost::num_vertices(g)); boost::dijkstra_shortest_paths( g, v_start, boost::weight_map(boost::get(&edge_properties_t::weight, g)) .distance_map(boost::make_iterator_property_map( distances.begin(), boost::get(boost::vertex_index, g))) .predecessor_map(boost::make_iterator_property_map( predecessors.begin(), boost::get(boost::vertex_index, g)))); // Extract the shortest path from v_end to v_start (reversed order) typedef std::vector path_t; path_t path; vertex_descriptor_t v = v_end; for (vertex_descriptor_t u = predecessors[v]; u != v; v = u, u = predecessors[v]) { std::pair edge_pair = boost::edge(u, v, g); path.push_back(edge_pair.first); } // compute average effective stress over the loop and discard if <1 Array facet_loop_tmp; for (auto edge : path) { facet_loop_tmp.push_back(g[edge].facet); } auto av_eff_stress = averageEffectiveStressInMultipleFacets(facet_loop_tmp); if (av_eff_stress < 1) return facets_in_loop; // check path for smoothness bool recompute_path{false}; do { recompute_path = false; if (path.size()) { for (UInt i = 0; i <= path.size(); i++) { // evaluate two neighboring facets Element facet1, facet2; if (i == 0) { facet1 = limit_facets(1); facet2 = g[path[i]].facet; } else if (i == path.size()) { facet1 = g[path[i - 1]].facet; facet2 = limit_facets(0); } else { facet1 = g[path[i - 1]].facet; facet2 = g[path[i]].facet; } // discard short distances auto facet1_indiam = MeshUtils::getInscribedCircleDiameter(model, facet1); auto facet2_indiam = MeshUtils::getInscribedCircleDiameter(model, facet2); auto hypotenuse = sqrt(facet1_indiam * facet1_indiam + facet2_indiam * facet2_indiam) / 2; auto dist = MeshUtils::distanceBetweenIncentersCorrected( mesh_facets, facet1, facet2); // get abs of dot product between two normals Real dot = MeshUtils::cosSharpAngleBetween2Facets(model, facet1, facet2); if (dist < hypotenuse or dot < min_dot) { if (i < path.size()) { boost::remove_edge(path[i], g); } else { boost::remove_edge(path[i - 1], g); } recompute_path = true; break; } } // find new shortest path without problematic edge if (recompute_path) { boost::dijkstra_shortest_paths( g, v_start, boost::weight_map(boost::get(&edge_properties_t::weight, g)) .distance_map(boost::make_iterator_property_map( distances.begin(), boost::get(boost::vertex_index, g))) .predecessor_map(boost::make_iterator_property_map( predecessors.begin(), boost::get(boost::vertex_index, g)))); path.clear(); v = v_end; for (vertex_descriptor_t u = predecessors[v]; u != v; v = u, u = predecessors[v]) { std::pair edge_pair = boost::edge(u, v, g); path.push_back(edge_pair.first); } } } } while (recompute_path); // add the shortest path facets for (typename path_t::reverse_iterator riter = path.rbegin(); riter != path.rend(); ++riter) { vertex_descriptor_t u_tmp = boost::source(*riter, g); vertex_descriptor_t v_tmp = boost::target(*riter, g); edge_descriptor_t e_tmp = boost::edge(u_tmp, v_tmp, g).first; facets_in_loop.push_back(g[e_tmp].facet); } return facets_in_loop; } /* ------------------------------------------------------------------- */ Array ASRTools::findSingleFacetLoop(const Array & limit_facets, const Array & limit_segments, const UInt & cent_node) { auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(limit_facets(0) != limit_facets(1), "Limit facets are equal"); AKANTU_DEBUG_ASSERT(limit_segments(0) != limit_segments(1), "Limit segments are equal"); AKANTU_DEBUG_ASSERT( dim == 3, "findFacetsLoopFromSegment2Segment currently supports only 3D"); // get the point el corresponding to the central node CSR points_to_nodes; MeshUtils::buildNode2Elements(mesh_facets, points_to_nodes, 0); auto && points_to_node = points_to_nodes.getRow(cent_node); auto point_el = *points_to_node.begin(); auto & segments_to_node = mesh_facets.getElementToSubelement(point_el); AKANTU_DEBUG_ASSERT(point_el.type == _point_1, "Provided node element type is wrong"); for (UInt i = 0; i < 2; i++) { // check the dimensions AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(limit_facets(i).type) == dim - 1, "Facet" << limit_facets(i) << " has wrong spatial dimension"); AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(limit_segments(i).type) == dim - 2, "Segment" << limit_segments(i) << " has wrong spatial dimension"); // check if the provided node is shared by two segments auto ret = std::find(segments_to_node.begin(), segments_to_node.end(), limit_segments(i)); AKANTU_DEBUG_ASSERT(ret != segments_to_node.end(), "Segment " << limit_segments(i) << " doesn't contain the central node"); // check if limit facets contain limit segments auto & facets_to_segment = mesh_facets.getElementToSubelement(limit_segments(i)); auto ret1 = std::find(facets_to_segment.begin(), facets_to_segment.end(), limit_facets(i)); AKANTU_DEBUG_ASSERT(ret1 != facets_to_segment.end(), "Facet " << limit_facets(i) << " doesn't contain the segment " << limit_segments(i)); } Array facets_in_loop; // check if two segments are not connected by a single facet auto st_segment_neighbors = mesh_facets.getElementToSubelement(limit_segments(0)); auto end_segment_neighbors = mesh_facets.getElementToSubelement(limit_segments(1)); std::sort(st_segment_neighbors.begin(), st_segment_neighbors.end()); std::sort(end_segment_neighbors.begin(), end_segment_neighbors.end()); Array shared_facets(st_segment_neighbors.size()); // An iterator to the end of the constructed range. auto it = std::set_intersection( st_segment_neighbors.begin(), st_segment_neighbors.end(), end_segment_neighbors.begin(), end_segment_neighbors.end(), shared_facets.begin()); shared_facets.resize(it - shared_facets.begin()); // if there is a single facet connecting two segments test it if (shared_facets.size() == 1) { bool bad_facet{false}; if (not isFacetAndNodesGood(shared_facets(0))) bad_facet = true; // check angle with starting and ending -> not to be sharp bool sharp_angle{false}; auto facet_indiam = MeshUtils::getInscribedCircleDiameter(model, shared_facets(0)); for (auto && limit_facet : limit_facets) { auto limit_facet_indiam = MeshUtils::getInscribedCircleDiameter(model, limit_facet); auto hypotenuse = sqrt(limit_facet_indiam * limit_facet_indiam + facet_indiam * facet_indiam) / 2; auto dist = MeshUtils::distanceBetweenIncentersCorrected( mesh_facets, limit_facet, shared_facets(0)); if (dist <= hypotenuse) { sharp_angle = true; break; } } if (not sharp_angle and not bad_facet) { facets_in_loop.push_back(shared_facets(0)); } } return facets_in_loop; } /* ------------------------------------------------------------------- */ bool ASRTools::isFacetAndNodesGood(const Element & facet, UInt material_id, bool check_asr_nodes) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto facet_type = facet.type; auto gt = facet.ghost_type; auto & mesh_facets = inserter.getMeshFacets(); Vector facet_nodes = mesh_facets.getConnectivity(facet); auto & check_facets = inserter.getCheckFacets(facet_type, gt); // check if the facet is not ghost if (gt == _ghost) return false; // check if the facet is allowed to be inserted if (check_facets(facet.element) == false) return false; for (auto & facet_node : facet_nodes) { // check if the facet's nodes are not pure ghost nodes if (model.getMesh().isPureGhostNode(facet_node)) return false; // check if the facet's nodes are asr nodes if (check_asr_nodes and this->ASR_nodes(facet_node)) return false; } // check if the facet material is good if (material_id != UInt(-1)) { const auto & facet_material = coh_model.getFacetMaterial(facet_type, gt)(facet.element); if (facet_material != material_id) return false; } return true; } /* ------------------------------------------------------------------- */ bool ASRTools::belong2SameElement(const Element & facet1, const Element & facet2) { auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); Array f1_neighbors, f2_neighbors; auto facet1_neighbors = mesh_facets.getElementToSubelement(facet1); for (auto && neighbor : facet1_neighbors) { if (mesh.getKind(neighbor.type) == _ek_regular) f1_neighbors.push_back(neighbor); if (mesh.getKind(neighbor.type) == _ek_cohesive) { auto && coh_neighbors = mesh_facets.getSubelementToElement(neighbor); for (auto && coh_neighbor : coh_neighbors) { if (coh_neighbor == facet1) continue; auto coh_facet_neighbors = mesh_facets.getElementToSubelement(coh_neighbor); for (auto && coh_facet_neighbor : coh_facet_neighbors) { if (mesh.getKind(coh_facet_neighbor.type) == _ek_regular) f1_neighbors.push_back(coh_facet_neighbor); } } } } auto facet2_neighbors = mesh_facets.getElementToSubelement(facet2); for (auto && neighbor : facet2_neighbors) { if (mesh.getKind(neighbor.type) == _ek_regular) f2_neighbors.push_back(neighbor); if (mesh.getKind(neighbor.type) == _ek_cohesive) { auto && coh_neighbors = mesh_facets.getSubelementToElement(neighbor); for (auto && coh_neighbor : coh_neighbors) { if (coh_neighbor == facet2) continue; auto coh_facet_neighbors = mesh_facets.getElementToSubelement(coh_neighbor); for (auto && coh_facet_neighbor : coh_facet_neighbors) { if (mesh.getKind(coh_facet_neighbor.type) == _ek_regular) f2_neighbors.push_back(coh_facet_neighbor); } } } } std::sort(f1_neighbors.begin(), f1_neighbors.end()); std::sort(f2_neighbors.begin(), f2_neighbors.end()); Array shared_neighbors(f1_neighbors.size()); // An iterator to the end of the constructed range. auto it = std::set_intersection(f1_neighbors.begin(), f1_neighbors.end(), f2_neighbors.begin(), f2_neighbors.end(), shared_neighbors.begin()); shared_neighbors.resize(it - shared_neighbors.begin()); if (shared_neighbors.size() == 0) return false; // check kind of shared neighbor auto shared_neighbor = shared_neighbors(0); if (mesh.getKind(shared_neighbor.type) != _ek_regular) return false; return true; } /* ------------------------------------------------------------------- */ bool ASRTools::isNodeWithinMaterial(Element & node, UInt material_id) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = coh_model.getSpatialDimension(); CSR nodes_to_facets; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_facets, dim - 1); for (auto & facet : nodes_to_facets.getRow(node.element)) { auto & facet_material = coh_model.getFacetMaterial(facet.type, facet.ghost_type)(facet.element); if (facet_material != material_id) return false; } return true; } /* ------------------------------------------------------------------- */ bool ASRTools::pickFacetNeighborsOld(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); two_neighbors.set(ElementNull); 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(); // 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 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; two_neighbors(node) = elem; } } } // insert neighbors only if two of them were identified if (two_neighbors.find(ElementNull) == 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; } /* ----------------------------------------------------------------------- */ 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(); auto facet_type = cent_facet.type; auto facet_gt = cent_facet.ghost_type; auto facet_nb = cent_facet.element; auto & doubled_facets = mesh_facets.getElementGroup("doubled_facets"); auto & doubled_nodes = mesh.getNodeGroup("doubled_nodes"); auto & facet_conn = mesh_facets.getConnectivity(facet_type, facet_gt); auto & facet_material_index = coh_model.getFacetMaterial(facet_type, facet_gt)(facet_nb); // get list of the subel connected to the facet (nodes in 2D, segments // in 3D) const Vector & subelements_to_element = mesh_facets.getSubelementToElement(cent_facet); Array neighbors(subelements_to_element.size()); neighbors.set(ElementNull); // max weight parameter computed for each subelement iteration // here this parameter = dot product between two normals // minimum 60 deg between normals -> 120 between dips Real weight_parameter_threshold(0.5); std::vector max_weight_parameters(subelements_to_element.size(), weight_parameter_threshold); // identify a facet's neighbor through each subelement for (UInt i : arange(subelements_to_element.size())) { auto & subel = subelements_to_element(i); auto & connected_elements = mesh_facets.getElementToSubelement( subel.type, subel.ghost_type)(subel.element); for (auto & connected_element : connected_elements) { // check all possible unsatisfactory conditions if (connected_element.type != facet_type) continue; if (connected_element.element == facet_nb) continue; auto & candidate_facet_material_index = coh_model.getFacetMaterial( connected_element.type, connected_element.ghost_type)(connected_element.element); if (candidate_facet_material_index != facet_material_index) continue; if (not inserter.getCheckFacets(connected_element.type, connected_element.ghost_type)( connected_element.element)) continue; // discard facets intersecting other ASR elements bool ASR_node{false}; for (UInt j : arange(facet_conn.getNbComponent())) { auto node = facet_conn(connected_element.element, j); if (this->ASR_nodes(node)) ASR_node = true; } if (ASR_node) continue; // get inscribed diameter Real facet_indiam = MeshUtils::getInscribedCircleDiameter(model, cent_facet); // get distance between two barycenters auto dist = MeshUtils::distanceBetweenIncentersCorrected( mesh_facets, cent_facet, connected_element); // ad-hoc rule on barycenters spacing // it should discard all elements under sharp angle if (dist < facet_indiam) continue; // get abs of dot product between two normals Real dot = MeshUtils::cosSharpAngleBetween2Facets(model, cent_facet, connected_element); auto weight_parameter = dot; if (weight_parameter > max_weight_parameters[i]) { max_weight_parameters[i] = weight_parameter; neighbors(i) = connected_element; } } } // different insertion procedures for 2D and 3D cases switch (dim) { case 2: { if (neighbors.find(ElementNull) == UInt(-1)) { for (auto & neighbor : neighbors) { doubled_facets.add(neighbor); for (UInt node : arange(2)) { this->ASR_nodes(facet_conn(neighbor.element, node)) = true; } } return true; } else return false; } case 3: { auto max_el_pos = std::max_element(max_weight_parameters.begin(), max_weight_parameters.end()); if (*max_el_pos > weight_parameter_threshold) { auto max_param_index = std::distance(max_weight_parameters.begin(), max_el_pos); auto & neighbor = neighbors(max_param_index); doubled_facets.add(neighbor); for (UInt node : arange(facet_conn.getNbComponent())) { doubled_nodes.add(facet_conn(neighbor.element, node)); this->ASR_nodes(facet_conn(neighbor.element, node)) = true; } return true; } else return false; } default: AKANTU_EXCEPTION("Provided dimension is not supported"); } } /* ------------------------------------------------------------------ */ 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 subfacet_type = Mesh::getFacetType(facet_type); // auto & facet_conn = mesh.getConnectivity(facet_type, gt); // auto facet_conn_it = facet_conn.begin(facet_conn.getNbComponent()); auto & subfacets_to_facets = mesh_facets.getSubelementToElement(facet_type, gt); auto nb_subfacet = subfacets_to_facets.getNbComponent(); auto subf_to_fac_it = subfacets_to_facets.begin(nb_subfacet); // 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; // } // } auto & el_group = mesh_facets.getElementGroup("doubled_facets"); Array element_ids = el_group.getElements(facet_type); for (UInt i : arange(element_ids.size() / 2)) { auto facet1 = element_ids(i); auto facet2 = element_ids(i + element_ids.size() / 2); std::vector subfacets1(nb_subfacet); std::vector subfacets2(nb_subfacet); for (UInt i : arange(nb_subfacet)) { subfacets1[i] = subf_to_fac_it[facet1](i); subfacets2[i] = subf_to_fac_it[facet2](i); } std::sort(subfacets1.begin(), subfacets1.end()); std::sort(subfacets2.begin(), subfacets2.end()); std::vector dif_subfacets(2 * nb_subfacet); std::vector::iterator it; // identify subfacets not shared by two facets it = std::set_difference(subfacets1.begin(), subfacets1.end(), subfacets2.begin(), subfacets2.end(), dif_subfacets.begin()); dif_subfacets.resize(it - dif_subfacets.begin()); // prevent insertion of cohesives in the facets including these subf for (auto & subfacet : dif_subfacets) { auto neighbor_facets = mesh_facets.getElementToSubelement(subfacet); for (auto & neighbor_facet : neighbor_facets) { inserter.getCheckFacets(facet_type, gt)(neighbor_facet.element) = false; } } } } /* ------------------------------------------------------------------- */ void ASRTools::insertOppositeFacetsAndCohesives() { 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; 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 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}; new_facets_event.getList().push_back(new_facet); } } MeshUtils::fillElementToSubElementsData(mesh); mesh.sendEvent(new_facets_event); // fill up ASR facets vector this->ASR_facets_from_mesh.resize(this->asr_facets.size()); for (UInt i = 0; i != this->asr_facets.size(); i++) { auto single_site_facets = this->asr_facets(i); for (auto & facet : single_site_facets) { Element cohesive{ElementNull}; auto & connected_els = mesh_facets.getElementToSubelement(facet); for (auto & connected_el : connected_els) { if (mesh.getKind(connected_el.type) == _ek_cohesive) { cohesive = connected_el; break; } } AKANTU_DEBUG_ASSERT(cohesive != ElementNull, "Connected cohesive element not identified"); auto && connected_subels = mesh.getSubelementToElement(cohesive); for (auto & connected_subel : connected_subels) { if (connected_subel.type != facet.type) continue; this->ASR_facets_from_mesh(i).push_back(connected_subel); } } } // 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); auto & crack_numbers = mesh.getDataPointer("crack_numbers", 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()); // 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 + 1; 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; // An iterator to the end of the constructed range. 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()); // switch (dim) { // case 2: { // // 1 common node between 2 segments // if (common_nodes.size() == 1) { // boost::add_edge(i, j, graph); // } // break; // } // case 3: { // 2 or 3 common nodes between 2 triangles if (common_nodes.size() > 1) { boost::add_edge(i, j, graph); } // break; // } // } } } /// 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())) { crack_numbers(coh_el) = 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 & el_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 : el_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 = el_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::insertGap3D(const Real gap_ratio) { AKANTU_DEBUG_IN(); if (gap_ratio == 0) return; auto & mesh = model.getMesh(); auto & mesh_facets = mesh.getMeshFacets(); MeshAccessor mesh_accessor(mesh); auto & pos2modify = mesh_accessor.getNodes(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto & el_group = mesh_facets.getElementGroup("doubled_facets"); auto & pos = mesh.getNodes(); auto pos_it = make_view(pos, dim).begin(); auto pos2modify_it = make_view(pos2modify, dim).begin(); for (auto & type : el_group.elementTypes(dim - 1)) { if ((dim == 3) and (type != _triangle_6)) AKANTU_EXCEPTION("The only facet type supported in 3D is _triangle_6"); const Array element_ids = el_group.getElements(type, gt); auto && fe_engine = model.getFEEngineBoundary(); auto nb_qpoints_per_facet = fe_engine.getNbIntegrationPoints(type, gt); // identify pairs of facets Array facets_considered; for (UInt i : arange(element_ids.size())) { auto el_id = element_ids(i); Element element{type, el_id, gt}; // if element was already considered -> skip if (facets_considered.find(el_id) != UInt(-1)) continue; Array element_ids_half(element_ids.size() / 2); UInt starting_pos; if (i < element_ids.size() / 2) { starting_pos = 0; } else { starting_pos = element_ids.size() / 2; } for (auto j : arange(element_ids_half.size())) { element_ids_half(j) = element_ids(j + starting_pos); } // find a neighbor const Vector & subelements_to_element = mesh_facets.getSubelementToElement(element); // identify a facet's neighbor through each subelement Element neighbor(ElementNull); Element border(ElementNull); for (UInt i : arange(subelements_to_element.size())) { auto & subel = subelements_to_element(i); auto & connected_elements = mesh_facets.getElementToSubelement(subel); for (auto & connected_element : connected_elements) { // check all possible unsatisfactory conditions if (connected_element.type != type) continue; if (connected_element.element == el_id) continue; // search for this neighbor in the doubled_nodes el group UInt neighbor_pos_in_el_ids = element_ids_half.find(connected_element.element); if (neighbor_pos_in_el_ids == UInt(-1)) continue; // if found -> subelement is the border neighbor = connected_element; border = subel; facets_considered.push_back(el_id); facets_considered.push_back(connected_element.element); break; } if (neighbor != ElementNull) break; } AKANTU_DEBUG_ASSERT(neighbor != ElementNull, "Neighbor for the facet was not identified for the " "purpose of the gap insertion"); // average the normal in between two neighbors const auto & normals_on_quad = fe_engine.getNormalsOnIntegrationPoints(type, gt); auto normals_it = make_view(normals_on_quad, dim).begin(); Vector normal = normals_it[el_id * nb_qpoints_per_facet]; normal += Vector(normals_it[neighbor.element * nb_qpoints_per_facet]); normal /= normal.norm(); // and correct it for the direction if (i < element_ids.size() / 2) { normal *= -1; } // compute segment length AKANTU_DEBUG_ASSERT( border.type == _segment_3, "The only supported segment type in 3D is _segment_3"); auto & segment_conn = mesh_facets.getConnectivity(border.type, border.ghost_type); const UInt nb_nodes_segment = segment_conn.getNbComponent(); auto segment_nodes_it = make_view(segment_conn, nb_nodes_segment).begin(); auto segment_nodes = segment_nodes_it[border.element]; // 2 end nodes and the middle node UInt A, B, middle_node; A = segment_nodes(0); B = segment_nodes(1); middle_node = segment_nodes(2); Vector AB = Vector(pos_it[B]) - Vector(pos_it[A]); Real correction = AB.norm() * gap_ratio / 2; Vector half_opening_vector = normal * correction; Vector node_pos(pos2modify_it[middle_node]); node_pos += half_opening_vector; this->modified_pos(middle_node) = 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 // 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->partition_border_nodes.resize(new_nb_nodes); this->nodes_eff_stress.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) + if (this->asr_central_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]; auto central_nodes = this->asr_central_nodes; 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->asr_central_nodes.push_back(new_node); + break; } } } - this->doubled_nodes_ready = true; + + this->asr_central_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::insertCohesiveElementsOnContour() { AKANTU_DEBUG_IN(); auto & coh_model = dynamic_cast(model); CohesiveElementInserter & inserter = coh_model.getElementInserter(); if (not coh_model.getIsExtrinsic()) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } coh_model.interpolateStress(); const Mesh & mesh_facets = coh_model.getMeshFacets(); UInt nb_new_elements(0); auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); // find crack topology within specific material and insert cohesives for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; auto data = mat_coh->determineCrackSurface(); auto & contour_subfacets_coh_el = std::get<0>(data); auto & surface_subfacets_crack_nb = std::get<1>(data); auto & surface_nodes = std::get<2>(data); auto & contour_nodes = std::get<3>(data); if (not contour_subfacets_coh_el.size()) { continue; } std::map facet_nbs_crack_nbs = mat_coh->findCriticalFacetsOnContour(contour_subfacets_coh_el, surface_subfacets_crack_nb, contour_nodes, surface_nodes); mat_coh->insertCohesiveElements(facet_nbs_crack_nbs, type_facet, false); } // insert the most stressed cohesive element nb_new_elements = inserter.insertElements(); // communicate crack numbers communicateCrackNumbers(); // // loop until no more holes are present // bool new_holes_filled{true}; // while (new_holes_filled) { // new_holes_filled = false; // fill in holes for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; auto data = mat_coh->determineCrackSurface(); auto & contour_subfacets_coh_el = std::get<0>(data); auto & surface_subfacets_crack_nb = std::get<1>(data); auto & surface_nodes = std::get<2>(data); // auto & contour_nodes = std::get<3>(data); if (not contour_subfacets_coh_el.size()) { continue; } std::map facet_nbs_crack_nbs = mat_coh->findHolesOnContour( contour_subfacets_coh_el, surface_subfacets_crack_nb, surface_nodes); // new_holes_filled = facet_nbs_crack_nbs.size(); // comm.allReduce(new_holes_filled, SynchronizerOperation::_lor); mat_coh->insertCohesiveElements(facet_nbs_crack_nbs, type_facet, false); } auto nb_new_holes = inserter.insertElements(); nb_new_elements += nb_new_holes; // communicate crack numbers communicateCrackNumbers(); // } AKANTU_DEBUG_OUT(); return nb_new_elements; } template UInt ASRTools::insertCohesiveElementsOnContour<2>(); template UInt ASRTools::insertCohesiveElementsOnContour<3>(); /* ------------------------------------------------------------------- */ template UInt ASRTools::insertCohesiveElementsOnContourByNodes( Real av_stress_threshold) { AKANTU_DEBUG_IN(); auto & coh_model = dynamic_cast(model); CohesiveElementInserter & inserter = coh_model.getElementInserter(); if (not coh_model.getIsExtrinsic()) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } coh_model.interpolateStress(); const Mesh & mesh_facets = coh_model.getMeshFacets(); UInt nb_new_elements(0); auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); // find global crack topology auto data = determineCrackSurface(); auto & contour_subfacets_coh_el = std::get<0>(data); auto & surface_subfacets_crack_nb = std::get<1>(data); auto & surface_nodes = std::get<2>(data); auto & contour_nodes_subfacets = std::get<3>(data); // compute effective stress in all cohesive material linear sequent. for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; mat_coh->computeEffectiveStresses(); } // identify facets on contour where cohesives should be inserted std::map> mat_index_facet_nbs_crack_nbs = findCriticalFacetsOnContourByNodes( contour_subfacets_coh_el, surface_subfacets_crack_nb, contour_nodes_subfacets, surface_nodes, av_stress_threshold); // insert cohesives depending on a material for (auto & pair : mat_index_facet_nbs_crack_nbs) { auto mat_index = pair.first; auto & facet_nbs_crack_nbs = pair.second; auto & mat = model.getMaterial(mat_index); auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; mat_coh->insertCohesiveElements(facet_nbs_crack_nbs, type_facet, false); } // insert the most stressed cohesive element nb_new_elements = inserter.insertElements(); // communicate crack numbers communicateCrackNumbers(); AKANTU_DEBUG_OUT(); return nb_new_elements; } template UInt ASRTools::insertCohesiveElementsOnContourByNodes<2>(Real av_stress_threshold); template UInt ASRTools::insertCohesiveElementsOnContourByNodes<3>(Real av_stress_threshold); /* ------------------------------------------------------------------- */ template UInt ASRTools::insertLoopsOfStressedCohesives(Real min_dot) { AKANTU_DEBUG_IN(); auto & coh_model = dynamic_cast(model); CohesiveElementInserter & inserter = coh_model.getElementInserter(); if (not coh_model.getIsExtrinsic()) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } coh_model.interpolateStress(); // compute effective stress in all cohesive material linear sequent. for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; mat_coh->computeEffectiveStresses(); } const Mesh & mesh_facets = coh_model.getMeshFacets(); UInt nb_new_elements(0); auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); // find global crack topology auto data = determineCrackSurface(); auto & contour_subfacets_coh_el = std::get<0>(data); auto & surface_subfacets_crack_nb = std::get<1>(data); auto & surface_nodes = std::get<2>(data); auto & contour_nodes_subfacets = std::get<3>(data); // identify facets on contour where cohesives should be inserted std::map> mat_index_facet_nbs_crack_nbs = findStressedFacetsLoops(contour_subfacets_coh_el, surface_subfacets_crack_nb, contour_nodes_subfacets, surface_nodes, min_dot); // insert cohesives depending on a material for (auto & pair : mat_index_facet_nbs_crack_nbs) { auto mat_index = pair.first; auto & facet_nbs_crack_nbs = pair.second; auto & mat = model.getMaterial(mat_index); auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; mat_coh->insertCohesiveElements(facet_nbs_crack_nbs, type_facet, false); } // insert the most stressed cohesive element nb_new_elements = inserter.insertElements(); // communicate crack numbers communicateCrackNumbers(); AKANTU_DEBUG_OUT(); return nb_new_elements; } template UInt ASRTools::insertLoopsOfStressedCohesives<2>(Real min_dot); template UInt ASRTools::insertLoopsOfStressedCohesives<3>(Real min_dot); /* --------------------------------------------------------------- */ std::map> ASRTools::findCriticalFacetsOnContourByNodes( const std::map & contour_subfacets_coh_el, const std::map & /*surface_subfacets_crack_nb*/, const std::map> & contour_nodes_subfacets, const std::set & /*surface_nodes*/, Real av_stress_threshold) { AKANTU_DEBUG_IN(); const Mesh & mesh = model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); std::map> mat_index_facet_nbs_crack_nbs; std::map contour_to_single_facet_map; // TODO : make iteration on multiple facet types auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); auto nb_facets = mesh.getNbElement(type_facet); // skip if no facets of this type are present if (not nb_facets) { return mat_index_facet_nbs_crack_nbs; } // create a copy to remove nodes where single facets were inserted auto contour_nodes_subfacets_copy = contour_nodes_subfacets; // first loop for single elements for (auto && pair : contour_nodes_subfacets) { auto cent_node = pair.first; auto subfacets = pair.second; if (subfacets.size() != 2) { std::cout << "Number of contour segments connected to the node " << cent_node << " is not equal to 2. Skipping" << std::endl; continue; } auto it = subfacets.begin(); auto starting_segment(*it); std::advance(it, 1); auto ending_segment(*it); auto starting_coh_el = contour_subfacets_coh_el.at(starting_segment); auto ending_coh_el = contour_subfacets_coh_el.at(ending_segment); auto & crack_numbers = mesh.getData( "crack_numbers", starting_coh_el.type, starting_coh_el.ghost_type); auto crack_nb = crack_numbers(starting_coh_el.element); auto starting_facet = mesh_facets.getSubelementToElement(starting_coh_el)(0); auto ending_facet = mesh_facets.getSubelementToElement(ending_coh_el)(0); // if an outstanding triangle element - skip if (starting_facet == ending_facet) { contour_nodes_subfacets_copy.erase(cent_node); continue; } Array limit_facets, limit_segments; limit_facets.push_back(starting_facet); limit_facets.push_back(ending_facet); limit_segments.push_back(starting_segment); limit_segments.push_back(ending_segment); auto single_facet_loop = findSingleFacetLoop(limit_facets, limit_segments, cent_node); if (not single_facet_loop.size()) continue; // evaluate average effective stress over facets of the loop auto av_eff_stress = averageEffectiveStressInMultipleFacets(single_facet_loop); // if average effective stress exceeds threshold->decompose per mat if (av_eff_stress >= av_stress_threshold) { decomposeFacetLoopPerMaterial(single_facet_loop, crack_nb, mat_index_facet_nbs_crack_nbs); // update the map for (auto && limit_segment : limit_segments) { contour_to_single_facet_map[limit_segment] = single_facet_loop(0); } // remove node from the looping map contour_nodes_subfacets_copy.erase(cent_node); } } // second loop for long facet loops on a reduced map for (auto && pair : contour_nodes_subfacets_copy) { auto cent_node = pair.first; auto subfacets = pair.second; if (subfacets.size() != 2) { std::cout << "Number of contour segments connected to the node " << cent_node << " is not equal to 2. Skipping" << std::endl; continue; } auto it = subfacets.begin(); auto starting_segment(*it); std::advance(it, 1); auto ending_segment(*it); auto starting_coh_el = contour_subfacets_coh_el.at(starting_segment); auto ending_coh_el = contour_subfacets_coh_el.at(ending_segment); auto & crack_numbers = mesh.getData( "crack_numbers", starting_coh_el.type, starting_coh_el.ghost_type); auto crack_nb = crack_numbers(starting_coh_el.element); auto starting_facet = mesh_facets.getSubelementToElement(starting_coh_el)(0); auto ending_facet = mesh_facets.getSubelementToElement(ending_coh_el)(0); Array limit_facets, limit_segments; limit_facets.push_back(starting_facet); limit_facets.push_back(ending_facet); limit_segments.push_back(starting_segment); limit_segments.push_back(ending_segment); Array preinserted_facets(2, 1, ElementNull); for (UInt i = 0; i < 2; i++) { auto it = contour_to_single_facet_map.find(limit_segments(i)); if (it != contour_to_single_facet_map.end()) { preinserted_facets(i) = it->second; } } auto facet_loop = findFacetsLoopByGraphByArea( limit_facets, limit_segments, preinserted_facets, cent_node); if (not facet_loop.size()) continue; // evaluate average effective stress over facets of the loop auto av_eff_stress = averageEffectiveStressInMultipleFacets(facet_loop); // if average effective stress exceeds threshold->decompose per mat if (av_eff_stress >= av_stress_threshold) { decomposeFacetLoopPerMaterial(facet_loop, crack_nb, mat_index_facet_nbs_crack_nbs); } } return mat_index_facet_nbs_crack_nbs; } /* --------------------------------------------------------------- */ std::map> ASRTools::findStressedFacetsLoops( const std::map & contour_subfacets_coh_el, const std::map & /*surface_subfacets_crack_nb*/, const std::map> & contour_nodes_subfacets, const std::set & /*surface_nodes*/, Real min_dot) { AKANTU_DEBUG_IN(); const Mesh & mesh = model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); auto dim = mesh.getSpatialDimension(); std::map> mat_index_facet_nbs_crack_nbs; std::map>> node_to_crack_nb_facet_loop; // auto && comm = akantu::Communicator::getWorldCommunicator(); // auto prank = comm.whoAmI(); // clear the nodes_eff_stress array this->nodes_eff_stress.zero(); // TODO : make iteration on multiple facet types auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); auto nb_facets = mesh.getNbElement(type_facet); // skip if no facets of this type are present if (nb_facets) { for (auto && pair : contour_nodes_subfacets) { auto cent_node = pair.first; auto subfacets = pair.second; if (subfacets.size() != 2) { // std::cout << "Number of contour segments connected to the node " // << cent_node << " is not equal to 2. Skipping" << // std::endl; continue; } auto it = subfacets.begin(); auto starting_segment(*it); std::advance(it, 1); auto ending_segment(*it); auto starting_coh_el = contour_subfacets_coh_el.at(starting_segment); auto ending_coh_el = contour_subfacets_coh_el.at(ending_segment); auto & crack_numbers = mesh.getData( "crack_numbers", starting_coh_el.type, starting_coh_el.ghost_type); auto crack_nb = crack_numbers(starting_coh_el.element); auto starting_facet = mesh_facets.getSubelementToElement(starting_coh_el)(0); auto ending_facet = mesh_facets.getSubelementToElement(ending_coh_el)(0); Array limit_facets, limit_segments; limit_facets.push_back(starting_facet); limit_facets.push_back(ending_facet); limit_segments.push_back(starting_segment); limit_segments.push_back(ending_segment); auto facet_loop = findStressedFacetLoopAroundNode( limit_facets, limit_segments, cent_node, min_dot); // store the facet loop in the map node_to_crack_nb_facet_loop[cent_node] = std::make_pair(crack_nb, facet_loop); // evaluate average effective stress over facets of the loop auto av_eff_stress = averageEffectiveStressInMultipleFacets(facet_loop); // assign stress value to the node this->nodes_eff_stress(cent_node) = av_eff_stress; } } // communicate effective stresses by erasing the lowest values communicateEffStressesOnNodes(); for (auto && pair : node_to_crack_nb_facet_loop) { auto && cent_node = pair.first; auto && crack_nb = pair.second.first; auto && facet_loop = pair.second.second; if (this->nodes_eff_stress(cent_node) > 0) { decomposeFacetLoopPerMaterial(facet_loop, crack_nb, mat_index_facet_nbs_crack_nbs); } } return mat_index_facet_nbs_crack_nbs; } /* --------------------------------------------------------------- */ Real ASRTools::averageEffectiveStressInMultipleFacets( Array & facet_loop) { if (not facet_loop.size()) return 0; auto & coh_model = dynamic_cast(model); const FEEngine & fe_engine = model.getFEEngine("FacetsFEEngine"); Array eff_stress_array; ElementType type_facet; GhostType gt; Array facet_nbs; for (auto & facet : facet_loop) { type_facet = facet.type; gt = facet.ghost_type; facet_nbs.push_back(facet.element); auto & facet_material_index = coh_model.getFacetMaterial(type_facet, gt)(facet.element); auto & mat = model.getMaterial(facet_material_index); auto * mat_coh = dynamic_cast(&mat); if (mat_coh == nullptr) return 0; auto && mat_facet_filter = mat_coh->getFacetFilter(type_facet, gt); UInt facet_local_num = mat_facet_filter.find(facet.element); AKANTU_DEBUG_ASSERT(facet_local_num != UInt(-1), "Local facet number was not identified"); auto & eff_stress = mat.getInternal("effective_stresses")( type_facet, gt)(facet_local_num); for (__attribute__((unused)) UInt quad : arange(fe_engine.getNbIntegrationPoints(type_facet))) eff_stress_array.push_back(eff_stress); } Real int_eff_stress = fe_engine.integrate(eff_stress_array, type_facet, gt, facet_nbs); Array dummy_array( facet_loop.size() * fe_engine.getNbIntegrationPoints(type_facet), 1, 1.); Real loop_area = fe_engine.integrate(dummy_array, type_facet, gt, facet_nbs); return int_eff_stress / loop_area; } /* --------------------------------------------------------------- */ void ASRTools::decomposeFacetLoopPerMaterial( const Array & facet_loop, UInt crack_nb, std::map> & mat_index_facet_nbs_crack_nbs) { auto & coh_model = dynamic_cast(model); for (auto & facet : facet_loop) { ElementType type_facet = facet.type; GhostType gt = facet.ghost_type; auto & facet_material_index = coh_model.getFacetMaterial(type_facet, gt)(facet.element); mat_index_facet_nbs_crack_nbs[facet_material_index][facet.element] = crack_nb; } } /* --------------------------------------------------------------- */ void ASRTools::applyEigenOpening(Real eigen_strain) { auto & coh_model = dynamic_cast(model); if (not coh_model.getIsExtrinsic()) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } const Mesh & mesh_facets = coh_model.getMeshFacets(); const UInt dim = model.getSpatialDimension(); for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast(&mat); if (mat_coh == nullptr) continue; for (auto gt : ghost_types) { 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); if (not model.getMesh().getNbElement(type_cohesive, gt)) continue; const Array & elem_filter = mat_coh->getElementFilter(type_cohesive, gt); if (not elem_filter.size()) { continue; } // access openings of cohesive elements auto & eig_opening = mat_coh->getEigenOpening(type_cohesive, gt); auto & normals = mat_coh->getNormals(type_cohesive, gt); auto eig_op_it = eig_opening.begin(dim); auto normals_it = normals.begin(dim); // apply eigen strain as opening for (auto el_order : arange(elem_filter.size())) { Element cohesive{type_cohesive, elem_filter(el_order), gt}; const Vector & facets_to_cohesive = mesh_facets.getSubelementToElement(cohesive); Real facet_indiam = MeshUtils::getInscribedCircleDiameter( model, facets_to_cohesive(0)); auto eigen_opening = eigen_strain * facet_indiam; Vector normal(normals_it[el_order * nb_quad_cohesive]); for (UInt i : arange(nb_quad_cohesive)) { auto eig_op = eig_op_it[el_order * nb_quad_cohesive + i]; eig_op = normal * eigen_opening; } } } } } } /* --------------------------------------------------------------- */ void ASRTools::applyEigenOpeningGelVolumeBased(Real gel_volume_ratio) { auto & coh_model = dynamic_cast(model); const Mesh & mesh_facets = coh_model.getMeshFacets(); const UInt dim = model.getSpatialDimension(); // compute total asr gel volume if (!this->volume) computeModelVolume(); Real gel_volume = this->volume * gel_volume_ratio; // compute total area of existing cracks auto data_agg = computeCrackData("agg-agg"); auto data_agg_mor = computeCrackData("agg-mor"); auto data_mor = computeCrackData("mor-mor"); auto area_agg = std::get<0>(data_agg); auto area_agg_mor = std::get<0>(data_agg_mor); auto area_mor = std::get<0>(data_mor); Real total_area = area_agg + area_agg_mor + area_mor; // compute necessary opening to cover fit gel volume by crack area Real average_opening = gel_volume / total_area; for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast(&mat); if (mat_coh == nullptr) continue; for (auto gt : ghost_types) { 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); if (not model.getMesh().getNbElement(type_cohesive, gt)) continue; const Array & elem_filter = mat_coh->getElementFilter(type_cohesive, gt); if (not elem_filter.size()) { continue; } // access openings of cohesive elements auto & eig_opening = mat_coh->getEigenOpening(type_cohesive, gt); auto & normals = mat_coh->getNormals(type_cohesive, gt); auto eig_op_it = eig_opening.begin(dim); auto normals_it = normals.begin(dim); // apply average opening directly as eigen for (auto el_order : arange(elem_filter.size())) { Vector normal(normals_it[el_order * nb_quad_cohesive]); for (UInt i : arange(nb_quad_cohesive)) { auto eig_op = eig_op_it[el_order * nb_quad_cohesive + i]; eig_op = normal * average_opening; } } } } } } /* --------------------------------------------------------------- */ void ASRTools::applyEigenOpeningToInitialCrack( Real gel_volume_ratio, const Array> & asr_facets) { auto & coh_model = dynamic_cast(model); const FEEngine & fe_engine = model.getFEEngine("CohesiveFEEngine"); const UInt dim = model.getSpatialDimension(); const Mesh & mesh_facets = coh_model.getMeshFacets(); auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); UInt nb_quad_cohesive = model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); const Array & material_index_vec = model.getMaterialByElement(type_cohesive); const Array & material_local_numbering_vec = model.getMaterialLocalNumbering(type_cohesive); // compute total asr gel volume if (!this->volume) computeModelVolume(); Real gel_volume = this->volume * gel_volume_ratio; Real gel_volume_per_crack = gel_volume / asr_facets.size(); // iterate through each crack for (auto & single_site_facets : asr_facets) { // compute area of current ASR site // form cohesive element filter std::set site_cohesives; for (auto & asr_facet : single_site_facets) { // find connected cohesive auto & connected_els = mesh_facets.getElementToSubelement(asr_facet); for (auto & connected_el : connected_els) { if (connected_el.type == type_cohesive) { site_cohesives.emplace(connected_el); break; } } } // integrate 1 over identified element filter Real site_area{0}; for (auto & coh_el : site_cohesives) { Array single_el_array; single_el_array.push_back(coh_el.element); Array area(fe_engine.getNbIntegrationPoints(type_cohesive), 1, 1.); site_area += fe_engine.integrate(area, coh_el.type, coh_el.ghost_type, single_el_array); } // compute necessary opening to cover fit gel volume by crack area Real average_opening = gel_volume_per_crack / site_area; // iterate through each cohesive and apply eigen opening for (auto && coh_el : site_cohesives) { Material & mat = model.getMaterial(material_index_vec(coh_el.element)); auto * mat_coh = dynamic_cast(&mat); if (mat_coh == nullptr) continue; UInt coh_el_local_num = material_local_numbering_vec(coh_el.element); // access openings of cohesive elements auto & eig_opening = mat_coh->getEigenOpening(coh_el.type, coh_el.ghost_type); auto & normals = mat_coh->getNormals(coh_el.type, coh_el.ghost_type); auto eig_op_it = eig_opening.begin(dim); auto normals_it = normals.begin(dim); Vector normal(normals_it[coh_el_local_num * nb_quad_cohesive]); for (UInt i : arange(nb_quad_cohesive)) { Vector eig_op(eig_op_it[coh_el_local_num * nb_quad_cohesive + i]); eig_op = normal * average_opening; } } } } /* --------------------------------------------------------------- */ +void ASRTools::applyPointForceToAsrCentralNodes(Real force_norm) { + auto && mesh = model.getMesh(); + auto dim = mesh.getSpatialDimension(); + auto & forces = model.getExternalForce(); + auto & pos = mesh.getNodes(); + auto it_force = make_view(forces, dim).begin(); + auto it_pos = make_view(pos, dim).begin(); + + for (auto && data : zip(asr_central_node_pairs, asr_normals_pairs)) { + auto && node_pair = std::get<0>(data); + auto && normals_pair = std::get<1>(data); + auto node1 = node_pair.first; + auto node2 = node_pair.second; + auto normal1 = normals_pair.first; + auto normal2 = normals_pair.second; + Vector force1 = it_force[node1]; + Vector force2 = it_force[node2]; + force1 = normal1 * force_norm / 2.; + force2 = normal2 * force_norm / 2.; + } +} +/* --------------------------------------------------------------- */ void ASRTools::outputCrackData(std::ofstream & file_output, Real time) { auto data_agg = computeCrackData("agg-agg"); auto data_agg_mor = computeCrackData("agg-mor"); auto data_mor = computeCrackData("mor-mor"); auto area_agg = std::get<0>(data_agg); auto vol_agg = std::get<1>(data_agg); auto asr_vol_agg = std::get<2>(data_agg); auto area_agg_mor = std::get<0>(data_agg_mor); auto vol_agg_mor = std::get<1>(data_agg_mor); auto asr_vol_agg_mor = std::get<2>(data_agg_mor); auto area_mor = std::get<0>(data_mor); auto vol_mor = std::get<1>(data_mor); auto asr_vol_mor = std::get<2>(data_mor); Real total_area = area_agg + area_agg_mor + area_mor; Real total_volume = vol_agg + vol_agg_mor + vol_mor; Real asr_volume = asr_vol_agg + asr_vol_agg_mor + asr_vol_mor; auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); - if (prank == 0) + if (prank == 0) { file_output << time << "," << asr_volume << "," << area_agg << "," << area_agg_mor << "," << area_mor << "," << vol_agg << "," << vol_agg_mor << "," << vol_mor << "," << total_area << "," << total_volume << std::endl; + } } /* --------------------------------------------------------------- */ std::tuple ASRTools::computeCrackData(const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); auto * mat_coh = dynamic_cast(&mat); AKANTU_DEBUG_ASSERT(mat_coh != nullptr, "Crack data can not be computed on material with id" + material_name); const ElementTypeMapArray & filter_map = mat_coh->getElementFilter(); const FEEngine & fe_engine = model.getFEEngine("CohesiveFEEngine"); Real crack_volume{0}; Real crack_area{0}; Real ASR_volume{0}; // Loop over the cohesive element types for (auto & element_type : filter_map.elementTypes(dim, gt, _ek_cohesive)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; /// compute normal norm of real opening = opening + eigen opening auto opening_it = mat_coh->getOpening(element_type, gt).begin(dim); auto eigen_opening = mat_coh->getEigenOpening(element_type, gt); auto eigen_opening_it = eigen_opening.begin(dim); Array normal_eigen_opening_norm( filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1); auto normal_eigen_opening_norm_it = normal_eigen_opening_norm.begin(); Array real_normal_opening_norm( filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1); auto real_normal_opening_norm_it = real_normal_opening_norm.begin(); auto normal_it = mat_coh->getNormals(element_type, gt).begin(dim); auto normal_end = mat_coh->getNormals(element_type, gt).end(dim); /// loop on each quadrature point for (; normal_it != normal_end; ++opening_it, ++eigen_opening_it, ++normal_it, ++real_normal_opening_norm_it, ++normal_eigen_opening_norm_it) { auto real_opening = *opening_it + *eigen_opening_it; *real_normal_opening_norm_it = real_opening.dot(*normal_it); Vector eig_opening(*eigen_opening_it); *normal_eigen_opening_norm_it = eig_opening.dot(*normal_it); } ASR_volume += fe_engine.integrate(normal_eigen_opening_norm, element_type, gt, filter); crack_volume += fe_engine.integrate(real_normal_opening_norm, element_type, gt, filter); Array area( filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); crack_area += fe_engine.integrate(area, element_type, gt, filter); } /// do not communicate if model is multi-scale auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(ASR_volume, SynchronizerOperation::_sum); comm.allReduce(crack_volume, SynchronizerOperation::_sum); comm.allReduce(crack_area, SynchronizerOperation::_sum); return std::make_tuple(crack_area, crack_volume, ASR_volume); } +/* --------------------------------------------------------------- */ +void ASRTools::outputCrackOpenings(std::ofstream & file_output, Real time) { + const Communicator & comm = Communicator::getWorldCommunicator(); + UInt prank = comm.whoAmI(); + UInt nb_proc = comm.getNbProc(); + Real global_nb_openings(asr_central_node_pairs.size()); + comm.allReduce(global_nb_openings, SynchronizerOperation::_sum); + /// shift crack numbers by the number of cracks on lower processors + UInt nb_cracks{asr_central_node_pairs.size()}; + comm.exclusiveScan(nb_cracks); + UInt max_nb_cracks{asr_central_node_pairs.size()}; + comm.allReduce(max_nb_cracks, SynchronizerOperation::_max); + Array openings(nb_proc, max_nb_cracks, -1.); + + auto && mesh = model.getMesh(); + auto dim = mesh.getSpatialDimension(); + auto & disp = model.getDisplacement(); + auto & pos = mesh.getNodes(); + auto it_disp = make_view(disp, dim).begin(); + auto it_pos = make_view(pos, dim).begin(); + auto it_open = make_view(openings, max_nb_cracks).begin(); + auto open_end = make_view(openings, max_nb_cracks).end(); + + for (UInt i = 0; i < asr_central_node_pairs.size(); i++) { + auto && node_pair = asr_central_node_pairs(i); + auto && normals_pair = asr_normals_pairs(i); + + auto node1 = node_pair.first; + auto node2 = node_pair.second; + auto normal1 = normals_pair.first; + Vector disp1 = it_disp[node1]; + Vector disp2 = it_disp[node2]; + Vector rel_disp = disp1 - disp2; + Real norm_opening = rel_disp.dot(normal1); + openings(prank, i) = norm_opening; + } + comm.allGather(openings); + if (prank == 0) { + file_output << time; + for (; it_open != open_end; ++it_open) { + Vector openings_per_proc = *it_open; + for (auto && opening : openings_per_proc) { + if (opening != -1.) { + file_output << "," << opening; + } + } + } + file_output << std::endl; + } +} + +/* --------------------------------------------------------------- */ +void ASRTools::identifyASRCentralNodePairsAndNormals() { + auto && mesh = model.getMesh(); + auto && mesh_facets = mesh.getMeshFacets(); + auto dim = mesh.getSpatialDimension(); + auto & forces = model.getExternalForce(); + auto & pos = mesh.getNodes(); + auto type_facet = *mesh_facets.elementTypes(dim - 1).begin(); + auto && fe_engine = model.getFEEngine("FacetsFEEngine"); + UInt nb_quad_points = fe_engine.getNbIntegrationPoints(type_facet); + + // get normal to the initial positions + const auto & init_pos = mesh.getNodes(); + Array quad_normals(0, dim); + fe_engine.computeNormalsOnIntegrationPoints(init_pos, quad_normals, + type_facet); + auto normals_it = quad_normals.begin(dim); + auto it_force = make_view(forces, dim).begin(); + auto it_pos = make_view(pos, dim).begin(); + + this->asr_central_node_pairs.resize(asr_facets.size()); + this->asr_normals_pairs.resize(asr_facets.size()); + for (auto && data : + zip(asr_facets, asr_central_node_pairs, asr_normals_pairs)) { + auto && facets_per_site = std::get<0>(data); + auto & node_pair = std::get<1>(data); + node_pair = std::make_pair(UInt(-1), UInt(-1)); + auto & asr_normals_pair = std::get<2>(data); + + // find first node + auto && facet_conn = mesh_facets.getConnectivity(facets_per_site(0)); + for (auto && node : facet_conn) { + auto ret = asr_central_nodes.find(node); + if (ret != UInt(-1)) { + node_pair.first = node; + break; + } + } + AKANTU_DEBUG_ASSERT(node_pair.first != UInt(-1), + "First node was not identified"); + // compute first normal + Vector normal(dim, 0); + for (auto && facet : facets_per_site) { + // skip out opposite side facets + auto && facet_nodes = mesh_facets.getConnectivity(facet); + auto ret = + std::find(facet_nodes.begin(), facet_nodes.end(), node_pair.first); + if (ret == facet_nodes.end()) { + continue; + } + + // find connected solid element + auto && facet_neighbors = mesh_facets.getElementToSubelement(facet); + Element solid_el{ElementNull}; + for (auto && neighbor : facet_neighbors) { + if (mesh.getKind(neighbor.type) == _ek_regular) { + solid_el = neighbor; + break; + } + } + AKANTU_DEBUG_ASSERT(solid_el != ElementNull, + "Couldn't identify neighboring solid el"); + + // find the out-of-plane solid node + auto && solid_nodes = mesh.getConnectivity(solid_el); + auto it = + std::find(solid_nodes.begin(), solid_nodes.end(), node_pair.first); + AKANTU_DEBUG_ASSERT(it != solid_nodes.end(), + "Solid el does not contain central node"); + + UInt out_node(-1); + for (auto && solid_node : solid_nodes) { + auto ret = + std::find(facet_nodes.begin(), facet_nodes.end(), solid_node); + if (ret == facet_nodes.end()) { + out_node = solid_node; + break; + } + } + AKANTU_DEBUG_ASSERT(out_node != UInt(-1), + "Couldn't identify out-of-plane node"); + + // build the reference vector + Vector ref_vector(dim); + mesh_facets.getBarycenter(facet, ref_vector); + Vector pos(mesh.getNodes().begin(dim)[out_node]); + ref_vector = pos - ref_vector; + + // check if ref vector and the normal are in the same half space + Vector facet_normal = normals_it[facet.element * nb_quad_points]; + if (ref_vector.dot(facet_normal) < 0) { + facet_normal *= -1; + } + normal += facet_normal; + } + asr_normals_pair.first = normal / normal.norm(); + asr_normals_pair.second = normal / normal.norm() * (-1.); + + // find the second one + Vector first_node_pos = it_pos[node_pair.first]; + for (auto && asr_central_node : asr_central_nodes) { + Vector node_pos = it_pos[asr_central_node]; + if ((first_node_pos == node_pos) && + (asr_central_node != node_pair.first)) { + node_pair.second = asr_central_node; + break; + } + } + AKANTU_DEBUG_ASSERT(node_pair.second != UInt(-1), + "Second node was not identified"); + } +} + /* ----------------------------------------------------------------- */ std::tuple, std::map, std::set, std::map>> ASRTools::determineCrackSurface() { AKANTU_DEBUG_IN(); Mesh & mesh = model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); const UInt dim = mesh.getSpatialDimension(); std::map contour_subfacets_coh_el; std::map surface_subfacets_crack_nb; std::set surface_nodes; std::map> contour_nodes_subfacets; std::set visited_subfacets; // // first loop on facets (doesn't include duplicated ones) // for (auto & type_facet : // mesh.elementTypes(dim - 1, _not_ghost, _ek_regular)) { // auto nb_facets = mesh.getNbElement(type_facet); // std::cout << "Nb facets = " << nb_facets << std::endl; // if (not nb_facets) { // continue; // } // for (auto facet_nb : arange(nb_facets)) { // Element facet{type_facet, facet_nb, _not_ghost}; // searchCrackSurfaceInSingleFacet( // facet, contour_subfacets_coh_el, surface_subfacets_crack_nb, // surface_nodes, contour_nodes_subfacets, visited_subfacets); // } // } // second loop on cohesives (includes duplicated facets) for (auto gt : ghost_types) { for (auto & type_cohesive : mesh.elementTypes(dim, gt, _ek_cohesive)) { auto nb_cohesives = mesh.getNbElement(type_cohesive, gt); if (not nb_cohesives) { continue; } for (auto cohesive_nb : arange(nb_cohesives)) { Element cohesive{type_cohesive, cohesive_nb, gt}; auto && facets_to_cohesive = mesh_facets.getSubelementToElement(cohesive); for (auto coh_facet : facets_to_cohesive) { searchCrackSurfaceInSingleFacet( coh_facet, contour_subfacets_coh_el, surface_subfacets_crack_nb, surface_nodes, contour_nodes_subfacets, visited_subfacets); } } } } // remove external nodes from the surface nodes list for (auto & pair : contour_nodes_subfacets) { surface_nodes.erase(pair.first); } return std::make_tuple(contour_subfacets_coh_el, surface_subfacets_crack_nb, surface_nodes, contour_nodes_subfacets); } /* ----------------------------------------------------------------- */ void ASRTools::searchCrackSurfaceInSingleFacet( Element & facet, std::map & contour_subfacets_coh_el, std::map & surface_subfacets_crack_nb, std::set & surface_nodes, std::map> & contour_nodes_subfacets, std::set & visited_subfacets) { Mesh & mesh = model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); // iterate through subfacets of a facet searching for cohesives auto && subfacets_to_facet = mesh_facets.getSubelementToElement(facet); for (auto & subfacet : subfacets_to_facet) { // discard subfacets on the border of ghost area if (subfacet == ElementNull) continue; // discard pure ghost subfacets bool pure_ghost_subf{false}; auto && subfacet_conn = mesh_facets.getConnectivity(subfacet); for (auto & subfacet_node : subfacet_conn) { if (mesh.isPureGhostNode(subfacet_node)) pure_ghost_subf = true; } if (pure_ghost_subf) continue; auto ret = visited_subfacets.emplace(subfacet); if (not ret.second) continue; std::map connected_cohesives; auto && facets_to_subfacet = mesh_facets.getElementToSubelement(subfacet); for (auto & neighbor_facet : facets_to_subfacet) { auto & connected_to_facet = mesh_facets.getElementToSubelement(neighbor_facet); for (auto & connected_el : connected_to_facet) { if (connected_el.type == _not_defined) goto nextfacet; if (mesh.getKind(connected_el.type) == _ek_cohesive) { connected_cohesives[connected_el]++; } } nextfacet:; } if (connected_cohesives.size() == 1 and connected_cohesives.begin()->second == 2) { contour_subfacets_coh_el[subfacet] = connected_cohesives.begin()->first; Vector contour_subfacet_nodes = mesh_facets.getConnectivity(subfacet); for (auto & contour_subfacet_node : contour_subfacet_nodes) { contour_nodes_subfacets[contour_subfacet_node].emplace(subfacet); } } else if (connected_cohesives.size() > 1) { auto coh_element = connected_cohesives.begin()->first; UInt crack_nb = mesh.getData("crack_numbers", coh_element.type, coh_element.ghost_type)(coh_element.element); surface_subfacets_crack_nb[subfacet] = crack_nb; // add subfacet nodes to the set for (auto & subfacet_node : mesh_facets.getConnectivity(subfacet)) surface_nodes.emplace(subfacet_node); } } } /* ------------------------------------------------------------------- */ template UInt ASRTools::updateDeltaMax() { AKANTU_DEBUG_IN(); auto & coh_model = dynamic_cast(model); if (not coh_model.getIsExtrinsic()) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } UInt nb_updated_quads(0); for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; for (auto type : mat_coh->getElementFilter().elementTypes(dim, _not_ghost, _ek_cohesive)) { nb_updated_quads += mat_coh->updateDeltaMax(type, _not_ghost); } } // summ all updated quads auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(nb_updated_quads, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_updated_quads; } template UInt ASRTools::updateDeltaMax<2>(); template UInt ASRTools::updateDeltaMax<3>(); /* ------------------------------------------------------------------- */ template std::tuple ASRTools::getMaxDeltaMaxExcess() { AKANTU_DEBUG_IN(); auto & coh_model = dynamic_cast(model); coh_model.interpolateStress(); coh_model.assembleInternalForces(); Real proc_max_delta_max_excess{0}; Element critical_coh_el{ElementNull}; UInt global_nb_penetr_quads{0}; for (auto && mat : model.getMaterials()) { auto * mat_coh = dynamic_cast *>(&mat); if (mat_coh == nullptr) continue; for (auto type : mat_coh->getElementFilter().elementTypes(dim, _not_ghost, _ek_cohesive)) { auto data = mat_coh->computeMaxDeltaMaxExcess(type, _not_ghost); auto mat_max_delta_max_excess = std::get<0>(data); global_nb_penetr_quads += std::get<2>(data); if (mat_max_delta_max_excess > proc_max_delta_max_excess) { proc_max_delta_max_excess = mat_max_delta_max_excess; critical_coh_el = std::get<1>(data); } } } // communicate the highest excess over different processors auto global_max_delta_max = proc_max_delta_max_excess; auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(global_max_delta_max, SynchronizerOperation::_max); comm.allReduce(global_nb_penetr_quads, SynchronizerOperation::_sum); // find the critical coh between processors UInt coh_nb = critical_coh_el.element; if (proc_max_delta_max_excess != global_max_delta_max) coh_nb = 0; comm.allReduce(coh_nb, SynchronizerOperation::_max); critical_coh_el.element = coh_nb; AKANTU_DEBUG_OUT(); return std::make_tuple(global_max_delta_max, critical_coh_el, global_nb_penetr_quads); } template std::tuple ASRTools::getMaxDeltaMaxExcess<2>(); template std::tuple ASRTools::getMaxDeltaMaxExcess<3>(); } // namespace akantu diff --git a/extra_packages/asr-tools/src/asr_tools.hh b/extra_packages/asr-tools/src/asr_tools.hh index 0c9dcdaed..cee6f3ff8 100644 --- a/extra_packages/asr-tools/src/asr_tools.hh +++ b/extra_packages/asr-tools/src/asr_tools.hh @@ -1,1253 +1,1329 @@ /** * @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 multiple blocks of cohesive elements void insertASRCohesivesRandomly(const UInt & nb_coh_elem, std::string matrix_mat_name, Real gap_ratio); /// insert multiple blocks of cohesive elements void insertASRCohesivesRandomly3D(const UInt & nb_coh_elem, std::string matrix_mat_name, Real gap_ratio); /// ASR insertion for 1st order 3D elements void insertASRCohesiveLoops3D(const UInt & nb_insertions, std::string facet_mat_name, Real gap_ratio); /// no facets are added into the mesh template void insertPreCracks(const UInt & nb_insertions, std::string facet_mat_name); /// insert block of cohesive elements based on the coord of the central void insertASRCohesivesByCoords(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(); /// on master or slave nodes with potential facet loops around /// put 0 eff stress for the smallest value void communicateEffStressesOnNodes(); /// pick facets by passed coordinates void pickFacetsByCoord(const Matrix & positions); /// pick facets randomly within specified material void pickFacetsRandomly(UInt nb_insertions, std::string facet_mat_name); /// build closed facets loop around random point of specified material and /// returns number of successfully inserted ASR sites UInt closedFacetsLoopAroundPoint(UInt nb_insertions, std::string mat_name); /// finds facet loops around a point and insert cohesives template UInt insertCohesiveLoops(UInt nb_insertions, std::string mat_name); /// builds a facet loop around a point from starting to ending segm-s Array findFacetsLoopFromSegment2Segment( Element starting_facet, Element starting_segment, Element ending_segment, UInt cent_node, Real max_dot, UInt material_id, bool check_asr_facets); /// builds a facet loop around a point from starting to ending segment using /// the Dijkstra shortest path algorithm in boost library /// uses sum of distances to 2 incenters as weight Array findFacetsLoopByGraphByDist(const Array & limit_facets, const Array & limit_segments, const UInt & cent_node); /// use area as a weight Array findFacetsLoopByGraphByArea(const Array & limit_facets, const Array & limit_segments, const Array & preinserted_facets, const UInt & cent_node); /// include only facets that have eff_stress >= 1 Array findStressedFacetLoopAroundNode(const Array & limit_facets, const Array & limit_segments, const UInt & cent_node, Real min_dot); /// insert single facets before searching for the long loops Array findSingleFacetLoop(const Array & limit_facets, const Array & limit_segments, const UInt & cent_node); /// check if the cohesive can be inserted and nodes are not on the partition /// border and ASR nodes bool isFacetAndNodesGood(const Element & facet, UInt material_id = UInt(-1), bool check_asr_facets = false); /// check if 2 facets are bounding a common solid element bool belong2SameElement(const Element & facet1, const Element & facet2); /// check if the node is surrounded by the material bool isNodeWithinMaterial(Element & node, UInt material_id); /// pick two neighbors of a central facet in 2D: returns true if success bool pickFacetNeighborsOld(Element & cent_facet); /// version working for both 2d and 3d bool pickFacetNeighbors(Element & cent_facet); /// optimise doubled facets group, insert facets, and cohesive elements, /// update connectivities void insertOppositeFacetsAndCohesives(); /// 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); /// same in 3D void insertGap3D(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 insertCohesiveElementsSelectively(); /// insert multiple cohesives on contour template UInt insertCohesiveElementsOnContour(); /// insert multiple cohesives on contour by looping on contour nodes template UInt insertCohesiveElementsOnContourByNodes(Real av_stress_threshold); /// insert cohesives on closed loops of stressed facets template UInt insertLoopsOfStressedCohesives(Real min_dot); /// find critical facets by looping on contour nodes std::map> findCriticalFacetsOnContourByNodes( const std::map & contour_subfacets_coh_el, const std::map & surface_subfacets_crack_nb, const std::map> & contour_nodes_subfacets, const std::set & surface_nodes, Real av_stress_threshold); /// find critical facets by looping on contour nodes std::map> findStressedFacetsLoops( const std::map & contour_subfacets_coh_el, const std::map & surface_subfacets_crack_nb, const std::map> & contour_nodes_subfacets, const std::set & surface_nodes, Real min_dot); /// average integrated eff stress over the area of the loop Real averageEffectiveStressInMultipleFacets(Array & facet_loop); void decomposeFacetLoopPerMaterial( const Array & facet_loop, UInt crack_nb, std::map> & mat_index_facet_nbs_crack_nbs); /// apply eigen opening at all cohesives (including ghosts) void applyEigenOpening(Real eigen_strain); /// distribute total ASR gel volume along existing crack surface void applyEigenOpeningGelVolumeBased(Real gel_volume_ratio); /// distribute total ASR gel volume along initial crack surface void applyEigenOpeningToInitialCrack(Real gel_volume_ratio, const Array> & ASR_facets); + void applyPointForceToAsrCentralNodes(Real force_norm); + /// outputs crack area, volume into a file void outputCrackData(std::ofstream & file_output, Real time); /// computes crack area and volume per material and ASR volume std::tuple computeCrackData(const ID & material_name); + /// computes normal openings between asr central nodes and outputs them + void outputCrackOpenings(std::ofstream & file_output, Real time); + + /// splits asr central nodes into pairs and computes averaged normals + void identifyASRCentralNodePairsAndNormals(); + /// compute crack contour segments, surface segments, contour nodes and /// surface nodes std::tuple, std::map, std::set, std::map>> determineCrackSurface(); /// search for crack contour,etc. in a single facet void searchCrackSurfaceInSingleFacet( Element & facet, std::map & contour_subfacets_coh_el, std::map & surface_subfacets_crack_nb, std::set & surface_nodes, std::map> & contour_nodes_subfacets, std::set & visited_subfacets); /// update delta_max values in all sequential cohesive mats template UInt updateDeltaMax(); /// max ratio traction norm over max possible traction (SLA) template std::tuple getMaxDeltaMaxExcess(); // /// 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); AKANTU_GET_MACRO(NodesEffStress, nodes_eff_stress, const Array &); /// get the corner nodes AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(ASRFacets, asr_facets, const Array> &); AKANTU_GET_MACRO(ASRFacetsFromMesh, ASR_facets_from_mesh, const Array> &); AKANTU_GET_MACRO(ASRCentralNodes, asr_central_nodes, const Array &); + AKANTU_GET_MACRO(ASRCentralNodePairs, asr_central_node_pairs, const auto &); + AKANTU_GET_MACRO(ASRNormalsPairs, asr_normals_pairs, const auto &); /* --------------------------------------------------------------------- */ /* 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; + bool asr_central_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 partition_border_nodes; /// average eff stress on facet loops around a node Array nodes_eff_stress; /// array to store flags on nodes where ASR elements are inserted Array ASR_nodes; /// Vector to store the initially inserted facets per asr site Array> asr_facets; Array> ASR_facets_from_mesh; + /// duplicated nodes. first 1/2 - initialnodes, second - duplicated Array asr_central_nodes; + /// asr central nodes devided into pairs + Array> asr_central_node_pairs; + /// normals per node averaged per asr_facets of each site + Array, Vector>> asr_normals_pairs; }; /* ------------------------------------------------------------------ */ /* 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_ratio, const std::string group_name, const Real 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(); auto && fem_boundary = model.getFEEngineBoundary(); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(type_facet, gt); 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"); // get normal to the current positions const auto & current_pos = model.getCurrentPosition(); Array quad_normals(0, dim); fem_boundary.computeNormalsOnIntegrationPoints(current_pos, quad_normals, type_facet, gt); auto normals_it = quad_normals.begin(dim); Vector normal_corrected( normals_it[quad_point.element * nb_quad_points + quad_point.num_point]); // 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 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---M---B \ | D--N--C */ UInt A, B, C, D; A = first_facet_nodes(0); B = first_facet_nodes(1); 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]); 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{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_ratio; const std::string group_name; const Real compressibility; }; /* ------------------------------------------------------------------- */ class PressureVolumeDependent3D : public BC::Neumann::NeumannFunctor { public: PressureVolumeDependent3D(SolidMechanicsModelCohesive & model, const Real ASR_volume_ratio, const Real compressibility, const Array> ASR_facets_from_mesh) : model(model), ASR_volume_ratio(ASR_volume_ratio), compressibility(compressibility), ASR_facets_from_mesh(ASR_facets_from_mesh) {} inline void operator()(const IntegrationPoint & quad_point, Vector & dual, const Vector & /*coord*/, const Vector & /*normal*/) const { // get element types auto && mesh = model.getMesh(); const FEEngine & fe_engine = model.getFEEngine("CohesiveFEEngine"); auto dim = mesh.getSpatialDimension(); Element facet{quad_point.type, quad_point.element, quad_point.ghost_type}; ElementType type_cohesive = FEEngine::getCohesiveElementType(facet.type); auto nb_quad_coh = fe_engine.getNbIntegrationPoints(type_cohesive); auto && facet_nodes = mesh.getConnectivity(facet); auto && fem_boundary = model.getFEEngineBoundary(); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(facet.type, facet.ghost_type); // get normal to the current positions const auto & current_pos = model.getCurrentPosition(); Array quad_normals(0, dim); fem_boundary.computeNormalsOnIntegrationPoints( current_pos, quad_normals, facet.type, facet.ghost_type); auto normals_it = quad_normals.begin(dim); Vector normal( normals_it[quad_point.element * nb_quad_points + quad_point.num_point]); // search for this facet in the ASR facets array UInt ASR_site_nb(-1); UInt id_in_array(-1); for (auto && data : enumerate(this->ASR_facets_from_mesh)) { auto & one_site_facets = std::get<1>(data); id_in_array = one_site_facets.find(facet); if (id_in_array != UInt(-1)) { ASR_site_nb = std::get<0>(data); break; } } AKANTU_DEBUG_ASSERT(ASR_site_nb != UInt(-1), "Quad point doesn't belong to the ASR facets"); // find connected solid element auto && facet_neighbors = mesh.getElementToSubelement(facet); Element solid_el{ElementNull}; for (auto && neighbor : facet_neighbors) { if (mesh.getKind(neighbor.type) == _ek_regular) { solid_el = neighbor; break; } } AKANTU_DEBUG_ASSERT(solid_el != ElementNull, "Couldn't identify neighboring solid el"); // find the out-of-plane solid node auto && solid_nodes = mesh.getConnectivity(solid_el); UInt out_node(-1); for (auto && solid_node : solid_nodes) { auto ret = std::find(facet_nodes.begin(), facet_nodes.end(), solid_node); if (ret == facet_nodes.end()) { out_node = solid_node; break; } } AKANTU_DEBUG_ASSERT(out_node != UInt(-1), "Couldn't identify out-of-plane node"); // build the reference vector Vector ref_vector(dim); mesh.getBarycenter(facet, ref_vector); Vector pos(mesh.getNodes().begin(dim)[out_node]); ref_vector = pos - ref_vector; // check if ref vector and the normal are in the same half space if (ref_vector.dot(normal) < 0) { normal *= -1; } // compute volume (area * normal_opening) of current ASR site // form cohesive element filter from the 1st half of facet filter std::set site_cohesives; for (auto & ASR_facet : this->ASR_facets_from_mesh(ASR_site_nb)) { if (ASR_facet.type == facet.type and ASR_facet.ghost_type == facet.ghost_type) { // find connected cohesive auto & connected_els = mesh.getElementToSubelement(ASR_facet); for (auto & connected_el : connected_els) { if (connected_el.type == type_cohesive) { site_cohesives.emplace(connected_el); break; } } } } // integrate normal opening over identified element filter Real site_volume{0}; Real site_area{0}; const Array & material_index_vec = model.getMaterialByElement(type_cohesive, facet.ghost_type); const Array & material_local_numbering_vec = model.getMaterialLocalNumbering(type_cohesive, facet.ghost_type); for (auto & coh_el : site_cohesives) { Material & material = model.getMaterial(material_index_vec(coh_el.element)); UInt material_local_num = material_local_numbering_vec(coh_el.element); Array single_el_array; single_el_array.push_back(coh_el.element); auto & opening_norm_array = material.getInternal( "normal_opening_norm")(coh_el.type, coh_el.ghost_type); Array opening_norm_el; for (UInt i = 0; i != nb_quad_coh; i++) { auto & opening_per_quad = opening_norm_array(material_local_num * nb_quad_coh + i); opening_norm_el.push_back(opening_per_quad); } site_volume += fe_engine.integrate(opening_norm_el, coh_el.type, coh_el.ghost_type, single_el_array); Array area(fe_engine.getNbIntegrationPoints(type_cohesive), 1, 1.); site_area += fe_engine.integrate(area, coh_el.type, coh_el.ghost_type, single_el_array); } Real ASR_volume = site_area * ASR_volume_ratio; Real volume_change = site_volume - ASR_volume; Real pressure_change{0}; if (volume_change < 0) { pressure_change = -volume_change / ASR_volume / this->compressibility; } std::cout << "ASR volume = " << site_area << " x " << ASR_volume_ratio << " = " << ASR_volume << " site volume " << site_volume << " volume change " << volume_change << " pressure delta " << pressure_change << std::endl; dual = pressure_change * normal; } protected: SolidMechanicsModelCohesive & model; const Real ASR_volume_ratio; const Real compressibility; const Array> ASR_facets_from_mesh; }; /* ------------------------------------------------------------------------ */ 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; }; +/* ------------------------------------------------------------------ */ +/* solid mechanics model cohesive + delta d at the ASR nodes */ +/* ------------------------------------------------------------------ */ +class SolidMechanicsModelCohesiveDelta : public SolidMechanicsModelCohesive { +public: + SolidMechanicsModelCohesiveDelta(Mesh & mesh) + : SolidMechanicsModelCohesive(mesh), mesh(mesh), delta_u(0.) {} + + /* ------------------------------------------------------------------*/ + void assembleInternalForces() { + // displacement correction + applyDisplacementDifference(); + + SolidMechanicsModelCohesive::assembleInternalForces(); + } + /* ----------------------------------------------------------------- */ + void applyDisplacementDifference() { + auto dim = this->mesh.getSpatialDimension(); + auto & disp = this->getDisplacement(); + auto & boun = this->getBlockedDOFs(); + + // get normal to the initial positions + auto it_disp = make_view(disp, dim).begin(); + auto it_boun = make_view(boun, dim).begin(); + + for (auto && data : zip(asr_central_node_pairs, asr_normals_pairs)) { + auto && node_pair = std::get<0>(data); + auto && normals_pair = std::get<1>(data); + + auto node1 = node_pair.first; + auto node2 = node_pair.second; + auto normal1 = normals_pair.first; + auto normal2 = normals_pair.second; + Vector displ1 = it_disp[node1]; + Vector displ2 = it_disp[node2]; + if (mesh.isPeriodicSlave(node1)) { + displ1.copy(displ2); + } else { + displ2.copy(displ1); + } + displ1 += normal1 * this->delta_u / 2.; + displ2 += normal2 * this->delta_u / 2.; + // Vector dof = it_boun[node]; + // dof.set(true); + } + } + /* ------------------------------------------------------------------*/ +public: + // Acessors + AKANTU_GET_MACRO_NOT_CONST(DeltaU, delta_u, Real &); + AKANTU_SET_MACRO(ASRNodePairs, asr_central_node_pairs, auto); + AKANTU_SET_MACRO(ASRNormalsPairs, asr_normals_pairs, auto); + +protected: + Mesh & mesh; + Real delta_u; + Array> asr_central_node_pairs; + Array, Vector>> asr_normals_pairs; +}; +/* ------------------------------------------------------------------ */ + } // namespace akantu #endif /* __AKANTU_ASR_TOOLS_HH__ */ diff --git a/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.hh b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.hh index 99b984b50..ca9a2be76 100644 --- a/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.hh +++ b/extra_packages/asr-tools/src/mesh_utils/crack_numbers_updater_inline_impl.hh @@ -1,98 +1,101 @@ /* -------------------------------------------------------------------------- */ #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::_crack_nb) { size += sizeof(int); 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(); size += nb_elements * sizeof(UInt); } } } return size; } /* -------------------------------------------------------------------------- */ inline void CrackNumbersUpdater::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { - if (tag != SynchronizationTag::_crack_nb) + if (tag != SynchronizationTag::_crack_nb) { return; + } auto & mesh = this->model.getMesh(); auto && comm = mesh.getCommunicator(); int prank = comm.whoAmI(); buffer << prank; 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)) { auto & coh_el_ids = elements_range.getElements(); for (UInt coh_el : coh_el_ids) { auto & crack_numbers = model.getMesh().getData("crack_numbers", type); buffer << crack_numbers(coh_el); } } } } /* -------------------------------------------------------------------------- */ inline void CrackNumbersUpdater::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { - if (tag != SynchronizationTag::_crack_nb) + if (tag != SynchronizationTag::_crack_nb) { 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)) { if (not mesh.hasData("crack_numbers", type, gt)) { mesh.getDataPointer("crack_numbers", type, gt); } auto & crack_numbers = mesh.getData("crack_numbers", type, gt); // // resize crack_nb array if size of ghost coh elements is > // if (crack_numbers.size() != ghost_coh_el_ids.size()) // crack_numbers.resize(ghost_coh_el_ids.size()); auto & ghost_coh_el_ids = elements_range.getElements(); for (UInt ghost_coh_el : ghost_coh_el_ids) { // resize crack_nb array on each iteration - if (crack_numbers.size() <= ghost_coh_el) + if (crack_numbers.size() <= ghost_coh_el) { crack_numbers.resize(ghost_coh_el + 1); + } buffer >> crack_numbers(ghost_coh_el); } } } } } // namespace akantu #endif /* __AKANTU_CRACK_NUMBERS_UPDATER_INLINE_IMPL_CC__ */ diff --git a/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc b/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc index c56e7ce53..38ca2d40e 100644 --- a/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc +++ b/extra_packages/asr-tools/src/solid_mechanics_model_RVE.cc @@ -1,273 +1,275 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" #include "aka_random_generator.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "non_local_manager.hh" #include "parser.hh" #include "sparse_matrix.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), ASRTools(dynamic_cast(*this)), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), stiffness_changed(true) { AKANTU_DEBUG_IN(); RandomGenerator::seed(1); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); mesh.makePeriodic(_x, "left", "right"); mesh.makePeriodic(_y, "bottom", "top"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 50); solver.set("threshold", 1e-4); solver.set("convergence_type", SolveConvergenceCriteria::_solution); // this->initMaterials(); /// compute the volume of the RVE this->computeModelVolume(); std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id; // << this->memory_id - 1; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); this->addDumpField("delta_T"); this->addDumpField("extra_volume"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::assembleInternalForces() { /// displacement correction auto disp_begin = make_view(*this->displacement, spatial_dimension).begin(); auto u_top_corr = Vector(disp_begin[this->corner_nodes(2)]) - Vector(disp_begin[this->corner_nodes(1)]); auto u_right_corr = Vector(disp_begin[this->corner_nodes(2)]) - Vector(disp_begin[this->corner_nodes(3)]); auto correct_disp = [&](auto && node, auto && correction) { const auto & pair = mesh.getPeriodicMasterSlaves().equal_range(node); for (auto && data : range(pair.first, pair.second)) { auto slave = data.second; auto slave_disp = Vector(disp_begin[slave]); slave_disp = Vector(disp_begin[node]); slave_disp += correction; } }; for (auto node : bottom_nodes) { correct_disp(node, u_top_corr); } for (auto node : left_nodes) { correct_disp(node, u_right_corr); } SolidMechanicsModel::assembleInternalForces(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); - for (; prestrain_it != prestrain_end; ++prestrain_it) + for (; prestrain_it != prestrain_end; ++prestrain_it) { (*prestrain_it) = prestrain; + } } /// advance the damage MaterialDamageIterativeInterface & mat_paste = dynamic_cast( this->getMaterial("paste")); MaterialDamageIterativeInterface & mat_aggregate = dynamic_cast( this->getMaterial("aggregate")); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; this->stiffness_changed = false; /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); /// save nodal fields (disp, boun & ext_force) this->storeNodalFields(); do { /// restore nodals and update grad_u accordingly this->restoreNodalFields(); /// restore historical internal fields (sigma_v for visc) this->restoreInternalFields(); this->solveStep(); /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; - if (max_eq_stress_aggregate > max_eq_stress_paste) + if (max_eq_stress_aggregate > max_eq_stress_paste) { nb_damaged_elements = mat_aggregate.updateDamage(); - else if (max_eq_stress_aggregate < max_eq_stress_paste) + } else if (max_eq_stress_aggregate < max_eq_stress_paste) { nb_damaged_elements = mat_paste.updateDamage(); - else + } else { nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); - + } /// mark the flag to update stiffness if elements were damaged if (nb_damaged_elements) { this->stiffness_changed = true; } std::cout << "Proc " << prank << " the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); // if (this->nb_dumps % 10 == 0) { // this->dump(); // } // this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated - if (!are_materials_instantiated) + if (!are_materials_instantiated) { instantiateMaterials(); + } if (use_RVE_mat_selector) { auto tmp = std::make_shared( *this, "gel", this->nb_gel_pockets, "aggregate"); tmp->setFallback(material_selector); material_selector = tmp; } this->assignMaterialToElements(); // synchronize the element material arrays this->synchronize(SynchronizationTag::_material_id); for (auto & material : materials) { /// init internals properties const auto type = material->getID(); if (type.find("material_FE2") != std::string::npos) continue; material->initMaterial(); } this->synchronize(SynchronizationTag::_smm_init_mat); if (this->non_local_manager) { this->non_local_manager->initialize(); } // SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } } // namespace akantu