diff --git a/extra_packages/asr-tools/src/asr_tools.cc b/extra_packages/asr-tools/src/asr_tools.cc index bf0b134ff..80f2de069 100644 --- a/extra_packages/asr-tools/src/asr_tools.cc +++ b/extra_packages/asr-tools/src/asr_tools.cc @@ -1,2753 +1,2752 @@ /** * @file ASR_tools.cc * @author Aurelia Cuba Ramos * @author Emil Gallyamov * * @brief implementation tools for the analysis of ASR samples * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include /* -------------------------------------------------------------------------- */ #include "aka_voigthelper.hh" #include "asr_tools.hh" #include "communicator.hh" #include "material_FE2.hh" #include "material_damage_iterative_orthotropic.hh" #include "material_iterative_stiffness_reduction.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" #include "solid_mechanics_model_RVE.hh" #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), disp_stored(0, model.getSpatialDimension()), ext_force_stored(0, model.getSpatialDimension()), boun_stored(0, model.getSpatialDimension()), tensile_homogenization(false) { // register event handler for asr tools auto & mesh = model.getMesh(); // auto const dim = model.getSpatialDimension(); mesh.registerEventHandler(*this, _ehp_lowest); if (mesh.hasMeshFacets()) mesh.getMeshFacets().registerEventHandler(*this, _ehp_lowest); /// find four corner nodes of the RVE findCornerNodes(); } /* -------------------------------------------------------------------------- */ void ASRTools::computePhaseVolumes() { /// compute volume of each phase and save it into a map for (auto && mat : model.getMaterials()) { this->phase_volumes[mat.getName()] = computePhaseVolume(mat.getName()); auto it = this->phase_volumes.find(mat.getName()); if (it == this->phase_volumes.end()) this->phase_volumes.erase(mat.getName()); } } /* -------------------------------------------------------------------------- */ void ASRTools::computeModelVolume() { auto const dim = model.getSpatialDimension(); auto & mesh = model.getMesh(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); GhostType gt = _not_ghost; this->volume = 0; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { Array Volume(mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume += fem.integrate(Volume, element_type); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(this->volume, SynchronizerOperation::_sum); } } /* ------------------------------------------------------------------------- */ void ASRTools::applyFreeExpansionBC() { /// boundary conditions const auto & mesh = model.getMesh(); const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); const UInt dim = mesh.getSpatialDimension(); const Array & pos = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & boun = model.getBlockedDOFs(); /// accessing bounds Real bottom = lowerBounds(1); Real left = lowerBounds(0); Real top = upperBounds(1); Real eps = std::abs((top - bottom) * 1e-6); switch (dim) { case 2: { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; if (std::abs(pos(i, 0) - left) < eps) { boun(i, 0) = true; disp(i, 0) = 0.0; } } } break; } case 3: { /// accessing bounds Real back = lowerBounds(2); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; if ((std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 2) = 0.0; } } } break; } } } /* -------------------------------------------------------------------------- */ void ASRTools::applyLoadedBC(const Vector & traction, const ID & element_group, bool multi_axial) { /// boundary conditions const auto & mesh = model.getMesh(); const UInt dim = mesh.getSpatialDimension(); const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real left = lowerBounds(0); Real top = upperBounds(1); Real eps = std::abs((top - bottom) * 1e-6); const Array & pos = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & boun = model.getBlockedDOFs(); // disp.clear(); // boun.clear(); switch (dim) { case 2: { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; } if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } if (multi_axial && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } } break; } case 3: { Real back = lowerBounds(2); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 1) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; disp(i, 2) = 0.0; } if (std::abs(pos(i, 1) - bottom) < eps) { boun(i, 1) = true; disp(i, 1) = 0.0; } if (multi_axial && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 2) = true; disp(i, 2) = 0.0; } if (multi_axial && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; disp(i, 0) = 0.0; } } } } // try { model.applyBC(BC::Neumann::FromTraction(traction), element_group); // } catch (...) { // } } /* -------------------------------------------------------------------------- */ void ASRTools::fillNodeGroup(NodeGroup & node_group, bool multi_axial) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); const Vector & upperBounds = mesh.getUpperBounds(); const Vector & lowerBounds = mesh.getLowerBounds(); Real top = upperBounds(1); Real right = upperBounds(0); Real bottom = lowerBounds(1); Real front; if (dim == 3) front = upperBounds(2); Real eps = std::abs((top - bottom) * 1e-6); const Array & pos = mesh.getNodes(); if (multi_axial) { /// fill the NodeGroup with the nodes on the left, bottom and back /// surface for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 0) - right) < eps) node_group.add(i); if (std::abs(pos(i, 1) - top) < eps) node_group.add(i); if (dim == 3 && std::abs(pos(i, 2) - front) < eps) node_group.add(i); } } /// fill the NodeGroup with the nodes on the bottom surface else { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - top) < eps) node_group.add(i); } } } /* -------------------------------------------------------------------------- */ Real ASRTools::computeAverageDisplacement(SpatialDirection direction) { Real av_displ = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); const Vector & upperBounds = mesh.getUpperBounds(); const Vector & lowerBounds = mesh.getLowerBounds(); Real right = upperBounds(0); Real left = lowerBounds(0); Real top = upperBounds(1); Real front; if (dim == 3) front = upperBounds(2); Real eps = std::abs((right - left) * 1e-6); const Array & pos = mesh.getNodes(); const Array & disp = model.getDisplacement(); UInt nb_nodes = 0; if (direction == _x) { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 0) - right) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 0); ++nb_nodes; } } } else if (direction == _y) { for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 1) - top) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 1); ++nb_nodes; } } } else if ((direction == _z) && (model.getSpatialDimension() == 3)) { AKANTU_DEBUG_ASSERT(model.getSpatialDimension() == 3, "no z-direction in 2D problem"); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { if (std::abs(pos(i, 2) - front) < eps && mesh.isLocalOrMasterNode(i)) { av_displ += disp(i, 2); ++nb_nodes; } } } else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(av_displ, SynchronizerOperation::_sum); comm.allReduce(nb_nodes, SynchronizerOperation::_sum); } return av_displ / nb_nodes; } /* -------------------------------------------------------------------------- */ Real ASRTools::computeVolumetricExpansion(SpatialDirection direction) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); UInt tensor_component = 0; if (dim == 2) { if (direction == _x) tensor_component = 0; else if (direction == _y) tensor_component = 3; else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); } else if (dim == 3) { if (direction == _x) tensor_component = 0; else if (direction == _y) tensor_component = 4; else if (direction == _z) tensor_component = 8; else AKANTU_EXCEPTION("The parameter for the testing direction is wrong!!!"); } else AKANTU_EXCEPTION("This example does not work for 1D!!!"); Real gradu_tot = 0; Real tot_volume = 0; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_regular)) { const FEEngine & fe_engine = model.getFEEngine(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { const ElementTypeMapUInt & element_filter_map = model.getMaterial(m).getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = model.getMaterial(m).getElementFilter(element_type); if (!elem_filter.size()) continue; // const Array & gradu_vec = // model.getMaterial(m).getGradU(element_type); Array gradu_vec(elem_filter.size() * fe_engine.getNbIntegrationPoints(element_type), dim * dim, 0.); fe_engine.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vec, dim, element_type, gt, elem_filter); Array int_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fe_engine.integrate(gradu_vec, int_gradu_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) { gradu_tot += int_gradu_vec(k, tensor_component); } } Array Volume(mesh.getNbElement(element_type) * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); Real int_volume = fe_engine.integrate(Volume, element_type); tot_volume += int_volume; } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(gradu_tot, SynchronizerOperation::_sum); comm.allReduce(tot_volume, SynchronizerOperation::_sum); } return gradu_tot / tot_volume; } /* -------------------------------------------------------------------------- */ Real ASRTools::computeDamagedVolume(const ID & mat_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Real total_damage = 0; GhostType gt = _not_ghost; const Material & mat = model.getMaterial(mat_name); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const FEEngine & fe_engine = model.getFEEngine(); const ElementTypeMapUInt & element_filter_map = mat.getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = mat.getElementFilter(element_type); if (!elem_filter.size()) continue; const Array & damage = mat.getInternal("damage")(element_type); total_damage += fe_engine.integrate(damage, element_type, gt, elem_filter); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(total_damage, SynchronizerOperation::_sum); } return total_damage; } /* -------------------------------------------------------------------------- */ void ASRTools::computeStiffnessReduction(std::ofstream & file_output, Real time, bool tension) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); /// save nodal values before test storeNodalFields(); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { Real int_residual_x = 0; Real int_residual_y = 0; // try { int_residual_x = performLoadingTest(_x, tension); int_residual_y = performLoadingTest(_y, tension); // } catch (...) { // } if (prank == 0) file_output << time << "," << int_residual_x << "," << int_residual_y << std::endl; } else { Real int_residual_x = performLoadingTest(_x, tension); Real int_residual_y = performLoadingTest(_y, tension); Real int_residual_z = performLoadingTest(_z, tension); if (prank == 0) file_output << time << "," << int_residual_x << "," << int_residual_y << "," << int_residual_z << std::endl; } /// return the nodal values restoreNodalFields(); } /* -------------------------------------------------------------------------- */ void ASRTools::storeNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); this->disp_stored.copy(disp); this->boun_stored.copy(boun); this->ext_force_stored.copy(ext_force); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); disp.copy(this->disp_stored); boun.copy(this->boun_stored); ext_force.copy(this->ext_force_stored); /// update grad_u // model.assembleInternalForces(); } /* ---------------------------------------------------------------------- */ void ASRTools::resetNodalFields() { auto & disp = this->model.getDisplacement(); auto & boun = this->model.getBlockedDOFs(); auto & ext_force = this->model.getExternalForce(); disp.clear(); boun.clear(); ext_force.clear(); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreInternalFields() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); mat.restorePreviousState(); } } /* -------------------------------------------------------------------------- */ void ASRTools::resetInternalFields() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); mat.resetInternalsWithHistory(); } } /* -------------------------------------------------------------------------- */ void ASRTools::storeDamageField() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.isInternal("damage_stored", _ek_regular) && mat.isInternal("reduction_step_stored", _ek_regular)) { for (auto & el_type : mat.getElementFilter().elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto & red_stored = mat.getInternal("reduction_step_stored")(el_type, _not_ghost); auto & red_current = mat.getInternal("reduction_step")(el_type, _not_ghost); red_stored.copy(red_current); auto & dam_stored = mat.getInternal("damage_stored")(el_type, _not_ghost); auto & dam_current = mat.getInternal("damage")(el_type, _not_ghost); dam_stored.copy(dam_current); } } } } /* -------------------------------------------------------------------------- */ void ASRTools::restoreDamageField() { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.isInternal("damage_stored", _ek_regular) && mat.isInternal("reduction_step_stored", _ek_regular)) { for (auto & el_type : mat.getElementFilter().elementTypes( _all_dimensions, _not_ghost, _ek_not_defined)) { auto & red_stored = mat.getInternal("reduction_step_stored")(el_type, _not_ghost); auto & red_current = mat.getInternal("reduction_step")(el_type, _not_ghost); red_current.copy(red_stored); auto & dam_stored = mat.getInternal("damage_stored")(el_type, _not_ghost); auto & dam_current = mat.getInternal("damage")(el_type, _not_ghost); dam_current.copy(dam_stored); } } } } /* -------------------------------------------------------------------------- */ Real ASRTools::performLoadingTest(SpatialDirection direction, bool tension) { UInt dir; if (direction == _x) dir = 0; if (direction == _y) dir = 1; if (direction == _z) { AKANTU_DEBUG_ASSERT(model.getSpatialDimension() == 3, "Error in problem dimension!!!"); dir = 2; } const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); /// indicate to material that its in the loading test UInt nb_materials = model.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", true); } } /// boundary conditions const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real left = lowerBounds(0); Real back; if (dim == 3) back = lowerBounds(2); Real eps = std::abs((top - bottom) * 1e-6); Real imposed_displacement = std::abs(lowerBounds(dir) - upperBounds(dir)); imposed_displacement *= 0.001; const auto & pos = mesh.getNodes(); auto & disp = model.getDisplacement(); auto & boun = model.getBlockedDOFs(); auto & ext_force = model.getExternalForce(); UInt nb_nodes = mesh.getNbNodes(); disp.clear(); boun.clear(); ext_force.clear(); if (dim == 3) { for (UInt i = 0; i < nb_nodes; ++i) { /// fix one corner node to avoid sliding if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps) && (std::abs(pos(i, 2) - back) < eps)) { boun(i, 0) = true; boun(i, 1) = true; boun(i, 2) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; disp(i, 2) = 0.0; } if ((std::abs(pos(i, dir) - lowerBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = 0.0; } if ((std::abs(pos(i, dir) - upperBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = (2 * tension - 1) * imposed_displacement; } } } else { for (UInt i = 0; i < nb_nodes; ++i) { /// fix one corner node to avoid sliding if ((std::abs(pos(i, 1) - bottom) < eps) && (std::abs(pos(i, 0) - left) < eps)) { boun(i, 0) = true; boun(i, 1) = true; disp(i, 0) = 0.0; disp(i, 1) = 0.0; } if ((std::abs(pos(i, dir) - lowerBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = 0.0; } if ((std::abs(pos(i, dir) - upperBounds(dir)) < eps)) { boun(i, dir) = true; disp(i, dir) = (2 * tension - 1) * imposed_displacement; } } } try { model.solveStep(); } catch (debug::Exception & e) { auto & solver = model.getNonLinearSolver("static"); int nb_iter = solver.get("nb_iterations"); std::cout << "Loading test did not converge in " << nb_iter << " iterations." << std::endl; throw e; } /// compute the force (residual in this case) along the edge of the /// imposed displacement Real int_residual = 0.; const Array & residual = model.getInternalForce(); for (UInt n = 0; n < nb_nodes; ++n) { if (std::abs(pos(n, dir) - upperBounds(dir)) < eps && mesh.isLocalOrMasterNode(n)) { int_residual += -residual(n, dir); } } /// indicate to material that its out of the loading test for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", false); } } /// restore historical internal fields restoreInternalFields(); /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(int_residual, SynchronizerOperation::_sum); } return int_residual; } /* -------------------------------------------------------------------------- */ void ASRTools::computeAverageProperties(std::ofstream & file_output) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real damage_agg, damage_paste, crack_agg, crack_paste; computeDamageRatioPerMaterial(damage_agg, "aggregate"); computeDamageRatioPerMaterial(damage_paste, "paste"); computeCrackVolumePerMaterial(crack_agg, "aggregate"); computeCrackVolumePerMaterial(crack_paste, "paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << damage_agg << "," << damage_paste << "," << crack_agg << "," << crack_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << "," << damage_agg << "," << damage_paste << "," << crack_agg << "," << crack_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAverageProperties(std::ofstream & file_output, Real time) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real damage_agg, damage_paste, crack_agg, crack_paste; computeDamageRatioPerMaterial(damage_agg, "aggregate"); computeDamageRatioPerMaterial(damage_paste, "paste"); computeCrackVolumePerMaterial(crack_agg, "aggregate"); computeCrackVolumePerMaterial(crack_paste, "paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << damage_agg << "," << damage_paste << "," << crack_agg << "," << crack_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << "," << damage_agg << "," << damage_paste << "," << crack_agg << "," << crack_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAveragePropertiesCohesiveModel( std::ofstream & file_output, Real time) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::computeAveragePropertiesFe2Material(std::ofstream & file_output, Real time) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim != 1, "Example does not work for 1D"); Real av_strain_x = computeVolumetricExpansion(_x); Real av_strain_y = computeVolumetricExpansion(_y); Real av_displ_x = computeAverageDisplacement(_x); Real av_displ_y = computeAverageDisplacement(_y); Real crack_agg = averageScalarField("crack_volume_ratio_agg"); Real crack_paste = averageScalarField("crack_volume_ratio_paste"); auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); if (dim == 2) { if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_displ_x << "," << av_displ_y << "," << crack_agg << "," << crack_paste << std::endl; } else { Real av_displ_z = computeAverageDisplacement(_z); Real av_strain_z = computeVolumetricExpansion(_z); if (prank == 0) file_output << time << "," << av_strain_x << "," << av_strain_y << "," << av_strain_z << "," << av_displ_x << "," << av_displ_y << "," << av_displ_z << "," << crack_agg << "," << crack_paste << std::endl; } } /* -------------------------------------------------------------------------- */ void ASRTools::saveState(UInt current_step) { /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); GhostType gt = _not_ghost; /// save the reduction step on each quad UInt nb_materials = model.getNbMaterials(); /// open the output file std::stringstream file_stream; file_stream << "./restart/reduction_steps_file_" << prank << "_" << current_step << ".txt"; std::string reduction_steps_file = file_stream.str(); std::ofstream file_output; file_output.open(reduction_steps_file.c_str()); if (!file_output.is_open()) AKANTU_EXCEPTION("Could not create the file " + reduction_steps_file + ", does its folder exist?"); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = model.getMaterial(m); if (mat.getName() == "gel") continue; /// get the reduction steps internal field const InternalField & reduction_steps = mat.getInternal("damage_step"); /// loop over all types in that material for (auto element_type : reduction_steps.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; const Array & reduction_step_array = reduction_steps(element_type, gt); for (UInt i = 0; i < reduction_step_array.size(); ++i) { file_output << reduction_step_array(i) << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + reduction_steps_file); } } } /// close the file file_output.close(); comm.barrier(); /// write the number of the last successfully saved step in a file if (prank == 0) { std::string current_step_file = "./restart/current_step.txt"; std::ofstream file_output; file_output.open(current_step_file.c_str()); file_output << current_step << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + current_step_file); file_output.close(); } } /* -------------------------------------------------------------------------- */ bool ASRTools::loadState(UInt & current_step) { current_step = 0; bool restart = false; /// variables for parallel execution auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); GhostType gt = _not_ghost; /// proc 0 has to save the current step if (prank == 0) { std::string line; std::string current_step_file = "./restart/current_step.txt"; std::ifstream file_input; file_input.open(current_step_file.c_str()); if (file_input.is_open()) { std::getline(file_input, line); std::stringstream sstr(line); sstr >> current_step; file_input.close(); } } comm.broadcast(current_step); if (!current_step) return restart; if (prank == 0) std::cout << "....Restarting simulation" << std::endl; restart = true; /// save the reduction step on each quad UInt nb_materials = model.getNbMaterials(); /// open the output file std::stringstream file_stream; file_stream << "./restart/reduction_steps_file_" << prank << "_" << current_step << ".txt"; std::string reduction_steps_file = file_stream.str(); std::ifstream file_input; file_input.open(reduction_steps_file.c_str()); if (!file_input.is_open()) AKANTU_EXCEPTION("Could not open file " + reduction_steps_file); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = model.getMaterial(m); if (mat.getName() == "gel") continue; /// get the material parameters Real E = mat.getParam("E"); Real max_damage = mat.getParam("max_damage"); Real max_reductions = mat.getParam("max_reductions"); /// get the internal field that need to be set InternalField & reduction_steps = const_cast &>(mat.getInternal("damage_step")); InternalField & Sc = const_cast &>(mat.getInternal("Sc")); InternalField & damage = const_cast &>(mat.getInternal("damage")); /// loop over all types in that material Real reduction_constant = mat.getParam("reduction_constant"); for (auto element_type : reduction_steps.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; Array & reduction_step_array = reduction_steps(element_type, gt); Array & Sc_array = Sc(element_type, gt); Array & damage_array = damage(element_type, gt); const Array & D_tangent = mat.getInternal("tangent")(element_type, gt); const Array & eps_u = mat.getInternal("ultimate_strain")(element_type, gt); std::string line; for (UInt i = 0; i < reduction_step_array.size(); ++i) { std::getline(file_input, line); if (file_input.fail()) AKANTU_EXCEPTION("Could not read data from file " + reduction_steps_file); std::stringstream sstr(line); sstr >> reduction_step_array(i); if (reduction_step_array(i) == 0) continue; if (reduction_step_array(i) == max_reductions) { damage_array(i) = max_damage; Real previous_damage = 1. - (1. / std::pow(reduction_constant, reduction_step_array(i) - 1)); Sc_array(i) = eps_u(i) * (1. - previous_damage) * E * D_tangent(i) / ((1. - previous_damage) * E + D_tangent(i)); } else { damage_array(i) = 1. - (1. / std::pow(reduction_constant, reduction_step_array(i))); Sc_array(i) = eps_u(i) * (1. - damage_array(i)) * E * D_tangent(i) / ((1. - damage_array(i)) * E + D_tangent(i)); } } } } /// close the file file_input.close(); return restart; } /* -------------------------------------------------------------------------- */ Real ASRTools::computePhaseVolume(const ID & mat_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Real total_volume = 0; GhostType gt = _not_ghost; const Material & mat = model.getMaterial(mat_name); for (auto element_type : mesh.elementTypes(dim, gt, _ek_regular)) { const FEEngine & fe_engine = model.getFEEngine(); const ElementTypeMapUInt & element_filter_map = mat.getElementFilter(); if (!element_filter_map.exists(element_type, gt)) continue; const Array & elem_filter = mat.getElementFilter(element_type); if (!elem_filter.size()) continue; Array volume(elem_filter.size() * fe_engine.getNbIntegrationPoints(element_type), 1, 1.); total_volume += fe_engine.integrate(volume, element_type, gt, elem_filter); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(total_volume, SynchronizerOperation::_sum); } return total_volume; } /* -------------------------------------------------------------------------- */ void ASRTools::applyEigenGradUinCracks( const Matrix prescribed_eigen_grad_u, const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); const Real max_damage = mat.getParam("max_damage"); const ElementTypeMapUInt & element_filter = mat.getElementFilter(); for (auto element_type : element_filter.elementTypes(dim, gt, _ek_not_defined)) { if (!element_filter(element_type, gt).size()) continue; const Array & damage = mat.getInternal("damage")(element_type, gt); const Real * dam_it = damage.storage(); Array & eigen_gradu = mat.getInternal("eigen_grad_u")(element_type, gt); Array::matrix_iterator eigen_it = eigen_gradu.begin(dim, dim); Array::matrix_iterator eigen_end = eigen_gradu.end(dim, dim); for (; eigen_it != eigen_end; ++eigen_it, ++dam_it) { if (Math::are_float_equal(max_damage, *dam_it)) { Matrix & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::applyEigenGradUinCracks( const Matrix prescribed_eigen_grad_u, const ElementTypeMapUInt & critical_elements, bool to_add) { GhostType gt = _not_ghost; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const Array & critical_elements_vec = critical_elements(element_type, gt); const Array & material_index_vec = model.getMaterialByElement(element_type, gt); const Array & material_local_numbering_vec = model.getMaterialLocalNumbering(element_type, gt); for (UInt e = 0; e < critical_elements_vec.size(); ++e) { UInt element = critical_elements_vec(e); Material & material = model.getMaterial(material_index_vec(element)); Array & eigen_gradu = material.getInternal("eigen_grad_u")(element_type, gt); UInt material_local_num = material_local_numbering_vec(element); Array::matrix_iterator eigen_it = eigen_gradu.begin(dim, dim); eigen_it += material_local_num; Matrix & current_eigengradu = *eigen_it; if (to_add) current_eigengradu += prescribed_eigen_grad_u; else current_eigengradu = prescribed_eigen_grad_u; } } } // /* // -------------------------------------------------------------------------- // */ void ASRTools::fillCracks(ElementTypeMapReal & saved_damage) { // const UInt dim = model.getSpatialDimension(); // const Material & mat_gel = model.getMaterial("gel"); // const Real E_gel = mat_gel.getParam("E"); // Real E_homogenized = 0.; // GhostType gt = _not_ghost; // const auto & mesh = model.getMesh(); // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // if (mat.getName() == "gel") // continue; // const Real E = mat.getParam("E"); // InternalField & damage = mat.getInternal("damage"); // for (auto element_type : mesh.elementTypes(dim, gt, _ek_regular)) { // const Array & elem_filter = // mat.getElementFilter(element_type, gt); if (!elem_filter.size()) // continue; // Array & damage_vec = damage(element_type, gt); // Array & saved_damage_vec = saved_damage(element_type, gt); // for (UInt i = 0; i < damage_vec.size(); ++i) { // saved_damage_vec(elem_filter(i)) = damage_vec(i); // E_homogenized = (E_gel - E) * damage_vec(i) + E; // damage_vec(i) = 1. - (E_homogenized / E); // } // } // } // } // /* // -------------------------------------------------------------------------- // */ void ASRTools::drainCracks(const ElementTypeMapReal & // saved_damage) { // // model.dump(); // const UInt dim = model.getSpatialDimension(); // const auto & mesh = model.getMesh(); // GhostType gt = _not_ghost; // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // if (mat.getName() == "gel") // continue; // else { // InternalField & damage = mat.getInternal("damage"); // for (auto element_type : mesh.elementTypes(dim, gt, // _ek_not_defined)) { // const Array & elem_filter = // mat.getElementFilter(element_type, gt); // if (!elem_filter.size()) // continue; // Array & damage_vec = damage(element_type, gt); // const Array & saved_damage_vec = // saved_damage(element_type, gt); for (UInt i = 0; i < // damage_vec.size(); ++i) { // damage_vec(i) = saved_damage_vec(elem_filter(i)); // } // } // } // } // // model.dump(); // } /* -------------------------------------------------------------------------- */ Real ASRTools::computeSmallestElementSize() { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // constexpr auto dim = model.getSpatialDimension(); /// compute smallest element size const Array & pos = mesh.getNodes(); Real el_h_min = std::numeric_limits::max(); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_element = mesh.getNbElement(element_type); Array X(0, nb_nodes_per_element * dim); model.getFEEngine().extractNodalToElementField(mesh, pos, X, element_type, _not_ghost); Array::matrix_iterator X_el = X.begin(dim, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el) { Real el_h = model.getFEEngine().getElementInradius(*X_el, element_type); if (el_h < el_h_min) el_h_min = el_h; } } /// find global Gauss point with highest stress /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(el_h_min, SynchronizerOperation::_min); } return el_h_min; } // /* // -------------------------------------------------------------------------- // */ // /// apply homogeneous temperature on the whole specimen // void ASRTools::applyTemperatureFieldToSolidmechanicsModel( // const Real & temperature) { // UInt dim = model.getMesh().getSpatialDimension(); // for (UInt m = 0; m < model.getNbMaterials(); ++m) { // Material & mat = model.getMaterial(m); // for (auto el_type : mat.getElementFilter().elementTypes(dim)) { // const auto & filter = mat.getElementFilter()(el_type); // if (filter.size() == 0) // continue; // auto & delta_T = mat.getArray("delta_T", el_type); // auto dt = delta_T.begin(); // auto dt_end = delta_T.end(); // for (; dt != dt_end; ++dt) { // *dt = temperature; // } // } // } // } /* -------------------------------------------------------------------------- */ void ASRTools::computeASRStrainLarive( const Real & delta_time_day, const Real & T, Real & ASRStrain, const Real & eps_inf, const Real & time_ch_ref, const Real & time_lat_ref, const Real & U_C, const Real & U_L, const Real & T_ref) { AKANTU_DEBUG_IN(); Real time_ch, time_lat, lambda, ksi, exp_ref; ksi = ASRStrain / eps_inf; if (T == 0) { ksi += 0; } else { time_ch = time_ch_ref * std::exp(U_C * (1. / T - 1. / T_ref)); time_lat = time_lat_ref * std::exp(U_L * (1. / T - 1. / T_ref)); exp_ref = std::exp(-time_lat / time_ch); lambda = (1 + exp_ref) / (ksi + exp_ref); ksi += delta_time_day / time_ch * (1 - ksi) / lambda; } ASRStrain = ksi * eps_inf; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real ASRTools::computeDeltaGelStrainThermal(const Real delta_time_day, const Real k, const Real activ_energy, const Real R, const Real T, Real & amount_reactive_particles, const Real saturation_const) { /// compute increase in gel strain value for interval of time delta_time /// as temperatures are stored in C, conversion to K is done Real delta_strain = amount_reactive_particles * k * std::exp(-activ_energy / (R * T)) * delta_time_day; amount_reactive_particles -= std::exp(-activ_energy / (R * T)) * delta_time_day / saturation_const; if (amount_reactive_particles < 0.) amount_reactive_particles = 0.; return delta_strain; } /* -----------------------------------------------------------------------*/ Real ASRTools::computeDeltaGelStrainLinear(const Real delta_time, const Real k) { /// compute increase in gel strain value for dt simply by deps = k * /// delta_time Real delta_strain = k * delta_time; return delta_strain; } /* ---------------------------------------------------------------------- */ void ASRTools::applyBoundaryConditionsRve( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // /// transform into Green strain to exclude rotations // auto F = displacement_gradient + Matrix::eye(dim); // auto E = 0.5 * (F.transpose() * F - Matrix::eye(dim)); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement /// that will be applied Vector x(dim); Vector appl_disp(dim); const auto & lower_bounds = mesh.getLowerBounds(); for (auto node : this->corner_nodes) { x(0) = pos(node, 0); x(1) = pos(node, 1); x -= lower_bounds; appl_disp.mul(displacement_gradient, x); (model.getBlockedDOFs())(node, 0) = true; (model.getDisplacement())(node, 0) = appl_disp(0); (model.getBlockedDOFs())(node, 1) = true; (model.getDisplacement())(node, 1) = appl_disp(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::applyHomogeneousTemperature(const Real & temperature) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); for (auto el_type : mat.getElementFilter().elementTypes(dim)) { const auto & filter = mat.getElementFilter()(el_type); if (filter.size() == 0) continue; auto & deltas_T = mat.getArray("delta_T", el_type); for (auto && delta_T : deltas_T) { delta_T = temperature; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::removeTemperature() { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); for (auto el_type : mat.getElementFilter().elementTypes(dim)) { const auto & filter = mat.getElementFilter()(el_type); if (filter.size() == 0) continue; auto & deltas_T = mat.getArray("delta_T", el_type); for (auto && delta_T : deltas_T) { delta_T = 0.; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::findCornerNodes() { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); switch (dim) { case 2: { corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, dim))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } break; } case 3: { corner_nodes.resize(8); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, dim))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), lower_bounds(_z))) { corner_nodes(3) = node; } // node 5 if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(4) = node; } // node 6 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(5) = node; } // node 7 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(6) = node; } // node 8 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y)) && Math::are_float_equal(X(_z), upper_bounds(_z))) { corner_nodes(7) = node; } } break; } } // AKANTU_DEBUG_ASSERT(dim == 2, "This is 2D only!"); // for (UInt i = 0; i < corner_nodes.size(); ++i) { // if (corner_nodes(i) == UInt(-1)) // AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real ASRTools::averageScalarField(const ID & field_name) { AKANTU_DEBUG_IN(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); Real average = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); if (!elem_filter.size()) continue; const auto & scalar_field = model.getMaterial(m).getInternal(field_name)(element_type); Array int_scalar_vec(elem_filter.size(), 1, "int_of_scalar"); fem.integrate(scalar_field, int_scalar_vec, 1, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_scalar_vec(k); } } /// compute total model volume if (!this->volume) computeModelVolume(); return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real ASRTools::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = model.getFEEngine("SolidMechanicsFEEngine"); Real average = 0; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & stress_vec = model.getMaterial(m).getStress(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), dim * dim, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec(k, row_index * dim + col_index); // 3 is the value // for the yy (in // 3D, the value // is 4) } } else if (field_type == "strain") { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & gradu_vec = model.getMaterial(m).getGradU(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and /// strain are equal average += 0.5 * (int_gradu_vec(k, row_index * dim + col_index) + int_gradu_vec(k, col_index * dim + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < model.getNbMaterials(); ++m) { const auto & eigen_gradu_vec = model.getMaterial(m).getInternal( "eigen_grad_u")(element_type); const auto & elem_filter = model.getMaterial(m).getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), dim * dim, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, dim * dim, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and /// strain are equal average += int_eigen_gradu_vec(k, row_index * dim + col_index); } } else { AKANTU_ERROR("Averaging not implemented for this field!!!"); } } /// compute total model volume if (!this->volume) computeModelVolume(); return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::homogenizeStressField(Matrix & stress) { AKANTU_DEBUG_IN(); stress(0, 0) = averageTensorField(0, 0, "stress"); stress(1, 1) = averageTensorField(1, 1, "stress"); stress(0, 1) = averageTensorField(0, 1, "stress"); stress(1, 0) = averageTensorField(1, 0, "stress"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::setStiffHomogenDir(Matrix & stress) { AKANTU_DEBUG_IN(); auto dim = stress.rows(); Vector eigenvalues(dim); stress.eig(eigenvalues); Real hydrostatic_stress = 0; UInt denominator = 0; for (UInt i = 0; i < dim; ++i, ++denominator) hydrostatic_stress += eigenvalues(i); hydrostatic_stress /= denominator; AKANTU_DEBUG_OUT(); this->tensile_homogenization = (hydrostatic_stress > 0); } /* -------------------------------------------------------------------------- */ void ASRTools::homogenizeStiffness(Matrix & C_macro, bool tensile_test) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim == 2, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (0.001;0;0) 2. eps_el = (0,0.001,0) 3. eps_el = (0,0,0.001) /// store and clear the eigenstrain ///(to exclude stresses due to internal pressure) Array stored_eig(0, dim * dim, 0); storeASREigenStrain(stored_eig); clearASREigenStrain(); /// save nodal values before tests storeNodalFields(); /// storage for results of 3 different loading states UInt voigt_size = 1; switch (dim) { case 2: { voigt_size = VoigtHelper<2>::size; break; } case 3: { voigt_size = VoigtHelper<3>::size; break; } } Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, // _default_value = 0); this->fillCracks(saved_damage); /// indicate to material that its in the loading test UInt nb_materials = model.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", true); } } /// virtual test 1: H(0, 0) = 0.001 * (2 * tensile_test - 1); performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.clear(); H(1, 1) = 0.001 * (2 * tensile_test - 1); performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.clear(); H(0, 1) = 0.001; H(1, 0) = 0.001; performVirtualTesting(H, stresses, strains, 2); /// indicate to material that its out of the loading test for (UInt m = 0; m < nb_materials; ++m) { Material & mat = model.getMaterial(m); if (aka::is_of_type>(mat)) { mat.setParam("loading_test", false); } } /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); /// Make C matrix symmetric Matrix C_direct(voigt_size, voigt_size); C_direct.mul(stresses, eps_inverse); for (UInt i = 0; i != voigt_size; ++i) { for (UInt j = 0; j != voigt_size; ++j) { C_macro(i, j) = 0.5 * (C_direct(i, j) + C_direct(j, i)); C_macro(j, i) = C_macro(i, j); } } /// return the nodal values and the ASR eigenstrain restoreNodalFields(); restoreASREigenStrain(stored_eig); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); auto & disp = model.getDisplacement(); auto & boun = model.getBlockedDOFs(); auto & ext_force = model.getExternalForce(); disp.clear(); boun.clear(); ext_force.clear(); applyBoundaryConditionsRve(H); model.solveStep(); /// get average stress and strain eff_stresses(0, test_no) = averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * averageTensorField(1, 0, "strain"); /// restore historical internal fields restoreInternalFields(); AKANTU_DEBUG_OUT(); } /* --------------------------------------------------- */ void ASRTools::homogenizeEigenGradU(Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = model.getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; Real E = mat.get("E"); auto & damage = mat.getInternal("damage"); GhostType gt = _not_ghost; for (auto && type : damage.elementTypes(gt)) { const auto & elem_filter = mat.getElementFilter(type); auto nb_integration_point = model.getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void ASRTools::drainCracks(const ElementTypeMapReal & saved_damage) { for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; auto & damage = mat.getInternal("damage"); GhostType gt = _not_ghost; for (auto && type : damage.elementTypes(gt)) { const auto & elem_filter = mat.getElementFilter(type); auto nb_integration_point = model.getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } /* -------------------------------------------------------------------------- */ void ASRTools::computeDamageRatio(Real & damage_ratio) { damage_ratio = 0.; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; GhostType gt = _not_ghost; const ElementTypeMapArray & filter_map = mat.getElementFilter(); // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; const FEEngine & fe_engine = model.getFEEngine(); auto & damage_array = mat.getInternal("damage")(element_type); damage_ratio += fe_engine.integrate(damage_array, element_type, gt, filter); } } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(damage_ratio, SynchronizerOperation::_sum); } /// compute total model volume if (!this->volume) computeModelVolume(); damage_ratio /= this->volume; } /* -------------------------------------------------------------------------- */ void ASRTools::computeDamageRatioPerMaterial(Real & damage_ratio, const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); const ElementTypeMapArray & filter_map = mat.getElementFilter(); const FEEngine & fe_engine = model.getFEEngine(); damage_ratio = 0.; // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; auto & damage_array = mat.getInternal("damage")(element_type); damage_ratio += fe_engine.integrate(damage_array, element_type, gt, filter); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(damage_ratio, SynchronizerOperation::_sum); } if (not this->phase_volumes.size()) computePhaseVolumes(); damage_ratio /= this->phase_volumes[material_name]; } /* -------------------------------------------------------------------------- */ void ASRTools::computeCrackVolume(Real & crack_volume_ratio) { crack_volume_ratio = 0.; const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); if (mat.getName() == "gel" || mat.getName() == "FE2_mat") continue; GhostType gt = _not_ghost; const ElementTypeMapArray & filter_map = mat.getElementFilter(); // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; auto extra_volume_copy = mat.getInternal("extra_volume")(element_type); crack_volume_ratio += Math::reduce(extra_volume_copy); } } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(crack_volume_ratio, SynchronizerOperation::_sum); } /// compute total model volume if (!this->volume) computeModelVolume(); crack_volume_ratio /= this->volume; } /* -------------------------------------------------------------------------- */ void ASRTools::computeCrackVolumePerMaterial(Real & crack_volume, const ID & material_name) { const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; Material & mat = model.getMaterial(material_name); const ElementTypeMapArray & filter_map = mat.getElementFilter(); crack_volume = 0.; // Loop over the boundary element types for (auto & element_type : filter_map.elementTypes(dim, gt)) { const Array & filter = filter_map(element_type); if (!filter_map.exists(element_type, gt)) continue; if (filter.size() == 0) continue; auto extra_volume_copy = mat.getInternal("extra_volume")(element_type); crack_volume += Math::reduce(extra_volume_copy); } /// do not communicate if model is multi-scale if (not aka::is_of_type(model)) { auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(crack_volume, SynchronizerOperation::_sum); } if (not this->phase_volumes.size()) computePhaseVolumes(); crack_volume /= this->phase_volumes[material_name]; } /* ------------------------------------------------------------------------ */ void ASRTools::dumpRve() { // if (this->nb_dumps % 10 == 0) { model.dump(); // } // this->nb_dumps += 1; } /* ------------------------------------------------------------------------ */ // void ASRTools::applyBodyForce(const Real gravity = 9.81) { // auto spatial_dimension = model.getSpatialDimension(); // model.assembleMassLumped(); // auto & mass = model.getMass(); // auto & force = model.getExternalForce(); // Vector gravity(spatial_dimension); // gravity(1) = -gravity; // for (auto && data : zip(make_view(mass, spatial_dimension), // make_view(force, spatial_dimension))) { // const auto & mass_vec = (std::get<0>(data)); // AKANTU_DEBUG_ASSERT(mass_vec.norm(), "Mass vector components are zero"); // auto & force_vec = (std::get<1>(data)); // force_vec += gravity * mass_vec; // } // } /* ------------------------------------------------------------------------ */ void ASRTools::insertCohElemByCoords(const Matrix & positions) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); const auto gt = _not_ghost; // insert cohesive elements auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh_facets = inserter.getMeshFacets(); Vector bary_facet(dim); auto & doubled_facets = mesh_facets.createElementGroup("doubled_facets", dim - 1); auto & doubled_nodes = mesh.createNodeGroup("doubled_nodes"); for (auto & position : positions) { Real min_dist = std::numeric_limits::max(); Element source_facet; for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto dist = bary_facet.distance(position); if (dist < min_dist) { min_dist = dist; source_facet = facet; } }, _spatial_dimension = dim - 1); inserter.getCheckFacets(source_facet.type, gt)(source_facet.element) = false; insertion(source_facet) = true; } inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertCohElemByLimits(const Matrix & insertion_limits, std::string coh_mat_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); auto pos_it = make_view(mesh.getNodes(), dim).begin(); // insert cohesive elements auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh_facets = inserter.getMeshFacets(); auto coh_material = model.getMaterialIndex(coh_mat_name); auto tolerance = Math::getTolerance(); Vector bary_facet(dim); for_each_element( mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); UInt coord_in_limit = 0; while (coord_in_limit < dim and bary_facet(coord_in_limit) > (insertion_limits(coord_in_limit, 0) - tolerance) and bary_facet(coord_in_limit) < (insertion_limits(coord_in_limit, 1) + tolerance)) ++coord_in_limit; if (coord_in_limit == dim) { coh_model.getFacetMaterial(facet.type, facet.ghost_type)( facet.element) = coh_material; inserter.getCheckFacets(facet.type, facet.ghost_type)(facet.element) = false; insertion(facet) = true; } }, _spatial_dimension = dim - 1); inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertCohElemRandomly(const UInt & nb_coh_elem, std::string coh_mat_name, std::string matrix_mat_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); auto & coh_model = dynamic_cast(model); // get element types const GhostType gt = _not_ghost; auto type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto type_facets = Mesh::getFacetType(type); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); Vector bary_facet(dim); auto coh_material = model.getMaterialIndex(coh_mat_name); auto & mesh_facets = inserter.getMeshFacets(); auto matrix_material_id = model.getMaterialIndex(matrix_mat_name); auto & matrix_elements = mesh_facets.createElementGroup("matrix_facets", dim - 1); for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto & facet_material = coh_model.getFacetMaterial( facet.type, facet.ghost_type)(facet.element); if (facet_material == matrix_material_id) matrix_elements.add(facet); }, _spatial_dimension = dim - 1); // Standard mersenne_twister_engine seeded with 0 UInt nb_element = matrix_elements.getElements(type_facets).size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_element - 1); for (auto _[[gnu::unused]] : arange(nb_coh_elem)) { auto id = dis(random_generator); Element facet; facet.type = type_facets; facet.element = matrix_elements.getElements(type_facets)(id); facet.ghost_type = gt; coh_model.getFacetMaterial(facet.type, facet.ghost_type)(facet.element) = coh_material; inserter.getCheckFacets(facet.type, facet.ghost_type)(facet.element) = false; insertion(facet) = true; } inserter.insertElements(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertFacetsRandomly(const UInt & nb_insertions, std::string facet_mat_name, bool add_neighbors, bool only_double_facets) { AKANTU_DEBUG_IN(); /// pick central facets (and neighbors if needed) pickFacetsRandomly(nb_insertions, facet_mat_name, add_neighbors); insertOppositeFacets(only_double_facets); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void ASRTools::insertFacetsByCoords(const Matrix & positions, bool add_neighbors, bool only_double_facets) { AKANTU_DEBUG_IN(); /// pick the facets to duplicate pickFacetsByCoord(positions, add_neighbors); insertOppositeFacets(only_double_facets); preventCohesiveInsertionInNeighbors(); AKANTU_DEBUG_OUT(); } /* ----------------------------------------------------------------------- */ void ASRTools::pickFacetsByCoord(const Matrix & positions, bool add_neighbors) { 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); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto facet_conn_it = make_view(facet_conn, nb_nodes_facet).begin(); Vector bary_facet(dim); for (auto & position : positions) { Real min_dist = std::numeric_limits::max(); Element cent_facet; cent_facet.ghost_type = gt; for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto dist = std::abs(bary_facet.distance(position)); if (dist < min_dist) { min_dist = dist; cent_facet = facet; } }, _spatial_dimension = dim - 1, _ghost_type = _not_ghost); /// eliminate posibility of a ghost neighbor bool ghost_neighbors{false}; Vector facet_nodes = facet_conn_it[cent_facet.element]; for (auto node : facet_nodes) { if (not mesh.isLocalNode(node)) ghost_neighbors = true; } if (ghost_neighbors) { std::cout << "Facet at the position " << position(0) << "," << position(1) << " is connected to a ghost element. Skipping this position" << std::endl; continue; } doubled_facets.add(cent_facet); // add all facet nodes to the group for (auto node : arange(nb_nodes_facet)) doubled_nodes.add(facet_conn(cent_facet.element, node)); if (add_neighbors) pickFacetNeighbors(cent_facet); } } /* ----------------------------------------------------------------------- */ void ASRTools::pickFacetsRandomly(UInt nb_insertions, std::string facet_mat_name, bool add_neighbors) { 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(); if (not nb_insertions) return; Vector bary_facet(dim); for_each_element(mesh_facets, [&](auto && facet) { mesh_facets.getBarycenter(facet, bary_facet); auto & facet_material = coh_model.getFacetMaterial( facet.type, facet.ghost_type)(facet.element); if (facet_material == facet_material_id) matrix_elements.add(facet); }, _spatial_dimension = dim - 1, _ghost_type = _not_ghost); UInt nb_element = matrix_elements.getElements(facet_type).size(); std::mt19937 random_generator(0); std::uniform_int_distribution<> dis(0, nb_element - 1); if (not nb_element) { auto && comm = akantu::Communicator::getWorldCommunicator(); auto prank = comm.whoAmI(); std::cout << "Proc " << prank << " couldn't place " << nb_insertions << " ASR sites" << std::endl; return; } for (UInt i = 0; i < nb_insertions; i++) { auto id = dis(random_generator); Element cent_facet; cent_facet.type = facet_type; cent_facet.element = matrix_elements.getElements(facet_type)(id); cent_facet.ghost_type = gt; /// eliminate posibility of a ghost neighbor bool ghost_neighbors{false}; Vector facet_nodes = facet_conn_it[cent_facet.element]; for (auto node : facet_nodes) { if (not mesh.isLocalNode(node)) ghost_neighbors = true; } if (ghost_neighbors) { i--; continue; } doubled_facets.add(cent_facet); /// add all facet nodes to the group for (auto node : arange(nb_nodes_facet)) { doubled_nodes.add(facet_conn(cent_facet.element, node)); } if (add_neighbors) pickFacetNeighbors(cent_facet); } } /* ----------------------------------------------------------------------- */ void ASRTools::pickFacetNeighbors(Element & cent_facet) { 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); CSR nodes_to_segments; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_segments, dim - 1); for (auto node : arange(2)) { /// vector of the central facet Vector cent_facet_dir(pos_it[facet_conn(cent_facet.element, !node)], true); cent_facet_dir -= Vector(pos_it[facet_conn(cent_facet.element, node)]); cent_facet_dir /= cent_facet_dir.norm(); Real min_dot = std::numeric_limits::max(); Element neighbor; 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; /// decide which node of the neighbor is the second UInt first_node{facet_conn(cent_facet.element, node)}; UInt second_node(-1); if (facet_conn(elem.element, 0) == first_node) { second_node = facet_conn(elem.element, 1); } else if (facet_conn(elem.element, 1) == first_node) { second_node = facet_conn(elem.element, 0); } else AKANTU_EXCEPTION( "Neighboring facet" << elem << " with nodes " << facet_conn(elem.element, 0) << " and " << facet_conn(elem.element, 1) << " doesn't have node in common with the central facet " << cent_facet << " with nodes " << facet_conn(cent_facet.element, 0) << " and " << facet_conn(cent_facet.element, 1)); // /// discard facets laying on the border of the partition // if (not mesh.isLocalNode(second_node)) { // std::cout << "Neighbor's node is on the border of the partition" // << std::endl; // continue; // } neighbor_facet_dir = pos_it[second_node]; neighbor_facet_dir -= Vector(pos_it[first_node]); neighbor_facet_dir /= neighbor_facet_dir.norm(); Real dot = cent_facet_dir.dot(neighbor_facet_dir); if (dot < min_dot) { min_dot = dot; neighbor = elem; } } doubled_facets.add(neighbor); } } /* ----------------------------------------------------------------------- */ void ASRTools::preventCohesiveInsertionInNeighbors() { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto el_type = *mesh.elementTypes(dim, gt, _ek_regular).begin(); auto facet_type = Mesh::getFacetType(el_type); auto & doubled_nodes = mesh.getNodeGroup("doubled_nodes"); CSR nodes_to_elements; MeshUtils::buildNode2Elements(mesh_facets, nodes_to_elements, dim - 1); for (auto node : doubled_nodes.getNodes()) { for (auto & elem : nodes_to_elements.getRow(node)) { if (elem.type != facet_type) continue; inserter.getCheckFacets(elem.type, gt)(elem.element) = false; } } } /* ----------------------------------------------------------------------- */ void ASRTools::insertOppositeFacets(bool only_double_facets) { auto & coh_model = dynamic_cast(model); auto & inserter = coh_model.getElementInserter(); auto & insertion = inserter.getInsertionFacetsByElement(); auto & mesh = model.getMesh(); auto & mesh_facets = inserter.getMeshFacets(); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto & crack_facets = mesh.createElementGroup("crack_facets"); const auto & pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); auto & doubled_facets = mesh_facets.getElementGroup("doubled_facets"); /// sort facet numbers in doubled_facets_group to comply with new_elements /// event passed by the inserter doubled_facets.optimize(); /// instruct the inserter which facets to duplicate for (auto & type : doubled_facets.elementTypes(dim - 1)) { const auto element_ids = doubled_facets.getElements(type, gt); /// iterate over facets to duplicate for (auto && el : element_ids) { inserter.getCheckFacets(type, gt)(el) = false; Element new_facet{type, el, gt}; insertion(new_facet) = true; } } /// duplicate facets (and insert coh els if needed) inserter.insertElements(only_double_facets); - MeshUtils::resetFacetToDouble(mesh_facets); /// add facets connectivity to the mesh and the element group NewElementsEvent new_facets_event; for (auto & type : doubled_facets.elementTypes(dim - 1)) { const auto element_ids = doubled_facets.getElements(type, gt); if (not mesh.getConnectivities().exists(type, gt)) mesh.addConnectivityType(type, gt); auto & facet_conn = mesh.getConnectivity(type, gt); auto & mesh_facet_conn = mesh_facets.getConnectivity(type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto mesh_facet_conn_it = make_view(mesh_facet_conn, nb_nodes_facet).begin(); /// iterate over duplicated facets for (auto && el : element_ids) { Vector facet_nodes = mesh_facet_conn_it[el]; facet_conn.push_back(facet_nodes); Element new_facet{type, facet_conn.size() - 1, gt}; crack_facets.add(new_facet); for (auto & facet_node : facet_nodes) { crack_facets.addNode(facet_node, true); } new_facets_event.getList().push_back(new_facet); } } MeshUtils::fillElementToSubElementsData(mesh); mesh.sendEvent(new_facets_event); /// create an element group with nodes to apply Dirichlet model.getMesh().createElementGroupFromNodeGroup("doubled_nodes", "doubled_nodes", dim - 1); /// update FEEngineBoundary with new elements model.getFEEngineBoundary().initShapeFunctions(_not_ghost); model.getFEEngineBoundary().initShapeFunctions(_ghost); model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_not_ghost); model.getFEEngineBoundary().computeNormalsOnIntegrationPoints(_ghost); } /* ----------------------------------------------------------------------- */ void ASRTools::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; /// function is activated only when expanding cohesive elements is on if (not this->cohesive_insertion) return; if (this->doubled_nodes_ready) return; auto & mesh = model.getMesh(); auto & node_group = mesh.getNodeGroup("doubled_nodes"); auto & central_nodes = node_group.getNodes(); auto pos_it = make_view(mesh.getNodes(), mesh.getSpatialDimension()).begin(); for (auto & new_node : new_nodes) { const Vector & new_node_coord = pos_it[new_node]; for (auto & central_node : central_nodes) { const Vector & central_node_coord = pos_it[central_node]; if (new_node_coord == central_node_coord) node_group.add(new_node); } } this->doubled_nodes_ready = true; AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ void ASRTools::updateElementGroup(const std::string group_name) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); AKANTU_DEBUG_ASSERT(mesh.elementGroupExists(group_name), "Element group is not registered in the mesh"); auto dim = mesh.getSpatialDimension(); const GhostType gt = _not_ghost; auto && group = mesh.getElementGroup(group_name); auto && pos = mesh.getNodes(); const auto pos_it = make_view(pos, dim).begin(); for (auto & type : group.elementTypes(dim - 1)) { auto & facet_conn = mesh.getConnectivity(type, gt); const UInt nb_nodes_facet = facet_conn.getNbComponent(); const auto facet_nodes_it = make_view(facet_conn, nb_nodes_facet).begin(); AKANTU_DEBUG_ASSERT( type == _segment_2, "Currently update group works only for el type _segment_2"); const auto element_ids = group.getElements(type, gt); for (auto && el_id : element_ids) { const auto connected_els = mesh.getElementToSubelement(type, gt)(el_id); for (auto && connected_el : connected_els) { auto type_solid = connected_el.type; auto & solid_conn = mesh.getConnectivity(type_solid, gt); const UInt nb_nodes_solid_el = solid_conn.getNbComponent(); const auto solid_nodes_it = make_view(solid_conn, nb_nodes_solid_el).begin(); Vector facet_nodes = facet_nodes_it[el_id]; Vector solid_nodes = solid_nodes_it[connected_el.element]; // /// to which of connected elements facet belongs // Array solid_nodes(nb_nodes_solid_el); // for (UInt node : arange(nb_nodes_solid_el)) { // solid_nodes(node) = solid_nodes_it[connected_el.element](node); // } // /// check for the central facet node (id doesn't change nb) // auto id = solid_nodes.find(facet_nodes(2)); // if (id == UInt(-1)) // continue; /// check only the corner nodes of facets - central will not change for (auto f : arange(2)) { auto facet_node = facet_nodes(f); const Vector & facet_node_coords = pos_it[facet_node]; for (auto s : arange(nb_nodes_solid_el)) { auto solid_node = solid_nodes(s); const Vector & solid_node_coords = pos_it[solid_node]; if (solid_node_coords == facet_node_coords) { if (solid_node != facet_node) { // group.removeNode(facet_node); facet_conn(el_id, f) = solid_node; // group.addNode(solid_node, true); } break; } } } } } } group.optimize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------*/ // void ASRTools::applyDeltaU(Real delta_u) { // // get element group with nodes to apply Dirichlet // auto & crack_facets = model.getMesh().getElementGroup("doubled_nodes"); // DeltaU delta_u_bc(model, delta_u, getNodePairs()); // model.applyBC(delta_u_bc, crack_facets); // } /* -------------------------------------------------------------------------- */ void ASRTools::applyGelStrain(const Matrix & prestrain) { AKANTU_DEBUG_IN(); auto & mesh = model.getMesh(); auto dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(dim == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(dim, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::clearASREigenStrain() { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); Matrix zero_eigenstrain(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigenstrain; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::storeASREigenStrain(Array & stored_eig) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; UInt nb_quads = 0; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { const UInt nb_gp = model.getFEEngine().getNbIntegrationPoints(element_type); nb_quads += model.getMaterial("gel").getElementFilter(element_type).size() * nb_gp; } stored_eig.resize(nb_quads); auto stored_eig_it = stored_eig.begin(dim, dim); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it, ++stored_eig_it) (*stored_eig_it) = (*prestrain_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ASRTools::restoreASREigenStrain(Array & stored_eig) { AKANTU_DEBUG_IN(); const auto & mesh = model.getMesh(); const auto dim = mesh.getSpatialDimension(); GhostType gt = _not_ghost; auto stored_eig_it = stored_eig.begin(dim, dim); for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(model.getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(dim, dim); auto prestrain_end = prestrain_vect.end(dim, dim); for (; prestrain_it != prestrain_end; ++prestrain_it, ++stored_eig_it) (*prestrain_it) = (*stored_eig_it); } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake index 5c41df20d..1fb88938f 100644 --- a/packages/cohesive_element.cmake +++ b/packages/cohesive_element.cmake @@ -1,137 +1,139 @@ #=============================================================================== # @file cohesive_element.cmake # # @author Mauro Corrado # @author Nicolas Richart # @author Marco Vocialta # # @date creation: Tue Oct 16 2012 # @date last modification: Tue Jan 12 2016 # # @brief package description for cohesive elements # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # #=============================================================================== package_declare(cohesive_element DESCRIPTION "Use cohesive_element package of Akantu" DEPENDS lapack solid_mechanics) package_declare_sources(cohesive_element fe_engine/cohesive_element.hh fe_engine/fe_engine_template_cohesive.cc fe_engine/shape_cohesive.hh fe_engine/shape_cohesive_inline_impl.cc mesh_utils/cohesive_element_inserter.cc mesh_utils/cohesive_element_inserter.hh mesh_utils/cohesive_element_inserter_inline_impl.cc mesh_utils/cohesive_element_inserter_parallel.cc + mesh_utils/cohesive_element_inserter_helper.cc + mesh_utils/cohesive_element_inserter_helper.hh model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear_friction.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear_friction.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_includes.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_inline_impl.cc model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc ) package_declare_elements(cohesive_element ELEMENT_TYPES _cohesive_1d_2 _cohesive_2d_4 _cohesive_2d_6 _cohesive_3d_12 _cohesive_3d_16 _cohesive_3d_6 _cohesive_3d_8 KIND cohesive GEOMETRICAL_TYPES _gt_cohesive_1d_2 _gt_cohesive_2d_4 _gt_cohesive_2d_6 _gt_cohesive_3d_12 _gt_cohesive_3d_16 _gt_cohesive_3d_6 _gt_cohesive_3d_8 FE_ENGINE_LISTS compute_normals_on_integration_points contains get_shapes_derivatives gradient_on_integration_points interpolate_on_integration_points inverse_map lagrange_base ) package_declare_material_infos(cohesive_element LIST AKANTU_COHESIVE_MATERIAL_LIST INCLUDE material_cohesive_includes.hh ) package_declare_documentation_files(cohesive_element manual-cohesive_elements.tex manual-cohesive_elements_insertion.tex manual-cohesive_laws.tex manual-appendix-materials-cohesive.tex figures/cohesive2d.pdf figures/cohesive_exponential.pdf figures/linear_cohesive_law.pdf figures/bilinear_cohesive_law.pdf ) package_declare_documentation(cohesive_element "This package activates the cohesive elements engine within Akantu." "It depends on:" "\\begin{itemize}" " \\item A fortran compiler." " \\item An implementation of BLAS/LAPACK." "\\end{itemize}" ) diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index 2bb99bf84..625878ad8 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1379 +1,1386 @@ /** * @file aka_array_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 * @date last modification: Tue Feb 20 2018 * * @brief Inline functions of the classes Array and ArrayBase * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" /* -------------------------------------------------------------------------- */ #include "aka_static_memory.hh" #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ARRAY_TMPL_HH__ #define __AKANTU_AKA_ARRAY_TMPL_HH__ namespace akantu { namespace debug { struct ArrayException : public Exception {}; } // namespace debug /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(UInt size, UInt nb_component, const ID & id) : ArrayBase(id) { allocate(size, nb_component); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(UInt size, UInt nb_component, const_reference value, const ID & id) : ArrayBase(id) { allocate(size, nb_component, value); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(const ArrayDataLayer & vect, const ID & id) : ArrayBase(vect, id) { this->data_storage = vect.data_storage; this->size_ = vect.size_; this->nb_component = vect.nb_component; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer( const std::vector & vect) { this->data_storage = vect; this->size_ = vect.size(); this->nb_component = 1; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template -ArrayDataLayer & ArrayDataLayer:: -operator=(const ArrayDataLayer & other) { +ArrayDataLayer & +ArrayDataLayer::operator=(const ArrayDataLayer & other) { if (this != &other) { this->data_storage = other.data_storage; this->nb_component = other.nb_component; this->size_ = other.size_; this->values = this->data_storage.data(); } return *this; } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(ArrayDataLayer && other) = default; /* -------------------------------------------------------------------------- */ template -ArrayDataLayer & ArrayDataLayer:: -operator=(ArrayDataLayer && other) = default; +ArrayDataLayer & +ArrayDataLayer::operator=(ArrayDataLayer && other) = + default; /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(UInt new_size, UInt nb_component) { this->nb_component = nb_component; this->resize(new_size); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(UInt new_size, UInt nb_component, const T & val) { this->nb_component = nb_component; this->resize(new_size, val); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(UInt new_size) { this->data_storage.resize(new_size * this->nb_component); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(UInt new_size, const T & value) { this->data_storage.resize(new_size * this->nb_component, value); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::reserve(UInt size, UInt new_size) { if (new_size != UInt(-1)) { this->data_storage.resize(new_size * this->nb_component); } this->data_storage.reserve(size * this->nb_component); this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ /** * append a tuple to the array with the value value for all components * @param value the new last tuple or the array will contain nb_component copies * of value */ template inline void ArrayDataLayer::push_back(const T & value) { this->data_storage.push_back(value); this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ /** * append a matrix or a vector to the array * @param new_elem a reference to a Matrix or Vector */ template template