diff --git a/extra_packages/asr-tools/src/material_FE2/material_FE2.cc b/extra_packages/asr-tools/src/material_FE2/material_FE2.cc index 5bd1f5ab2..3106da55e 100644 --- a/extra_packages/asr-tools/src/material_FE2/material_FE2.cc +++ b/extra_packages/asr-tools/src/material_FE2/material_FE2.cc @@ -1,579 +1,578 @@ /** * @file material_FE2.cc * @author Emil Gallyamov * @author Aurelia Isabel Cuba Ramos * @date Tue Feb 8 2022 * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * @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 "material_FE2.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "solid_mechanics_model_RVE.hh" #include /* ------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialFE2::MaterialFE2(SolidMechanicsModel & model, const ID & id) : Parent(model, id), C("material_stiffness", *this), gelstrain("gelstrain", *this), /*non_reacted_gel("non_reacted_gel", *this),*/ crack_volume_ratio("crack_volume_ratio", *this), crack_volume_ratio_paste("crack_volume_ratio_paste", *this), crack_volume_ratio_agg("crack_volume_ratio_agg", *this) { AKANTU_DEBUG_IN(); this->C.initialize(voigt_h::size * voigt_h::size); this->gelstrain.initialize(spatial_dimension * spatial_dimension); this->crack_volume_ratio.initialize(1); this->crack_volume_ratio_paste.initialize(1); this->crack_volume_ratio_agg.initialize(1); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialFE2::~MaterialFE2() = default; /* -------------------------------------------------------------------------- */ template void MaterialFE2::initialize() { this->registerParam("element_type", el_type, _triangle_3, _pat_parsmod, "element type in RVE mesh"); this->registerParam("mesh_file", mesh_file, _pat_parsmod, "the mesh file for the RVE"); this->registerParam("nb_gel_pockets", nb_gel_pockets, _pat_parsmod, "the number of gel pockets in each RVE"); this->registerParam("eps_inf", eps_inf, _pat_parsmod, "asymptotic value of ASR expansion"); this->registerParam("U_C", U_C, _pat_parsmod, "thermal activation energy C"); this->registerParam("U_L", U_L, _pat_parsmod, "thermal activation energy L"); this->registerParam("T_ref", T_ref, _pat_parsmod, "reference temperature"); this->registerParam("time_lat_ref", time_lat_ref, _pat_parsmod, "latency time at the reference temperature"); this->registerParam("time_ch_ref", time_ch_ref, _pat_parsmod, "characteristic time at the reference temperature"); this->registerParam("reset_damage", reset_damage, false, _pat_parsmod, "reset damage"); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); // create a Mesh and SolidMechanicsModel on each integration point of the // material auto & mesh = this->model.getMesh(); auto & g_ids = mesh.template getData("global_ids", this->el_type); - // AKANTU_DEBUG_ASSERT(g_ids.size() > 0, "Global numbering array is empty"); auto const & element_filter = this->getElementFilter()(this->el_type); for (auto && data : zip(element_filter, make_view(C(this->el_type), voigt_h::size, voigt_h::size))) { UInt proc_el_id = std::get<0>(data); UInt gl_el_id = g_ids(proc_el_id); auto & C = std::get<1>(data); meshes.emplace_back(std::make_unique( spatial_dimension, "RVE_mesh_" + std::to_string(gl_el_id))); auto & mesh = *meshes.back(); mesh.read(mesh_file); RVEs.emplace_back(std::make_unique( mesh, true, this->nb_gel_pockets, _all_dimensions, "SMM_RVE_" + std::to_string(gl_el_id))); auto & RVE = *RVEs.back(); RVE.initFull(_analysis_method = _static); // compute intial stiffness of the RVE RVE.homogenizeStiffness(C, false); } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------*/ template void MaterialFE2::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Compute thermal stresses first Parent::computeStress(el_type, ghost_type); for (auto && data : zip(RVEs, make_view(this->gradu(el_type), spatial_dimension, spatial_dimension), make_view(this->stress(el_type), spatial_dimension, spatial_dimension), this->sigma_th(el_type), make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension), this->delta_T(this->el_type))) { auto & RVE = *(std::get<0>(data)); // reset nodal and internal fields RVE.resetNodalFields(); RVE.resetInternalFields(); // reset the damage to the previously converged state if (reset_damage) RVE.restoreDamageField(); // apply boundary conditions based on the current macroscopic displ. // gradient RVE.applyBoundaryConditionsRve(std::get<1>(data)); // advance the ASR in every RVE based on the new gel strain RVE.advanceASR(std::get<5>(data)); // compute the average average rVE stress RVE.homogenizeStressField(std::get<2>(data)); // compute the new effective stiffness of the RVE // decide whether stiffness homogenization is done via tension RVE.setStiffHomogenDir(std::get<2>(data)); // compute the new effective stiffness of the RVE RVE.homogenizeStiffness(std::get<4>(data), RVE.isTensileHomogen()); } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------*/ template void MaterialFE2::increaseGelStrain(Real & dt_day) { for (auto && data : zip(RVEs, this->delta_T(this->el_type), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension))) { auto & RVE = *(std::get<0>(data)); auto & strain_matrix = std::get<2>(data); Real ASRStrain = strain_matrix(0, 0); RVE.computeASRStrainLarive( dt_day, std::get<1>(data), ASRStrain, this->eps_inf, this->time_ch_ref, this->time_lat_ref, this->U_C, this->U_L, this->T_ref); for (UInt i = 0; i != spatial_dimension; ++i) strain_matrix(i, i) = ASRStrain; } } /* ------------------------------------------------------------------*/ template void MaterialFE2::increaseGelStrainArrhenius( Real & dt_day, const Real & k, const Real & Ea) { for (auto && data : zip(RVEs, this->delta_T(this->el_type), make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension))) { auto & RVE = *(std::get<0>(data)); auto & strain_matrix = std::get<2>(data); Real ASRStrain = strain_matrix(0, 0); RVE.computeASRStrainArrhenius(dt_day, std::get<1>(data), ASRStrain, k, Ea); for (UInt i = 0; i != spatial_dimension; ++i) strain_matrix(i, i) = ASRStrain; } } /* ------------------------------------------------------------------*/ template void MaterialFE2::updateStiffness() { AKANTU_DEBUG_IN(); for (const auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { Array & elem_filter = this->element_filter(type, _not_ghost); if (elem_filter.size() == 0) return; } for (auto && data : zip(RVEs, make_view(this->C(this->el_type), voigt_h::size, voigt_h::size), make_view(this->stress(el_type), spatial_dimension, spatial_dimension), this->delta_T(this->el_type), this->crack_volume_ratio(this->el_type), this->crack_volume_ratio_paste(this->el_type), this->crack_volume_ratio_agg(this->el_type))) { auto & RVE = *(std::get<0>(data)); if (reset_damage) RVE.storeDamageField(); // decide whether stiffness homogenization is done via tension RVE.setStiffHomogenDir(std::get<2>(data)); // compute the new effective stiffness of the RVE RVE.homogenizeStiffness(std::get<1>(data), RVE.isTensileHomogen()); // compute crack volume ratio in each RVE RVE.computeCrackVolume(std::get<4>(data)); // compute crack volume ratio per material RVE.computeCrackVolumePerMaterial(std::get<5>(data), "paste"); RVE.computeCrackVolumePerMaterial(std::get<6>(data), "aggregate"); } AKANTU_DEBUG_OUT(); } /* ---------------------------------------------------------------------- */ template void MaterialFE2::computeTangentModuli( ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::const_matrix_iterator C_it = this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); tangent.copy(*C_it); ++C_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* ----------------------------------------------------------- */ template Real MaterialFE2::computeAverageGelStrain() { AKANTU_DEBUG_IN(); Real av_gelstrain = 0; UInt nb_RVEs = 0; for (auto && data : enumerate(make_view(this->gelstrain(this->el_type), spatial_dimension, spatial_dimension))) { av_gelstrain += std::get<1>(data)(0, 0); nb_RVEs = std::get<0>(data) + 1; } auto && comm = akantu::Communicator::getWorldCommunicator(); comm.allReduce(av_gelstrain, SynchronizerOperation::_sum); comm.allReduce(nb_RVEs, SynchronizerOperation::_sum); return av_gelstrain /= nb_RVEs; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialFE2::setDirectoryToRveDumper( const std::string & directory) { AKANTU_DEBUG_IN(); for (auto && data : RVEs) { // set default directory to all RVEs data->setDirectory(directory); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialFE2::getNbRVEs() { AKANTU_DEBUG_IN(); UInt nb_RVEs = 0; for (auto && data : enumerate(RVEs)) { nb_RVEs = std::get<0>(data) + 1; } return nb_RVEs; AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------ */ template void MaterialFE2::dump() { AKANTU_DEBUG_IN(); for (auto && RVE : RVEs) { // update stress field before dumping RVE->assembleInternalForces(); // dump all the RVEs RVE->dumpRve(); } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------ */ template void MaterialFE2::dump(UInt dump_nb) { AKANTU_DEBUG_IN(); for (auto && RVE : RVEs) { // update stress field before dumping RVE->assembleInternalForces(); // dump all the RVEs RVE->dumpRve(dump_nb); } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------*/ template void MaterialFE2::setTimeStep(Real time_step) { AKANTU_DEBUG_IN(); for (auto && RVE : RVEs) { // set time step to all the RVEs RVE->setTimeStep(time_step); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------*/ template void MaterialFE2::saveRVEsState(std::string & restart_dir) { AKANTU_DEBUG_IN(); auto & mesh = this->model.getMesh(); auto & g_ids = mesh.template getData("global_ids", this->el_type); GhostType gt = _not_ghost; for (auto && data : enumerate(RVEs)) { auto proc_RVE_id = std::get<0>(data); UInt gl_RVE_id = g_ids(proc_RVE_id); auto & RVE = *(std::get<1>(data)); std::stringstream file_stream; file_stream << restart_dir << "/RVE_" << std::to_string(gl_RVE_id) << "_state.txt"; std::string RVE_state_file = file_stream.str(); std::ofstream file_output; file_output.open(RVE_state_file.c_str()); if (!file_output.is_open()) AKANTU_EXCEPTION("Could not create the file " + RVE_state_file + ", does its folder exist?"); UInt nb_materials = RVE.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = RVE.getMaterial(m); for (auto pair : mat.getInternalVectorsReal()) { // set up precision constexpr auto max_digits10 = std::numeric_limits::max_digits10; file_output << std::setprecision(max_digits10); // get the internal field const InternalField & internal_real = *pair.second; // loop over all types in that material for (auto element_type : internal_real.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; auto && internal_field_array = internal_real(element_type, gt); for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { file_output << internal_field_array(i, j) << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + RVE_state_file); } } } } for (auto pair : mat.getInternalVectorsUInt()) { // set up precision constexpr auto max_digits10 = std::numeric_limits::max_digits10; file_output << std::setprecision(max_digits10); // get the internal field const InternalField & internal_uint = *pair.second; // loop over all types in that material for (auto element_type : internal_uint.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; const Array & internal_field_array = internal_uint(element_type, gt); for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { file_output << internal_field_array(i, j) << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + RVE_state_file); } } } } for (auto pair : mat.getInternalVectorsBool()) { // set up precision constexpr auto max_digits10 = std::numeric_limits::max_digits10; file_output << std::setprecision(max_digits10); // get the internal field const InternalField & internal_bool = *pair.second; // loop over all types in that material for (auto element_type : internal_bool.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; const Array & internal_field_array = internal_bool(element_type, gt); for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { file_output << internal_field_array(i, j) << std::endl; if (file_output.fail()) AKANTU_EXCEPTION("Error in writing data to file " + RVE_state_file); } } } } } file_output.close(); AKANTU_DEBUG_OUT(); } } /* ------------------------------------------------------------------ */ template void MaterialFE2::loadRVEsState(std::string & restart_dir) { // variables for parallel execution auto & mesh = this->model.getMesh(); auto & g_ids = mesh.template getData("global_ids", this->el_type); GhostType gt = _not_ghost; for (auto && data : enumerate(RVEs)) { auto proc_RVE_id = std::get<0>(data); UInt gl_RVE_id = g_ids(proc_RVE_id); auto & RVE = *(std::get<1>(data)); // open the input file std::stringstream file_stream; file_stream << restart_dir << "/RVE_" << std::to_string(gl_RVE_id) << "_state.txt"; std::string RVE_state_file = file_stream.str(); std::ifstream file_input; file_input.open(RVE_state_file.c_str()); if (!file_input.is_open()) AKANTU_EXCEPTION("Could not open file " + RVE_state_file); auto nb_materials = RVE.getNbMaterials(); for (UInt m = 0; m < nb_materials; ++m) { const Material & mat = RVE.getMaterial(m); // get the internal field that need to be set for (auto pair : mat.getInternalVectorsReal()) { // get the internal field InternalField & internal_real = *pair.second; // loop over all types in that material for (auto element_type : internal_real.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; Array & internal_field_array = internal_real(element_type, gt); std::string line; for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { std::getline(file_input, line); if (file_input.fail()) AKANTU_EXCEPTION("Could not read data from file " + RVE_state_file); std::stringstream sstr(line); sstr >> internal_field_array(i, j); } } } } // get the internal field that need to be set for (auto pair : mat.getInternalVectorsUInt()) { // get the internal field InternalField & internal_uint = *pair.second; // loop over all types in that material for (auto element_type : internal_uint.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; Array & internal_field_array = internal_uint(element_type, gt); std::string line; for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { std::getline(file_input, line); if (file_input.fail()) AKANTU_EXCEPTION("Could not read data from file " + RVE_state_file); std::stringstream sstr(line); sstr >> internal_field_array(i, j); } } } } // get the internal field that need to be set for (auto pair : mat.getInternalVectorsBool()) { // get the internal field InternalField & internal_bool = *pair.second; // loop over all types in that material for (auto element_type : internal_bool.elementTypes(gt)) { const Array & elem_filter = mat.getElementFilter(element_type, gt); if (!elem_filter.size()) continue; Array & internal_field_array = internal_bool(element_type, gt); std::string line; for (UInt i = 0; i < internal_field_array.size(); ++i) { for (UInt j = 0; j < internal_field_array.getNbComponent(); ++j) { std::getline(file_input, line); if (file_input.fail()) AKANTU_EXCEPTION("Could not read data from file " + RVE_state_file); std::stringstream sstr(line); sstr >> internal_field_array(i, j); } } } } } // close the file file_input.close(); } } INSTANTIATE_MATERIAL(material_FE2, MaterialFE2); } // namespace akantu diff --git a/packages/fluid_diffusion.cmake b/packages/fluid_diffusion.cmake deleted file mode 100644 index b64028200..000000000 --- a/packages/fluid_diffusion.cmake +++ /dev/null @@ -1,41 +0,0 @@ -#=============================================================================== -# @file heat_transfer.cmake -# -# @author Guillaume Anciaux -# @author Nicolas Richart -# -# @date creation: Mon Nov 21 2011 -# @date last modification: Mon Mar 30 2015 -# -# @brief package description for heat transfer -# -# @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(fluid_diffusion - DESCRIPTION "Use Fluid Diffusion package of Akantu") - -package_declare_sources(fluid_diffusion - model/fluid_diffusion/fluid_diffusion_model.cc - model/fluid_diffusion/fluid_diffusion_model.hh - model/fluid_diffusion/fluid_diffusion_model_inline_impl.hh - model/fluid_diffusion/fluid_diffusion_model_event_handler.hh - ) diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 653e2224e..ae263b50f 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,691 +1,683 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_COMMON_HH_ #define AKANTU_COMMON_HH_ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Constants */ /* -------------------------------------------------------------------------- */ namespace { [[gnu::unused]] constexpr UInt _all_dimensions{ std::numeric_limits::max()}; #ifdef AKANTU_NDEBUG [[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.}; #else [[gnu::unused]] constexpr Real REAL_INIT_VALUE{ std::numeric_limits::quiet_NaN()}; #endif } // namespace /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "aka_enum_macros.hh" /* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of /// events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ - (fluid_diffusion_model) \ (structural_mechanics_model) \ (embedded_model) \ (phase_field_model) \ (coupler_solid_phasefield) // clang-format on /// enum ModelType defines which type of physics is solved AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { model, solid_mechanics_model, solid_mechanics_model_cohesive, heat_transfer_model, structural_mechanics_model, embedded_model, }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_NON_LINEAR_SOLVER_TYPES \ (linear) \ (newton_raphson) \ (newton_raphson_modified) \ (lumped) \ (gmres) \ (bfgs) \ (cg) \ (auto) // clang-format on AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES) #else /// Type of non linear resolution available in akantu enum class NonLinearSolverType { _linear, ///< No non linear convergence loop _newton_raphson, ///< Regular Newton-Raphson _newton_raphson_modified, ///< Newton-Raphson with initial tangent _lumped, ///< Case of lumped mass or equivalent matrix _gmres, _bfgs, _cg, _auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_TIME_STEP_SOLVER_TYPE \ (static) \ (dynamic) \ (dynamic_lumped) \ (not_defined) // clang-format on AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE) #else /// Type of time stepping solver enum class TimeStepSolverType { _static, ///< Static solution _dynamic, ///< Dynamic solver _dynamic_lumped, ///< Dynamic solver with lumped mass _not_defined, ///< For not defined cases }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_INTEGRATION_SCHEME_TYPE \ (pseudo_time) \ (forward_euler) \ (trapezoidal_rule_1) \ (backward_euler) \ (central_difference) \ (fox_goodwin) \ (trapezoidal_rule_2) \ (linear_acceleration) \ (newmark_beta) \ (generalized_trapezoidal) // clang-format on AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE) #else /// Type of integration scheme enum class IntegrationSchemeType { _pseudo_time, ///< Pseudo Time _forward_euler, ///< GeneralizedTrapezoidal(0) _trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _backward_euler, ///< GeneralizedTrapezoidal(1) _central_difference, ///< NewmarkBeta(0, 1/2) _fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; #endif #if !defined(DOXYGEN) // clang-format off #define AKANTU_SOLVE_CONVERGENCE_CRITERIA \ (residual) \ (solution) \ (residual_mass_wgh) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA) #else /// enum SolveConvergenceCriteria different convergence criteria enum class SolveConvergenceCriteria { _residual, ///< Use residual to test the convergence _solution, ///< Use solution to test the convergence _residual_mass_wgh ///< Use residual weighted by inv. nodal mass to ///< testb }; #endif /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum MatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; #if !defined(DOXYGEN) // clang-format off #define AKANTU_SYNCHRONIZATION_TAG \ (whatever) \ (update) \ (ask_nodes) \ (size) \ (smm_mass) \ (smm_for_gradu) \ (smm_boundary) \ (smm_uv) \ (smm_res) \ (smm_init_mat) \ (smm_stress) \ (smmc_facets) \ (smmc_facets_conn) \ (smmc_facets_stress) \ (smmc_damage) \ (giu_global_conn) \ (ce_groups) \ (ce_insertion_order) \ (gm_clusters) \ (htm_temperature) \ (htm_gradient_temperature) \ - (fdm_pressure) \ - (fdm_gradient_pressure) \ (htm_phi) \ (htm_gradient_phi) \ (pfm_damage) \ (pfm_driving) \ (pfm_history) \ (pfm_energy) \ (csp_damage) \ (csp_strain) \ (mnl_for_average) \ (mnl_weight) \ (nh_criterion) \ (test) \ (user_1) \ (user_2) \ (material_id) \ (for_dump) \ (cf_nodal) \ (cf_incr) \ (solver_solution) \ (crack_nb) \ (border_nodes) // clang-format on AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG) #else /// @enum SynchronizationTag type of synchronizations enum class SynchronizationTag { //--- Generic tags --- _whatever, _update, _ask_nodes, _size, //--- SolidMechanicsModel tags --- _smm_mass, ///< synchronization of the SolidMechanicsModel.mass _smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _smm_uv, ///< synchronization of the nodal velocities and displacement _smm_res, ///< synchronization of the nodal residual _smm_init_mat, ///< synchronization of the data to initialize materials _smm_stress, ///< synchronization of the stresses to compute the ///< internal /// forces _smmc_facets, ///< synchronization of facet data to setup facet synch _smmc_facets_conn, ///< synchronization of facet global connectivity _smmc_facets_stress, ///< synchronization of facets' stress to setup ///< facet /// synch _smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups _ce_insertion_order, ///< synchronization of the order of insertion of /// cohesive elements // --- GroupManager tags --- _gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _htm_temperature, ///< synchronization of the nodal temperature _htm_gradient_temperature, ///< synchronization of the element gradient /// temperature - // --- FluidDiffusion tags --- - _htm_pressure, ///< synchronization of the nodal pressure - _htm_gradient_pressure, ///< synchronization of the element gradient - /// pressure - // --- PhaseFieldModel tags --- _pfm_damage, ///< synchronization of the nodal damage _pfm_driving, ///< synchronization of the driving forces to /// compute the internal _pfm_history, ///< synchronization of the damage history to /// compute the internal _pfm_energy, ///< synchronization of the damage energy /// density to compute the internal // --- CouplerSolidPhaseField tags --- _csp_damage, ///< synchronization of the damage from phase /// model to solid model _csp_strain, ///< synchronization of the strain from solid /// model to phase model // --- LevelSet tags --- _htm_phi, ///< synchronization of the nodal level set value phi _htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _mnl_for_average, ///< synchronization of data to average in non local /// material _mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _nh_criterion, // --- General tags --- _test, ///< Test tag _user_1, ///< tag for user simulations _user_2, ///< tag for user simulations _material_id, ///< synchronization of the material ids _for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _cf_nodal, ///< synchronization of disp, velo, and current position _cf_incr, ///< synchronization of increment // --- Solver tags --- _solver_solution, ///< synchronization of the solution obained with the /// PETSc solver // --- ASRTools tags --- _crack_nb ///< synchronization for crack numbers in asr tools _border_nodes ///< synchronization for the nodes at the partition border }; #endif /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } std::ostream & operator<<(std::ostream & stream, NodeFlag flag); } // namespace akantu AKANTU_ENUM_HASH(GhostType) namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; namespace { constexpr ghost_type_t ghost_types{_casper}; } /// standard output stream operator for GhostType // inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT ' ' #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_DEREF_PTR(name, ptr) \ inline const auto & get##name() const { \ if (not(ptr)) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*(ptr)); \ } #define AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(name, ptr) \ inline auto & get##name() { \ if (not(ptr)) { \ AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \ } \ return (*(ptr)); \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name(const support & el_type, \ GhostType ghost_type = _not_ghost) \ con { /* NOLINT */ \ return variable(el_type, ghost_type); \ } // NOLINT #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ struct TensorTrait {}; struct TensorProxyTrait {}; } // namespace akantu /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ namespace aka { /* ------------------------------------------------------------------------ */ template using is_tensor = std::is_base_of; template using is_tensor_proxy = std::is_base_of; /* ------------------------------------------------------------------------ */ template using is_scalar = std::is_arithmetic; /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> bool is_of_type(T && t) { return ( dynamic_cast>::value, std::add_const_t, R>>>(&t) != nullptr); } /* -------------------------------------------------------------------------- */ template bool is_of_type(std::unique_ptr & t) { return ( dynamic_cast::value, std::add_const_t, R>>>( t.get()) != nullptr); } /* ------------------------------------------------------------------------ */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { static_assert( disjunction< std::is_base_of, std::decay_t>, // down-cast std::is_base_of, std::decay_t> // up-cast >::value, "Type T and R are not valid for a as_type conversion"); return dynamic_cast>::value, std::add_const_t, R>>>(t); } /* -------------------------------------------------------------------------- */ template ::value> * = nullptr> decltype(auto) as_type(T && t) { return &as_type(*t); } /* -------------------------------------------------------------------------- */ template decltype(auto) as_type(const std::shared_ptr & t) { return std::dynamic_pointer_cast(t); } } // namespace aka #include "aka_common_inline_impl.hh" #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); #define AKANTU_CURRENT_FUNCTION \ (std::string(__func__) + "():" + std::to_string(__LINE__)) } // namespace akantu /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic * number is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif // AKANTU_COMMON_HH_ diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in index 51e8731aa..32662ed34 100644 --- a/src/common/aka_config.hh.in +++ b/src/common/aka_config.hh.in @@ -1,101 +1,102 @@ /** * @file aka_config.hh.in * * @author Nicolas Richart * * @date creation: Sun Sep 26 2010 * @date last modification: Thu Jan 25 2018 * * @brief Compilation time configuration of Akantu * * * 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_AKA_CONFIG_HH_ #define AKANTU_AKA_CONFIG_HH_ +// clang-format off #define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@ #define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@ #define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@ #define AKANTU_VERSION \ (AKANTU_VERSION_MAJOR * 10000 + AKANTU_VERSION_MINOR * 100 + \ AKANTU_VERSION_PATCH) @AKANTU_TYPES_EXTRA_INCLUDES@ namespace akantu { using Real = @AKANTU_FLOAT_TYPE@; using Int = @AKANTU_SIGNED_INTEGER_TYPE@; using UInt = @AKANTU_UNSIGNED_INTEGER_TYPE@; } // akantu #define AKANTU_INTEGER_SIZE @AKANTU_INTEGER_SIZE@ #define AKANTU_FLOAT_SIZE @AKANTU_FLOAT_SIZE@ #cmakedefine AKANTU_HAS_STRDUP #cmakedefine AKANTU_USE_BLAS #cmakedefine AKANTU_USE_LAPACK #cmakedefine AKANTU_PARALLEL #cmakedefine AKANTU_USE_MPI #cmakedefine AKANTU_USE_SCOTCH #cmakedefine AKANTU_USE_PTSCOTCH #cmakedefine AKANTU_SCOTCH_NO_EXTERN #cmakedefine AKANTU_IMPLICIT #cmakedefine AKANTU_USE_MUMPS #cmakedefine AKANTU_USE_PETSC #cmakedefine AKANTU_USE_IOHELPER #cmakedefine AKANTU_USE_QVIEW #cmakedefine AKANTU_USE_BLACKDYNAMITE #cmakedefine AKANTU_USE_PYBIND11 #cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY #cmakedefine AKANTU_EXTRA_MATERIALS #cmakedefine AKANTU_STUDENTS_EXTRA_PACKAGE #cmakedefine AKANTU_DAMAGE_NON_LOCAL #cmakedefine AKANTU_SOLID_MECHANICS #cmakedefine AKANTU_STRUCTURAL_MECHANICS #cmakedefine AKANTU_HEAT_TRANSFER -#cmakedefine AKANTU_FLUID_DIFFUSION #cmakedefine AKANTU_PHASE_FIELD #cmakedefine AKANTU_PYTHON_INTERFACE #cmakedefine AKANTU_COHESIVE_ELEMENT #cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT #cmakedefine MODEL_COUPLERS #cmakedefine AKANTU_IGFEM #cmakedefine AKANTU_USE_CGAL #cmakedefine AKANTU_EMBEDDED // Debug tools //#cmakedefine AKANTU_NDEBUG #cmakedefine AKANTU_DEBUG_TOOLS #cmakedefine READLINK_COMMAND @READLINK_COMMAND@ #cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@ +// clang-format on #endif /* AKANTU_AKA_CONFIG_HH_ */ diff --git a/src/fe_engine/shape_lagrange.hh b/src/fe_engine/shape_lagrange.hh index eb8c44aad..d31f6d526 100644 --- a/src/fe_engine/shape_lagrange.hh +++ b/src/fe_engine/shape_lagrange.hh @@ -1,183 +1,172 @@ /** * @file shape_lagrange.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Feb 15 2011 * @date last modification: Mon Jan 29 2018 * * @brief lagrangian shape functions class * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "shape_lagrange_base.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SHAPE_LAGRANGE_HH_ #define AKANTU_SHAPE_LAGRANGE_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template class ShapeCohesive; class ShapeIGFEM; template class ShapeLagrange : public ShapeLagrangeBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeLagrange(const Mesh & mesh, UInt spatial_dimension, const ID & id = "shape_lagrange"); ~ShapeLagrange() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialization function for structural elements not yet implemented inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, - ElementType type, - GhostType ghost_type); + ElementType type, GhostType ghost_type); /// computes the shape functions derivatives for given interpolation points template void computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, GhostType ghost_type, const Array & filter_elements = empty_filter) const; void computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, - Array & shape_derivatives, ElementType type, - GhostType ghost_type, + Array & shape_derivatives, ElementType type, GhostType ghost_type, const Array & filter_elements) const override; - /// computes the shape functions variations for spatial dim != natural dim - template - void computeShapeDerivativesOnIntegrationPoints1DIn2D( - const Array & nodes, const Matrix & integration_points, - Array & shape_derivatives, const GhostType & ghost_type, - const Array & filter_elements = empty_filter) const; - /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element integration points from /// natural coordinates template - void - precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, - GhostType ghost_type); + void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, + GhostType ghost_type); /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; template void interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, const Array & shapes, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// interpolate on physical point template void interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, - Vector & interpolated, - GhostType ghost_type) const; + Vector & interpolated, GhostType ghost_type) const; /// compute the gradient of u on the integration points template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; template void computeBtD(const Array & Ds, Array & BtDs, GhostType ghost_type, const Array & filter_elements) const; template void computeBtDB(const Array & Ds, Array & BtDBs, UInt order_d, GhostType ghost_type, const Array & filter_elements) const; /// multiply a field by shape functions @f$ fts_{ij} = f_i * \varphi_j @f$ template void computeNtb(const Array & bs, Array & Ntbs, GhostType ghost_type, const Array & filter_elements = empty_filter) const; template void computeNtbN(const Array & bs, Array & NtbNs, GhostType ghost_type, const Array & filter_elements) const; /// find natural coords from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, GhostType ghost_type = _not_ghost) const; /// return true if the coordinates provided are inside the element, false /// otherwise template bool contains(const Vector & real_coords, UInt elem, GhostType ghost_type) const; /// compute the shape on a provided point template void computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, GhostType ghost_type) const; /// compute the shape derivatives on a provided point template void computeShapeDerivatives(const Matrix & real_coords, UInt elem, Tensor3 & shapes, GhostType ghost_type) const; protected: /// compute the shape derivatives on integration points for a given element template inline void computeShapeDerivativesOnCPointsByElement(const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "shape_lagrange_inline_impl.hh" #endif /* AKANTU_SHAPE_LAGRANGE_HH_ */ diff --git a/src/fe_engine/shape_lagrange_inline_impl.hh b/src/fe_engine/shape_lagrange_inline_impl.hh index 86d581981..b97a78761 100644 --- a/src/fe_engine/shape_lagrange_inline_impl.hh +++ b/src/fe_engine/shape_lagrange_inline_impl.hh @@ -1,654 +1,577 @@ /** * @file shape_lagrange_inline_impl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Oct 27 2010 * @date last modification: Tue Feb 20 2018 * * @brief ShapeLagrange inline implementation * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "aka_voigthelper.hh" #include "fe_engine.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_ #define AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ if (ElementClass::getNaturalSpaceDimension() == \ mesh.getSpatialDimension() || \ kind != _ek_regular) \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); -// TODO: at this stage precomputing derivatives should not be done for the fluid -// diffusion model -// id.find("fluid_diffusion_model") != std::string::npos) + template inline void ShapeLagrange::initShapeFunctions( const Array & nodes, const Matrix & integration_points, ElementType type, GhostType ghost_type) { AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } #undef INIT_SHAPE_FUNCTIONS /* -------------------------------------------------------------------------- */ template template inline void ShapeLagrange::computeShapeDerivativesOnCPointsByElement( const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const { AKANTU_DEBUG_IN(); // compute dnds Tensor3 dnds(node_coords.rows(), node_coords.cols(), natural_coords.cols()); ElementClass::computeDNDS(natural_coords, dnds); // compute jacobian Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute dndx ElementClass::computeShapeDerivatives(J, dnds, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ElementClass::inverseMap(real_coords, nodes_coord, natural_coords); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template bool ShapeLagrange::contains(const Vector & real_coords, UInt elem, GhostType ghost_type) const { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); return ElementClass::contains(natural_coords); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, Vector & interpolated, GhostType ghost_type) const { UInt nb_shapes = ElementClass::getShapeSize(); Vector shapes(nb_shapes); computeShapes(real_coords, elem, shapes, ghost_type); ElementClass::interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); ElementClass::computeShapes(natural_coords, shapes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivatives( const Matrix & real_coords, UInt elem, Tensor3 & shapesd, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_points = real_coords.cols(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == shapesd.size(0) && nb_nodes_per_element == shapesd.size(1), "Shape size doesn't match"); AKANTU_DEBUG_ASSERT(nb_points == shapesd.size(2), "Number of points doesn't match shapes size"); Matrix natural_coords(spatial_dimension, nb_points); // Creates the matrix of natural coordinates for (UInt i = 0; i < nb_points; i++) { Vector real_point = real_coords(i); Vector natural_point = natural_coords(i); inverseMap(real_point, elem, natural_point, ghost_type); } UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); computeShapeDerivativesOnCPointsByElement(nodes_coord, natural_coords, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template ShapeLagrange::ShapeLagrange(const Mesh & mesh, UInt spatial_dimension, const ID & id) : ShapeLagrangeBase(mesh, spatial_dimension, kind, id) {} /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd, "The shapes_derivatives array does not have the correct " << "number of component"); shape_derivatives.resize(nb_element * nb_points); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, filter_elements); Real * shapesd_val = shape_derivatives.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { if (filter_elements != empty_filter) { shapesd_val = shape_derivatives.storage() + filter_elements(elem) * size_of_shapesd * nb_points; } Matrix & X = *x_it; Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, nb_points); computeShapeDerivativesOnCPointsByElement(X, integration_points, B); if (filter_elements == empty_filter) { shapesd_val += size_of_shapesd * nb_points; } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -template -template -void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints1DIn2D( - const Array & nodes, const Matrix & integration_points, - Array & shape_derivatives, const GhostType & ghost_type, - const Array & filter_elements) const { - AKANTU_DEBUG_IN(); - - UInt spatial_dimension = mesh.getSpatialDimension(); - UInt element_dimension = ElementClass::getNaturalSpaceDimension(); - - UInt nb_nodes_per_element = - ElementClass::getNbNodesPerInterpolationElement(); - - UInt nb_points = integration_points.cols(); - UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); - - UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); - AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd, - "The shapes_derivatives array does not have the correct " - << "number of component"); - shape_derivatives.resize(nb_element * nb_points); - - Array x_el(0, spatial_dimension * nb_nodes_per_element); - FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, - filter_elements); - - Real * shapesd_val = shape_derivatives.storage(); - Array::matrix_iterator x_it = - x_el.begin(spatial_dimension, nb_nodes_per_element); - - /// array of equivalent coordinates (works for 1D element in 2D only) - Array x_el_equiv(nb_element, element_dimension * nb_nodes_per_element); - for (auto && data : - zip(make_view(x_el, spatial_dimension, nb_nodes_per_element), - make_view(x_el_equiv, element_dimension, nb_nodes_per_element))) { - const Matrix & x = std::get<0>(data); - auto & x_eq = std::get<1>(data); - for (auto node_nb : arange(1, nb_nodes_per_element)) { - Vector delta_x = Vector(x(node_nb)) - Vector(x(0)); - x_eq(element_dimension - 1, node_nb) = delta_x.norm(); - } - } - - if (filter_elements != empty_filter) - nb_element = filter_elements.size(); - - for (auto && data : enumerate( - make_view(x_el_equiv, element_dimension, nb_nodes_per_element))) { - auto & elem = std::get<0>(data); - const auto & X = std::get<1>(data); - if (filter_elements != empty_filter) - shapesd_val = shape_derivatives.storage() + - filter_elements(elem) * size_of_shapesd * nb_points; - - Tensor3 B(shapesd_val, element_dimension, nb_nodes_per_element, - nb_points); - computeShapeDerivativesOnCPointsByElement(X, integration_points, B); - - if (filter_elements == empty_filter) - shapesd_val += size_of_shapesd * nb_points; - } - - AKANTU_DEBUG_OUT(); -} - /* -------------------------------------------------------------------------- */ template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, ElementType type, GhostType ghost_type, const Array & filter_elements) const { #define AKANTU_COMPUTE_SHAPES(type) \ computeShapeDerivativesOnIntegrationPoints( \ nodes, integration_points, shape_derivatives, ghost_type, \ filter_elements); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES); #undef AKANTU_COMPUTE_SHAPES } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapesOnIntegrationPoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapes = ElementClass::getShapeSize(); Array & shapes_tmp = shapes.alloc(0, size_of_shapes, itp_type, ghost_type); this->computeShapesOnIntegrationPoints(nodes, natural_coords, shapes_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); - UInt spatial_dimension = mesh.getSpatialDimension(); - UInt natural_dimension = ElementClass::getNaturalSpaceDimension(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); - if (spatial_dimension != natural_dimension) { - this->computeShapeDerivativesOnIntegrationPoints1DIn2D( - nodes, natural_coords, shapes_derivatives_tmp, ghost_type); - } else { - this->computeShapeDerivativesOnIntegrationPoints( - nodes, natural_coords, shapes_derivatives_tmp, ghost_type); - } + this->computeShapeDerivativesOnIntegrationPoints( + nodes, natural_coords, shapes_derivatives_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, const Array & shapes, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type), "No shapes for the type " << shapes.printType(itp_type, ghost_type)); this->interpolateOnIntegrationPoints(in_u, out_uq, nb_degree_of_freedom, shapes(itp_type, ghost_type), ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->gradientElementalFieldOnIntegrationPoints( u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtD( const Array & Ds, Array & BtDs, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); - // auto spatial_dimension = mesh.getSpatialDimension(); - auto spatial_dimension = ElementClass::getNaturalSpaceDimension(); + auto spatial_dimension = mesh.getSpatialDimension(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, spatial_dimension, nb_nodes_per_element); auto B_it = view.begin(); auto B_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes_derivatives, shapes_derivatives_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_derivatives_filtered, spatial_dimension, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } for (auto && values : zip(range(B_it, B_end), make_view(Ds, Ds.getNbComponent() / spatial_dimension, spatial_dimension), make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D = std::get<2>(values); // transposed due to the storage layout of B Bt_D.template mul(D, B); } } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtDB( const Array & Ds, Array & BtDBs, UInt order_d, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); constexpr auto dim = ElementClass::getSpatialDimension(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, dim, nb_nodes_per_element); auto B_it = view.begin(); auto B_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes_derivatives, shapes_derivatives_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_derivatives_filtered, dim, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } if (order_d == 4) { UInt tangent_size = VoigtHelper::size; Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); for (auto && values : zip(range(B_it, B_end), make_view(Ds, tangent_size, tangent_size), make_view(BtDBs, dim * nb_nodes_per_element, dim * nb_nodes_per_element))) { const auto & Bfull = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); VoigtHelper::transferBMatrixToSymVoigtBMatrix(Bfull, B, nb_nodes_per_element); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } else if (order_d == 2) { Matrix Bt_D(nb_nodes_per_element, dim); for (auto && values : zip(range(B_it, B_end), make_view(Ds, dim, dim), make_view(BtDBs, nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } } template <> template <> inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>( const Array & /*Ds*/, Array & /*BtDBs*/, UInt /*order_d*/, GhostType /*ghost_type*/, const Array & /*filter_elements*/) const { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeNtbN( const Array & bs, Array & NtbNs, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; auto size_of_shapes = ElementClass::getShapeSize(); auto nb_degree_of_freedom = bs.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_filtered(0, size_of_shapes); auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes); auto N_it = view.begin(); auto N_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type), shapes_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_filtered, 1, size_of_shapes); N_it = view.begin(); N_end = view.end(); } Matrix Nt_b(nb_nodes_per_element, nb_degree_of_freedom); for (auto && values : zip(range(N_it, N_end), make_view(bs, nb_degree_of_freedom, 1), make_view(NtbNs, nb_nodes_per_element, nb_nodes_per_element))) { const auto & N = std::get<0>(values); const auto & b = std::get<1>(values); auto & Nt_b_N = std::get<2>(values); Nt_b.template mul(N, b); Nt_b_N.template mul(Nt_b, N); } } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeNtb( const Array & bs, Array & Ntbs, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); Ntbs.resize(bs.size()); auto size_of_shapes = ElementClass::getShapeSize(); auto itp_type = ElementClassProperty::interpolation_type; auto nb_degree_of_freedom = bs.getNbComponent(); Array shapes_filtered(0, size_of_shapes); auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes); auto N_it = view.begin(); auto N_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type), shapes_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_filtered, 1, size_of_shapes); N_it = view.begin(); N_end = view.end(); } for (auto && values : zip(make_view(bs, nb_degree_of_freedom, 1), range(N_it, N_end), make_view(Ntbs, nb_degree_of_freedom, size_of_shapes))) { const auto & b = std::get<0>(values); const auto & N = std::get<1>(values); auto & Ntb = std::get<2>(values); Ntb.template mul(b, N); } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_ */ diff --git a/src/mesh/element_type_map.cc b/src/mesh/element_type_map.cc index fbf44aaa0..34bc20875 100644 --- a/src/mesh/element_type_map.cc +++ b/src/mesh/element_type_map.cc @@ -1,73 +1,71 @@ /** * @file element_type_map.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief A Documented file. * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, UInt nb_component, UInt spatial_dimension, GhostType ghost_type, ElementKind element_kind) : MeshElementTypeMapArrayInitializer( fe_engine.getMesh(), nb_component, spatial_dimension == UInt(-2) - ? fe_engine.getElementDimension() + ? fe_engine.getMesh().getSpatialDimension() : spatial_dimension, ghost_type, element_kind, true, false), fe_engine(fe_engine) {} FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, const ElementTypeMapArrayInitializer::CompFunc & nb_component, - UInt spatial_dimension, GhostType ghost_type, - ElementKind element_kind) + UInt spatial_dimension, GhostType ghost_type, ElementKind element_kind) : MeshElementTypeMapArrayInitializer( fe_engine.getMesh(), nb_component, spatial_dimension == UInt(-2) - ? fe_engine.getElementDimension() + ? fe_engine.getMesh().getSpatialDimension() : spatial_dimension, ghost_type, element_kind, true, false), fe_engine(fe_engine) {} -UInt FEEngineElementTypeMapArrayInitializer::size( - ElementType type) const { +UInt FEEngineElementTypeMapArrayInitializer::size(ElementType type) const { return MeshElementTypeMapArrayInitializer::size(type) * fe_engine.getNbIntegrationPoints(type, this->ghost_type); } FEEngineElementTypeMapArrayInitializer::ElementTypesIteratorHelper FEEngineElementTypeMapArrayInitializer::elementTypes() const { return this->fe_engine.elementTypes(spatial_dimension, ghost_type, element_kind); } } // namespace akantu diff --git a/src/model/common/integration_scheme/integration_scheme.cc b/src/model/common/integration_scheme/integration_scheme.cc index 5a68bf6aa..ae1fb6508 100644 --- a/src/model/common/integration_scheme/integration_scheme.cc +++ b/src/model/common/integration_scheme/integration_scheme.cc @@ -1,96 +1,92 @@ /** * @file integration_scheme.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Jan 31 2018 * * @brief Common interface to all interface schemes * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "integration_scheme.hh" #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ IntegrationScheme::IntegrationScheme(DOFManager & dof_manager, const ID & dof_id, UInt order) : Parsable(ParserType::_integration_scheme, dof_id), dof_manager(dof_manager), dof_id(dof_id), order(order), u_store(order + 1) {} /* -------------------------------------------------------------------------- */ /// standard input stream operator for SolutionType std::istream & operator>>(std::istream & stream, IntegrationScheme::SolutionType & type) { std::string str; stream >> str; if (str == "displacement") { type = IntegrationScheme::_displacement; } else if (str == "temperature") { type = IntegrationScheme::_temperature; } else if (str == "velocity") { type = IntegrationScheme::_velocity; } else if (str == "temperature_rate") { type = IntegrationScheme::_temperature_rate; - } else if (str == "pressure") { - type = IntegrationScheme::_pressure; - } else if (str == "pressure_rate") { - type = IntegrationScheme::_pressure_rate; } else if (str == "acceleration") { type = IntegrationScheme::_acceleration; } else if (str == "damage") { type = IntegrationScheme::_damage; } else { stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ void IntegrationScheme::store() { for (auto data : enumerate(u_store)) { auto o = std::get<0>(data); auto & u_store = std::get<1>(data); auto & u_o = dof_manager.getDOFsDerivatives(dof_id, o); if (not u_store) { u_store = std::make_unique>( u_o, "integration_scheme_store:" + dof_id + ":" + std::to_string(o)); } else { u_store->copy(u_o); } } } /* -------------------------------------------------------------------------- */ void IntegrationScheme::restore() { for (auto o : arange(order)) { auto & u_o = dof_manager.getDOFsDerivatives(dof_id, o); u_o.copy(*u_store[o]); } } } // namespace akantu diff --git a/src/model/common/integration_scheme/integration_scheme.hh b/src/model/common/integration_scheme/integration_scheme.hh index f9c7eecd8..c6706df0a 100644 --- a/src/model/common/integration_scheme/integration_scheme.hh +++ b/src/model/common/integration_scheme/integration_scheme.hh @@ -1,125 +1,123 @@ /** * @file integration_scheme.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Jan 31 2018 * * @brief This class is just a base class for the integration schemes * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_INTEGRATION_SCHEME_HH_ #define AKANTU_INTEGRATION_SCHEME_HH_ namespace akantu { class DOFManager; } namespace akantu { class IntegrationScheme : public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: enum SolutionType { _not_defined = -1, _displacement = 0, _temperature = 0, - _pressure = 0, _damage = 0, _velocity = 1, _temperature_rate = 1, - _pressure_rate = 1, _acceleration = 2, }; IntegrationScheme(DOFManager & dof_manager, const ID & dof_id, UInt order); ~IntegrationScheme() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// generic interface of a predictor virtual void predictor(Real delta_t) = 0; /// generic interface of a corrector virtual void corrector(const SolutionType & type, Real delta_t) = 0; /// assemble the jacobian matrix virtual void assembleJacobian(const SolutionType & type, Real delta_t) = 0; /// assemble the residual virtual void assembleResidual(bool is_lumped) = 0; /// returns a list of needed matrices virtual std::vector getNeededMatrixList() = 0; /// store dofs info (beginning of steps) virtual void store(); /// restore dofs (solve failed) virtual void restore(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the order of the integration scheme UInt getOrder() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// The underlying DOFManager DOFManager & dof_manager; /// The id of the dof treated by this integration scheme. ID dof_id; /// The order of the integrator UInt order; /// last release of M matrix UInt m_release{UInt(-1)}; /// stores the values at begining of solve std::vector>> u_store; }; /* -------------------------------------------------------------------------- */ // std::ostream & operator<<(std::ostream & stream, // const IntegrationScheme::SolutionType & type); std::istream & operator>>(std::istream & stream, IntegrationScheme::SolutionType & type); /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* AKANTU_INTEGRATION_SCHEME_HH_ */ diff --git a/src/model/fluid_diffusion/fluid_diffusion_model.cc b/src/model/fluid_diffusion/fluid_diffusion_model.cc deleted file mode 100644 index c876c1f7c..000000000 --- a/src/model/fluid_diffusion/fluid_diffusion_model.cc +++ /dev/null @@ -1,1273 +0,0 @@ -/** - * @file fluid_diffusion_model.hh - * - * @author Emil Gallyamov - * - * @date creation: Aug 2019 - * - * @brief Implementation of Transient fluid diffusion model - * - * @section LICENSE - * - * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -#include "fluid_diffusion_model.hh" -#include "dumpable_inline_impl.hh" -#include "element_synchronizer.hh" -#include "fe_engine_template.hh" -#include "generalized_trapezoidal.hh" -#include "group_manager_inline_impl.hh" -#include "integrator_gauss.hh" -#include "mesh.hh" -#include "parser.hh" -#include "shape_lagrange.hh" - -#ifdef AKANTU_USE_IOHELPER -#include "dumper_element_partition.hh" -#include "dumper_elemental_field.hh" -#include "dumper_internal_material_field.hh" -#include "dumper_iohelper_paraview.hh" -#endif - -/* -------------------------------------------------------------------------- */ -namespace akantu { - -namespace fluid_diffusion { - namespace details { - class ComputeRhoFunctor { - public: - ComputeRhoFunctor(const FluidDiffusionModel & model) : model(model){}; - - void operator()(Matrix & rho, const Element & element) const { - const auto & type = element.type; - const auto & gt = element.ghost_type; - const auto & aperture = model.getApertureOnQpoints(type, gt); - auto nb_quad_points = aperture.getNbComponent(); - Real aperture_av = 0.; - for (auto q : arange(nb_quad_points)) { - aperture_av += aperture(element.element, q); - } - aperture_av /= nb_quad_points; - if (model.isModelVelocityDependent()) { - rho.set(model.getCompressibility() * aperture_av); - } else { - rho.set((model.getCompressibility() + model.getPushability()) * - aperture_av); - } - } - - private: - const FluidDiffusionModel & model; - }; - } // namespace details -} // namespace fluid_diffusion - -/* -------------------------------------------------------------------------- */ -FluidDiffusionModel::FluidDiffusionModel(Mesh & mesh, UInt dim, const ID & id) - : Model(mesh, ModelType::_fluid_diffusion_model, dim, id), - pressure_gradient("pressure_gradient", id), - pressure_on_qpoints("pressure_on_qpoints", id), - delta_pres_on_qpoints("delta_pres_on_qpoints", id), - permeability_on_qpoints("permeability_on_qpoints", id), - aperture_on_qpoints("aperture_on_qpoints", id), - prev_aperture_on_qpoints("prev_aperture_on_qpoints", id), - k_gradp_on_qpoints("k_gradp_on_qpoints", id) { - AKANTU_DEBUG_IN(); - - // mesh.registerEventHandler(*this, akantu::_ehp_lowest); - - this->spatial_dimension = std::max(1, int(mesh.getSpatialDimension()) - 1); - - this->initDOFManager(); - - this->registerDataAccessor(*this); - - if (this->mesh.isDistributed()) { - auto & synchronizer = this->mesh.getElementSynchronizer(); - this->registerSynchronizer(synchronizer, SynchronizationTag::_fdm_pressure); - this->registerSynchronizer(synchronizer, - SynchronizationTag::_fdm_gradient_pressure); - } - - registerFEEngineObject(id + ":fem", mesh, spatial_dimension); - -#ifdef AKANTU_USE_IOHELPER - this->mesh.registerDumper("fluid_diffusion", id, true); - this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); -#endif - - this->registerParam("viscosity", viscosity, 1., _pat_parsmod); - this->registerParam("compressibility", compressibility, 1., _pat_parsmod); - this->registerParam("default_aperture", default_aperture, 1.e-5, - _pat_parsmod); - this->registerParam("use_aperture_speed", use_aperture_speed, _pat_parsmod); - this->registerParam("pushability", pushability, 1., _pat_parsmod); - this->registerParam("insertion_damage", insertion_damage, 0., _pat_parsmod); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::initModel() { - auto & fem = this->getFEEngine(); - fem.initShapeFunctions(_not_ghost); - fem.initShapeFunctions(_ghost); - - pressure_on_qpoints.initialize(fem, _nb_component = 1); - delta_pres_on_qpoints.initialize(fem, _nb_component = 1); - pressure_gradient.initialize(fem, _nb_component = 1); - permeability_on_qpoints.initialize(fem, _nb_component = 1); - aperture_on_qpoints.initialize(fem, _nb_component = 1); - if (use_aperture_speed) { - prev_aperture_on_qpoints.initialize(fem, _nb_component = 1); - } - k_gradp_on_qpoints.initialize(fem, _nb_component = 1); -} - -/* -------------------------------------------------------------------------- */ -FEEngine & FluidDiffusionModel::getFEEngineBoundary(const ID & name) { - return aka::as_type(getFEEngineClassBoundary(name)); -} -/* -------------------------------------------------------------------------- */ -FluidDiffusionModel::~FluidDiffusionModel() = default; - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleCapacityLumped(const GhostType & ghost_type) { - AKANTU_DEBUG_IN(); - - auto & fem = getFEEngineClass(); - fluid_diffusion::details::ComputeRhoFunctor compute_rho(*this); - - for (const auto & type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - fem.assembleFieldLumped(compute_rho, "M", "pressure", this->getDOFManager(), - type, ghost_type); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -MatrixType FluidDiffusionModel::getMatrixType(const ID & matrix_id) { - if (matrix_id == "K" or matrix_id == "M") { - return _symmetric; - } - - return _mt_not_defined; -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleMatrix(const ID & matrix_id) { - if (matrix_id == "K") { - this->assemblePermeabilityMatrix(); - } else if (matrix_id == "M" and need_to_reassemble_capacity) { - this->assembleCapacity(); - } -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleLumpedMatrix(const ID & matrix_id) { - if (matrix_id == "M" and need_to_reassemble_capacity) { - this->assembleCapacityLumped(); - } -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleResidual() { - AKANTU_DEBUG_IN(); - - this->assembleInternalFlux(); - - this->getDOFManager().assembleToResidual("pressure", *this->external_flux, 1); - this->getDOFManager().assembleToResidual("pressure", *this->internal_flux, 1); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::predictor() { ++pressure_release; } - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::afterSolveStep(bool converged) { - AKANTU_DEBUG_IN(); - - if (not converged) { - return; - } - - ++solution_release; - need_to_reassemble_capacity_lumped = true; - need_to_reassemble_capacity = true; - - /// interpolating pressures on integration points for further use by BC - const GhostType gt = _not_ghost; - for (auto && type : mesh.elementTypes(spatial_dimension, gt, _ek_regular)) { - Array tmp(pressure_on_qpoints(type, gt)); - this->getFEEngine().interpolateOnIntegrationPoints( - *pressure, pressure_on_qpoints(type, gt), 1, type, gt); - delta_pres_on_qpoints(type, gt) = pressure_on_qpoints(type, gt); - delta_pres_on_qpoints(type, gt) -= tmp; - } - AKANTU_DEBUG_OUT(); -} -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleCapacityLumped() { - AKANTU_DEBUG_IN(); - - if (!this->getDOFManager().hasLumpedMatrix("M")) { - this->getDOFManager().getNewLumpedMatrix("M"); - } - - this->getDOFManager().getLumpedMatrix("M").zero(); - - assembleCapacityLumped(_not_ghost); - assembleCapacityLumped(_ghost); - - need_to_reassemble_capacity_lumped = false; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::initSolver( - TimeStepSolverType time_step_solver_type, - NonLinearSolverType /*non_linear_solver_type*/) { - DOFManager & dof_manager = this->getDOFManager(); - - this->allocNodalField(this->pressure, 1, "pressure"); - this->allocNodalField(this->external_flux, 1, "external_flux"); - this->allocNodalField(this->internal_flux, 1, "internal_flux"); - this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); - // this->allocNodalField(this->aperture, "aperture"); - - if (!dof_manager.hasDOFs("pressure")) { - dof_manager.registerDOFs("pressure", *this->pressure, _dst_nodal); - dof_manager.registerBlockedDOFs("pressure", *this->blocked_dofs); - } - - if (time_step_solver_type == TimeStepSolverType::_dynamic || - time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { - this->allocNodalField(this->pressure_rate, 1, "pressure_rate"); - - if (!dof_manager.hasDOFsDerivatives("pressure", 1)) { - dof_manager.registerDOFsDerivative("pressure", 1, *this->pressure_rate); - } - } -} -/* -------------------------------------------------------------------------- */ -std::tuple -FluidDiffusionModel::getDefaultSolverID(const AnalysisMethod & method) { - switch (method) { - case _explicit_lumped_mass: { - return std::make_tuple("explicit_lumped", - TimeStepSolverType::_dynamic_lumped); - } - case _static: { - return std::make_tuple("static", TimeStepSolverType::_static); - } - case _implicit_dynamic: { - return std::make_tuple("implicit", TimeStepSolverType::_dynamic); - } - default: - return std::make_tuple("unknown", TimeStepSolverType::_not_defined); - } -} - -/* -------------------------------------------------------------------------- */ -ModelSolverOptions FluidDiffusionModel::getDefaultSolverOptions( - const TimeStepSolverType & type) const { - ModelSolverOptions options; - - switch (type) { - case TimeStepSolverType::_dynamic_lumped: { - options.non_linear_solver_type = NonLinearSolverType::_lumped; - options.integration_scheme_type["pressure"] = - IntegrationSchemeType::_forward_euler; - options.solution_type["pressure"] = IntegrationScheme::_pressure_rate; - break; - } - case TimeStepSolverType::_static: { - options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; - options.integration_scheme_type["pressure"] = - IntegrationSchemeType::_pseudo_time; - options.solution_type["pressure"] = IntegrationScheme::_not_defined; - break; - } - case TimeStepSolverType::_dynamic: { - if (this->method == _explicit_consistent_mass) { - options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; - options.integration_scheme_type["pressure"] = - IntegrationSchemeType::_forward_euler; - options.solution_type["pressure"] = IntegrationScheme::_pressure_rate; - } else { - options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; - options.integration_scheme_type["pressure"] = - IntegrationSchemeType::_backward_euler; - options.solution_type["pressure"] = IntegrationScheme::_pressure; - } - break; - } - default: - AKANTU_EXCEPTION(type << " is not a valid time step solver type"); - } - - return options; -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assemblePermeabilityMatrix() { - AKANTU_DEBUG_IN(); - - this->computePermeabilityOnQuadPoints(_not_ghost); - if (permeability_release[_not_ghost] == permeability_matrix_release) { - return; - } - - if (!this->getDOFManager().hasMatrix("K")) { - this->getDOFManager().getNewMatrix("K", getMatrixType("K")); - } - this->getDOFManager().getMatrix("K").zero(); - - this->assemblePermeabilityMatrix(_not_ghost); - // switch (mesh.getSpatialDimension()) { - // case 1: - // this->assemblePermeabilityMatrix<1>(_not_ghost); - // break; - // case 2: - // this->assemblePermeabilityMatrix<2>(_not_ghost); - // break; - // case 3: - // this->assemblePermeabilityMatrix<3>(_not_ghost); - // break; - // } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assemblePermeabilityMatrix( - const GhostType & ghost_type) { - AKANTU_DEBUG_IN(); - - auto & fem = this->getFEEngine(); - - for (auto && type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - auto nb_element = mesh.getNbElement(type, ghost_type); - auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); - - auto bt_d_b = std::make_unique>( - nb_element * nb_quadrature_points, - nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); - - fem.computeBtDB(permeability_on_qpoints(type, ghost_type), *bt_d_b, 2, type, - ghost_type); - - /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ - auto K_e = std::make_unique>( - nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); - - fem.integrate(*bt_d_b, *K_e, nb_nodes_per_element * nb_nodes_per_element, - type, ghost_type); - - this->getDOFManager().assembleElementalMatricesToMatrix( - "K", "pressure", *K_e, type, ghost_type, _symmetric); - } - - permeability_matrix_release = permeability_release[ghost_type]; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::computePermeabilityOnQuadPoints( - const GhostType & ghost_type) { - - // if already computed once check if need to compute - if (not initial_permeability[ghost_type]) { - // if aperture did not change, permeability will not vary - if (solution_release == permeability_release[ghost_type]) { - return; - } - } - - for (const auto & type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - - // auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); - // auto nb_element = mesh.getNbElement(type); - // Array aperture_interpolated(nb_element * nb_nodes_per_element, - // 1); - - // // compute the pressure on quadrature points - // this->getFEEngine().interpolateOnIntegrationPoints( - // *aperture, aperture_interpolated, 1, type, ghost_type); - - auto & perm = permeability_on_qpoints(type, ghost_type); - auto & aperture = aperture_on_qpoints(type, ghost_type); - for (auto && tuple : zip(perm, aperture)) { - auto & k = std::get<0>(tuple); - auto & aper = std::get<1>(tuple); - k = aper * aper * aper / (12 * this->viscosity); - } - } - - permeability_release[ghost_type] = solution_release; - initial_permeability[ghost_type] = false; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::computeKgradP(const GhostType & ghost_type) { - computePermeabilityOnQuadPoints(ghost_type); - - for (const auto & type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - auto & gradient = pressure_gradient(type, ghost_type); - this->getFEEngine().gradientOnIntegrationPoints(*pressure, gradient, 1, - type, ghost_type); - - for (auto && values : zip(permeability_on_qpoints(type, ghost_type), - gradient, k_gradp_on_qpoints(type, ghost_type))) { - const auto & k = std::get<0>(values); - const auto & BP = std::get<1>(values); - auto & k_BP = std::get<2>(values); - - k_BP = k * BP; - } - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleInternalFlux() { - AKANTU_DEBUG_IN(); - - this->internal_flux->clear(); - - this->synchronize(SynchronizationTag::_fdm_pressure); - auto & fem = this->getFEEngine(); - - for (auto ghost_type : ghost_types) { - // compute k \grad P - computeKgradP(ghost_type); - - for (auto type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - - auto & k_gradp_on_qpoints_vect = k_gradp_on_qpoints(type, ghost_type); - - UInt nb_quad_points = k_gradp_on_qpoints_vect.size(); - Array bt_k_gP(nb_quad_points, nb_nodes_per_element); - fem.computeBtD(k_gradp_on_qpoints_vect, bt_k_gP, type, ghost_type); - - if (this->use_aperture_speed) { - auto aper_speed = aperture_on_qpoints(type, ghost_type); - aper_speed -= prev_aperture_on_qpoints(type, ghost_type); - aper_speed *= 1 / this->time_step; - Array aper_speed_by_shapes(nb_quad_points, nb_nodes_per_element); - - fem.computeNtb(aper_speed, aper_speed_by_shapes, type, ghost_type); - bt_k_gP += aper_speed_by_shapes; - } - UInt nb_elements = mesh.getNbElement(type, ghost_type); - Array int_bt_k_gP(nb_elements, nb_nodes_per_element); - - fem.integrate(bt_k_gP, int_bt_k_gP, nb_nodes_per_element, type, - ghost_type); - - this->getDOFManager().assembleElementalArrayLocalArray( - int_bt_k_gP, *this->internal_flux, type, ghost_type, -1); - } - } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -Real FluidDiffusionModel::getStableTimeStep() { - AKANTU_DEBUG_IN(); - - Real el_size; - Real min_el_size = std::numeric_limits::max(); - Real max_aper_on_qpoint = std::numeric_limits::min(); - - for (const auto & type : - mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { - - auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); - auto & aper = aperture_on_qpoints(type, _not_ghost); - auto nb_quad_per_element = aper.getNbComponent(); - auto mesh_dim = this->mesh.getSpatialDimension(); - Array coord(0, nb_nodes_per_element * mesh_dim); - FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, - _not_ghost); - - for (auto && data : zip(make_view(coord, mesh_dim, nb_nodes_per_element), - make_view(aper, nb_quad_per_element))) { - Matrix & el_coord = std::get<0>(data); - Vector & el_aper = std::get<1>(data); - el_size = getFEEngine().getElementInradius(el_coord, type); - min_el_size = std::min(min_el_size, el_size); - max_aper_on_qpoint = std::max(max_aper_on_qpoint, el_aper.norm()); - } - - AKANTU_DEBUG_INFO("The minimum element size : " - << min_el_size - << " and the max aperture is : " << max_aper_on_qpoint); - } - - Real min_dt = 2. * min_el_size * min_el_size / 4 * this->compressibility / - max_aper_on_qpoint / max_aper_on_qpoint; - - mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); - - AKANTU_DEBUG_OUT(); - - return min_dt; -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::setTimeStep(Real dt, const ID & solver_id) { - Model::setTimeStep(time_step, solver_id); - this->time_step = dt; -#if defined(AKANTU_USE_IOHELPER) - // this->mesh.getDumper("fluid_diffusion").setTimeStep(time_step); -#endif -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::readMaterials() { - auto sect = this->getParserSection(); - - if (not std::get<1>(sect)) { - const auto & section = std::get<0>(sect); - this->parseSection(section); - } -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::initFullImpl(const ModelOptions & options) { - readMaterials(); - Model::initFullImpl(options); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::assembleCapacity() { - AKANTU_DEBUG_IN(); - auto ghost_type = _not_ghost; - - this->getDOFManager().getMatrix("M").zero(); - - auto & fem = getFEEngineClass(); - - fluid_diffusion::details::ComputeRhoFunctor rho_functor(*this); - - for (auto && type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { - fem.assembleFieldMatrix(rho_functor, "M", "pressure", this->getDOFManager(), - type, ghost_type); - } - - need_to_reassemble_capacity = false; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::computeRho(Array & rho, ElementType type, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - rho.copy(aperture_on_qpoints(type, ghost_type)); - if (this->use_aperture_speed) { - rho *= this->compressibility; - } else { - rho *= (this->compressibility + this->pushability); - } - - AKANTU_DEBUG_OUT(); -} // namespace akantu - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::onElementsAdded(const Array & element_list, - const NewElementsEvent & /*event*/) { - AKANTU_DEBUG_IN(); - if (element_list.empty()) { - return; - } - /// TODO had to do this ugly way because the meeting is in 3 days - // removing it will make the assemble mass matrix fail. jacobians are not - // updated - auto & fem = this->getFEEngine(); - fem.initShapeFunctions(_not_ghost); - fem.initShapeFunctions(_ghost); - - this->resizeFields(); - - auto type_fluid = - *mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular).begin(); - const auto nb_quad_elem = - getFEEngine().getNbIntegrationPoints(type_fluid, _not_ghost); - auto aperture_it = - make_view(getApertureOnQpoints(type_fluid), nb_quad_elem).begin(); - - /// assign default aperture to newly added nodes - for (auto && element : element_list) { - auto elem_id = element.element; - for (auto ip : arange(nb_quad_elem)) { - aperture_it[elem_id](ip) = this->default_aperture; - } - } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::resizeFields() { - - AKANTU_DEBUG_IN(); - - /// elemental fields - auto & fem = this->getFEEngine(); - pressure_on_qpoints.initialize(fem, _nb_component = 1); - delta_pres_on_qpoints.initialize(fem, _nb_component = 1); - pressure_gradient.initialize(fem, _nb_component = 1); - permeability_on_qpoints.initialize(fem, _nb_component = 1); - aperture_on_qpoints.initialize(fem, _nb_component = 1); - if (use_aperture_speed) { - prev_aperture_on_qpoints.initialize(fem, _nb_component = 1); - } - k_gradp_on_qpoints.initialize(fem, _nb_component = 1); - - /// nodal fields - UInt nb_nodes = mesh.getNbNodes(); - - if (pressure != nullptr) { - pressure->resize(nb_nodes, 0.); - ++solution_release; - } - if (external_flux != nullptr) { - external_flux->resize(nb_nodes, 0.); - } - if (internal_flux != nullptr) { - internal_flux->resize(nb_nodes, 0.); - } - if (blocked_dofs != nullptr) { - blocked_dofs->resize(nb_nodes, false); - } - if (pressure_rate != nullptr) { - pressure_rate->resize(nb_nodes, 0.); - } - need_to_reassemble_capacity_lumped = true; - need_to_reassemble_capacity = true; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -#if defined(AKANTU_COHESIVE_ELEMENT) -void FluidDiffusionModel::getApertureOnQpointsFromCohesive( - const SolidMechanicsModelCohesive & coh_model, bool first_time) { - AKANTU_DEBUG_IN(); - - // get element types - auto & coh_mesh = coh_model.getMesh(); - const auto dim = coh_mesh.getSpatialDimension(); - const GhostType gt = _not_ghost; - auto type = *coh_mesh.elementTypes(dim, gt, _ek_regular).begin(); - const ElementType type_facets = Mesh::getFacetType(type); - const ElementType typecoh = FEEngine::getCohesiveElementType(type_facets); - // const auto nb_coh_elem = coh_mesh.getNbElement(typecoh); - const auto nb_quad_coh_elem = coh_model.getFEEngine("CohesiveFEEngine") - .getNbIntegrationPoints(typecoh, gt); - - auto type_fluid = - *mesh.elementTypes(spatial_dimension, gt, _ek_regular).begin(); - const auto nb_fluid_elem = mesh.getNbElement(type_fluid); - const auto nb_quad_elem = - getFEEngine().getNbIntegrationPoints(type_fluid, gt); - // AKANTU_DEBUG_ASSERT(nb_quad_coh_elem == nb_quad_elem, - // "Different number of integration points per cohesive " - // "and fluid elements"); - - // save previous aperture - if (not first_time and use_aperture_speed) { - prev_aperture_on_qpoints.copy(aperture_on_qpoints); - } - - // AKANTU_DEBUG_ASSERT(nb_coh_elem == nb_fluid_elem, - // "Different number of cohesive and fluid elements"); - - auto aperture_it = - make_view(getApertureOnQpoints(type_fluid), nb_quad_elem).begin(); - - /// loop on each segment element - for (auto && element_nb : arange(nb_fluid_elem)) { - Element element{type_fluid, element_nb, gt}; - auto global_coh_elem = - mesh.getElementalData("cohesive_elements")(element); - auto ge = global_coh_elem.element; - auto le = coh_model.getMaterialLocalNumbering( - global_coh_elem.type, global_coh_elem.ghost_type)(ge); - auto model_mat_index = coh_model.getMaterialByElement( - global_coh_elem.type, global_coh_elem.ghost_type)(ge); - const auto & material = coh_model.getMaterial(model_mat_index); - - // /// access cohesive material - // const auto & mat = dynamic_cast(material); - // const auto & normal_open_norm = mat.getNormalOpeningNorm(typecoh, gt); - const auto normal_open_norm_it = - make_view(material.getArray("normal_opening_norm", typecoh), - nb_quad_coh_elem) - .begin(); - if (nb_quad_elem == 1) { - /// coh el quad points are averaged to assign single opening to flow el - Real aperture_av = 0.; - for (UInt ip : arange(nb_quad_coh_elem)) { - auto normal_open_norm_ip = normal_open_norm_it[le](ip); - aperture_av += normal_open_norm_ip; - } - aperture_av /= nb_quad_coh_elem; - auto & aperture_ip = aperture_it[ge](0); - if (aperture_av < this->default_aperture) { - aperture_ip = this->default_aperture; - } else { - aperture_ip = aperture_av; - } - } else { - /// each flow el quad point gets aperture from corresponding coh el qp - for (UInt ip : arange(nb_quad_coh_elem)) { - auto normal_open_norm_ip = normal_open_norm_it[le](ip); - auto & aperture_ip = aperture_it[ge](ip); - if (normal_open_norm_ip < this->default_aperture) { - aperture_ip = this->default_aperture; - } else { - aperture_ip = normal_open_norm_ip; - } - } - } - } - - // save previous aperture - if (first_time && use_aperture_speed) { - prev_aperture_on_qpoints.copy(aperture_on_qpoints); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::updateFluidElementsFromCohesive( - const SolidMechanicsModelCohesive & coh_model) { - AKANTU_DEBUG_IN(); - - // get element types - auto & coh_mesh = coh_model.getMesh(); - const auto dim = coh_mesh.getSpatialDimension(); - const GhostType gt = _not_ghost; - auto type = *coh_mesh.elementTypes(dim, gt, _ek_regular).begin(); - const ElementType type_facets = Mesh::getFacetType(type); - const ElementType typecoh = FEEngine::getCohesiveElementType(type_facets); - const auto nb_coh_elem = coh_mesh.getNbElement(typecoh); - const auto nb_quad_coh_elem = coh_model.getFEEngine("CohesiveFEEngine") - .getNbIntegrationPoints(typecoh, gt); - - auto type_fluid = - *mesh.elementTypes(spatial_dimension, gt, _ek_regular).begin(); - const auto nb_fluid_elem = mesh.getNbElement(type_fluid); - - // check if there is need to add more fluid elements - if (nb_coh_elem == nb_fluid_elem) { - return; - } - - NewElementsEvent new_facets; - NewNodesEvent new_nodes; - NewElementsEvent new_fluid_facets; - - /// loop on each cohesive element and check its damage - for (auto && coh_el : arange(nb_coh_elem)) { - - // check if this cohesive element is present in fluid mesh - auto & existing_coh_elements = - mesh.getData("cohesive_elements", type_facets, gt); - Element coh_element = {typecoh, coh_el, gt}; - auto it = existing_coh_elements.find(coh_element); - if (it != UInt(-1)) { - continue; - } - - // check damage at the cohesive element - // Element element{type_fluid, element_nb, gt}; - auto le = coh_model.getMaterialLocalNumbering(typecoh, gt)(coh_el); - auto model_mat_index = coh_model.getMaterialByElement(typecoh, gt)(coh_el); - const auto & material = coh_model.getMaterial(model_mat_index); - const auto damage_it = - make_view(material.getArray("damage", typecoh), nb_quad_coh_elem) - .begin(); - Real damage_av = 0.; - for (auto ip : arange(nb_quad_coh_elem)) { - auto damage_ip = damage_it[le](ip); - damage_av += damage_ip; - } - damage_av /= nb_quad_coh_elem; - - if (damage_av < this->insertion_damage) { - continue; - } - - // add fluid element to connectivity and corresponding nodes - Array & nodes_added = new_nodes.getList(); - auto & crack_facets = coh_mesh.getElementGroup("crack_facets"); - auto & nodes = mesh.getNodes(); - auto && coh_nodes = coh_mesh.getNodes(); - auto && coh_conn = coh_mesh.getConnectivity(typecoh, gt); - Vector conn(coh_conn.getNbComponent() / 2); - auto coh_facet_conn_it = make_view(coh_mesh.getConnectivity(typecoh, gt), - coh_conn.getNbComponent() / 2) - .begin(); - auto cohesive_nodes_first_half = coh_facet_conn_it[2 * coh_el]; - for (auto c : akantu::arange(conn.size())) { - auto coh_node = cohesive_nodes_first_half(c); - Vector pos(mesh.getSpatialDimension()); - for (auto && data : enumerate(pos)) { - std::get<1>(data) = coh_nodes(coh_node, std::get<0>(data)); - } - auto idx = nodes.find(pos); - if (idx == UInt(-1)) { - nodes.push_back(pos); - conn(c) = nodes.size() - 1; - nodes_added.push_back(nodes.size() - 1); - } else { - conn(c) = idx; - } - } - mesh.getData("cohesive_elements", type_facets, gt) - .push_back(coh_element); - mesh.getConnectivity(type_facets, gt).push_back(conn); - Element fluid_facet{type_facets, - mesh.getConnectivity(type_facets, gt).size() - 1, gt}; - new_fluid_facets.getList().push_back(fluid_facet); - - // instantiate connectivity type for this facet type if not yet done - if (not coh_mesh.getConnectivities().exists(type_facets, gt)) { - coh_mesh.addConnectivityType(type_facets, gt); - } - - // add two facets per cohesive element into cohesive mesh - auto & facet_conn_mesh = coh_mesh.getConnectivity(type_facets, gt); - for (auto f : akantu::arange(2)) { - Vector coh_nodes_one_side = coh_facet_conn_it[2 * coh_el + f]; - facet_conn_mesh.push_back(coh_nodes_one_side); - Element new_facet{type_facets, facet_conn_mesh.size() - 1, gt}; - crack_facets.add(new_facet); - new_facets.getList().push_back(new_facet); - } - } - mesh.sendEvent(new_nodes); - mesh.sendEvent(new_fluid_facets); - MeshUtils::fillElementToSubElementsData(coh_mesh); - coh_mesh.sendEvent(new_facets); - - AKANTU_DEBUG_OUT(); -} - -#endif - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::applyExternalFluxAtElementGroup( - const Real & rate, const ElementGroup & source_facets, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - for (auto && type : - source_facets.elementTypes(spatial_dimension, ghost_type)) { - const auto & element_ids = source_facets.getElements(type, ghost_type); - - const auto nb_fluid_elem = mesh.getNbElement(type); - const auto type_size = element_ids.size(); - AKANTU_DEBUG_ASSERT(type_size <= nb_fluid_elem, - "Number of provided source facets exceeds total number " - "of flow elements"); - const auto fluid_conn_it = - make_view(mesh.getConnectivity(type, ghost_type), - mesh.getConnectivity(type, ghost_type).getNbComponent()) - .begin(); - auto & source = getExternalFlux(); - source.clear(); - - /// loop on each element of the group - for (auto el : element_ids) { - AKANTU_DEBUG_ASSERT(el <= nb_fluid_elem, - "Element number in the source facets group exceeds " - "maximum number of flow elements"); - - const Vector fluid_conn = fluid_conn_it[el]; - for (const auto & node : fluid_conn) { - source(node) += rate / (fluid_conn.size() * element_ids.size()); - } - } - } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -Array FluidDiffusionModel::getQpointsCoord() { - AKANTU_DEBUG_IN(); - - const auto & nodes_coords = mesh.getNodes(); - const GhostType gt = _not_ghost; - auto type_fluid = - *mesh.elementTypes(spatial_dimension, gt, _ek_regular).begin(); - const auto nb_fluid_elem = mesh.getNbElement(type_fluid); - auto & fem = this->getFEEngine(); - const auto nb_quad_points = fem.getNbIntegrationPoints(type_fluid, gt); - - Array quad_coords(nb_fluid_elem * nb_quad_points, - spatial_dimension + 1); - fem.interpolateOnIntegrationPoints(nodes_coords, quad_coords, - spatial_dimension + 1, type_fluid, gt); - AKANTU_DEBUG_OUT(); - return quad_coords; -} - -/* -------------------------------------------------------------------------- */ -/* -------------------------------------------------------------------------- */ -#ifdef AKANTU_USE_IOHELPER - -std::shared_ptr FluidDiffusionModel::createNodalFieldBool( - const std::string & field_name, const std::string & group_name, - __attribute__((unused)) bool padding_flag) { - - std::map *> uint_nodal_fields; - uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); - - auto field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); - return field; -} - -/* -------------------------------------------------------------------------- */ -std::shared_ptr FluidDiffusionModel::createNodalFieldReal( - const std::string & field_name, const std::string & group_name, - __attribute__((unused)) bool padding_flag) { - - if (field_name == "capacity_lumped") { - AKANTU_EXCEPTION( - "Capacity lumped is a nodal field now stored in the DOF manager." - "Therefore it cannot be used by a dumper anymore"); - } - - std::map *> real_nodal_fields; - real_nodal_fields["pressure"] = pressure.get(); - real_nodal_fields["pressure_rate"] = pressure_rate.get(); - real_nodal_fields["external_flux"] = external_flux.get(); - real_nodal_fields["internal_flux"] = internal_flux.get(); - // real_nodal_fields["increment"] = increment.get(); - - std::shared_ptr field = - mesh.createNodalField(real_nodal_fields[field_name], group_name); - - return field; -} - -/* -------------------------------------------------------------------------- */ -std::shared_ptr FluidDiffusionModel::createElementalField( - const std::string & field_name, const std::string & group_name, - bool /*padding_flag*/, UInt /*spatial_dimension*/, - ElementKind element_kind) { - - std::shared_ptr field; - - if (field_name == "partitions") { - field = mesh.createElementalField( - mesh.getConnectivities(), group_name, this->spatial_dimension, - element_kind); - } else if (field_name == "pressure_gradient") { - ElementTypeMap nb_data_per_elem = - this->mesh.getNbDataPerElem(pressure_gradient); - - field = mesh.createElementalField( - pressure_gradient, group_name, this->spatial_dimension, element_kind, - nb_data_per_elem); - } else if (field_name == "permeability") { - ElementTypeMap nb_data_per_elem = - this->mesh.getNbDataPerElem(permeability_on_qpoints); - - field = mesh.createElementalField( - permeability_on_qpoints, group_name, this->spatial_dimension, - element_kind, nb_data_per_elem); - } else if (field_name == "aperture") { - ElementTypeMap nb_data_per_elem = - this->mesh.getNbDataPerElem(aperture_on_qpoints); - - field = mesh.createElementalField( - aperture_on_qpoints, group_name, this->spatial_dimension, element_kind, - nb_data_per_elem); - } else if (field_name == "pressure_on_qpoints") { - ElementTypeMap nb_data_per_elem = - this->mesh.getNbDataPerElem(pressure_on_qpoints); - - field = mesh.createElementalField( - pressure_on_qpoints, group_name, this->spatial_dimension, element_kind, - nb_data_per_elem); - } - - return field; -} - -/* -------------------------------------------------------------------------- */ -#else -/* -------------------------------------------------------------------------- */ -std::shared_ptr FluidDiffusionModel::createElementalField( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag, - __attribute__((unused)) const ElementKind & element_kind) { - return nullptr; -} - -/* -------------------------------------------------------------------------- */ -std::shared_ptr FluidDiffusionModel::createNodalFieldBool( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag) { - return nullptr; -} - -/* -------------------------------------------------------------------------- */ -std::shared_ptr FluidDiffusionModel::createNodalFieldReal( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag) { - return nullptr; -} -#endif - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::dump(const std::string & dumper_name) { - mesh.dump(dumper_name); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::dump(const std::string & dumper_name, UInt step) { - mesh.dump(dumper_name, step); -} - -/* ------------------------------------------------------------------------- - */ -void FluidDiffusionModel::dump(const std::string & dumper_name, Real time, - UInt step) { - mesh.dump(dumper_name, time, step); -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::dump() { mesh.dump(); } - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::dump(UInt step) { mesh.dump(step); } - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::dump(Real time, UInt step) { mesh.dump(time, step); } - -/* -------------------------------------------------------------------------- */ -inline UInt -FluidDiffusionModel::getNbData(const Array & indexes, - const SynchronizationTag & tag) const { - AKANTU_DEBUG_IN(); - - UInt size = 0; - UInt nb_nodes = indexes.size(); - - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - size += nb_nodes * sizeof(Real); - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } - - AKANTU_DEBUG_OUT(); - return size; -} - -/* -------------------------------------------------------------------------- */ -inline void -FluidDiffusionModel::packData(CommunicationBuffer & buffer, - const Array & indexes, - const SynchronizationTag & tag) const { - AKANTU_DEBUG_IN(); - - for (auto index : indexes) { - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - buffer << (*pressure)(index); - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } - } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -inline void FluidDiffusionModel::unpackData(CommunicationBuffer & buffer, - const Array & indexes, - const SynchronizationTag & tag) { - AKANTU_DEBUG_IN(); - - for (auto index : indexes) { - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - buffer >> (*pressure)(index); - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -inline UInt -FluidDiffusionModel::getNbData(const Array & elements, - const SynchronizationTag & tag) const { - AKANTU_DEBUG_IN(); - - UInt size = 0; - UInt nb_nodes_per_element = 0; - Array::const_iterator it = elements.begin(); - Array::const_iterator end = elements.end(); - for (; it != end; ++it) { - const Element & el = *it; - nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); - } - - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - size += nb_nodes_per_element * sizeof(Real); // pressure - break; - } - case SynchronizationTag::_fdm_gradient_pressure: { - // pressure gradient - size += getNbIntegrationPoints(elements) * sizeof(Real); - size += nb_nodes_per_element * sizeof(Real); // nodal pressures - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } - - AKANTU_DEBUG_OUT(); - return size; -} - -/* -------------------------------------------------------------------------- */ -inline void -FluidDiffusionModel::packData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) const { - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - packNodalDataHelper(*pressure, buffer, elements, mesh); - break; - } - case SynchronizationTag::_fdm_gradient_pressure: { - packElementalDataHelper(pressure_gradient, buffer, elements, true, - getFEEngine()); - packNodalDataHelper(*pressure, buffer, elements, mesh); - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } -} - -/* -------------------------------------------------------------------------- */ -inline void FluidDiffusionModel::unpackData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) { - switch (tag) { - case SynchronizationTag::_fdm_pressure: { - unpackNodalDataHelper(*pressure, buffer, elements, mesh); - break; - } - case SynchronizationTag::_fdm_gradient_pressure: { - unpackElementalDataHelper(pressure_gradient, buffer, elements, true, - getFEEngine()); - unpackNodalDataHelper(*pressure, buffer, elements, mesh); - - break; - } - default: { - AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); - } - } -} - -/* -------------------------------------------------------------------------- */ -void FluidDiffusionModel::injectIntoFacetsByCoord(const Vector & position, - const Real & injection_rate) { - AKANTU_DEBUG_IN(); - - auto dim = mesh.getSpatialDimension(); - const auto gt = _not_ghost; - const auto & pos = mesh.getNodes(); - const auto pos_it = make_view(pos, dim).begin(); - auto & source = *external_flux; - - Vector bary_facet(dim); - - Real min_dist = std::numeric_limits::max(); - Element inj_facet; - for_each_element( - mesh, - [&](auto && facet) { - mesh.getBarycenter(facet, bary_facet); - auto dist = std::abs(bary_facet.distance(position)); - if (dist < min_dist) { - min_dist = dist; - inj_facet = facet; - } - }, - _spatial_dimension = dim - 1); - - auto & facet_conn = mesh.getConnectivity(inj_facet.type, gt); - auto nb_nodes_per_elem = facet_conn.getNbComponent(); - for (auto node : arange(nb_nodes_per_elem)) { - source(facet_conn(inj_facet.element, node)) = - injection_rate / nb_nodes_per_elem; - } - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -} // namespace akantu diff --git a/src/model/fluid_diffusion/fluid_diffusion_model.hh b/src/model/fluid_diffusion/fluid_diffusion_model.hh deleted file mode 100644 index bcf576dbb..000000000 --- a/src/model/fluid_diffusion/fluid_diffusion_model.hh +++ /dev/null @@ -1,381 +0,0 @@ -/** - * @file fluid_diffusion_model.hh - * - * @author Emil Gallyamov - * - * @date creation: Aug 2019 - * - * @brief Transient fluid diffusion model based on the Heat Transfer Model - * - * @section LICENSE - * - * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Akantu. If not, see . - * - */ - -/* -------------------------------------------------------------------------- */ -#include "data_accessor.hh" -#include "fe_engine.hh" -#include "fluid_diffusion_model_event_handler.hh" -#include "model.hh" -#if defined(AKANTU_COHESIVE_ELEMENT) -#include "material_cohesive.hh" -#include "solid_mechanics_model_cohesive.hh" -#endif - -/* -------------------------------------------------------------------------- */ -#include -/* -------------------------------------------------------------------------- */ - -#ifndef __AKANTU_FLUID_DIFFUSION_MODEL_HH__ -#define __AKANTU_FLUID_DIFFUSION_MODEL_HH__ - -namespace akantu { -template -class IntegratorGauss; -template class ShapeLagrange; -} // namespace akantu - -namespace akantu { - -class FluidDiffusionModel - : public Model, - public DataAccessor, - public DataAccessor, - public EventHandlerManager { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - using FEEngineType = FEEngineTemplate; - - FluidDiffusionModel(Mesh & mesh, UInt dim = _all_dimensions, - const ID & id = "fluid_diffusion_model"); - - ~FluidDiffusionModel() override; - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -protected: - /// generic function to initialize everything ready for explicit dynamics - void initFullImpl(const ModelOptions & options) override; - - /// read one material file to instantiate all the materials - void readMaterials(); - - /// allocate all vectors - void initSolver(TimeStepSolverType, NonLinearSolverType) override; - - /// initialize the model - void initModel() override; - - void predictor() override; - - /// callback for the solver, this is called at end of solve - void afterSolveStep(bool is_converged = false) override; - - /// compute the heat flux - void assembleResidual() override; - - /// get the type of matrix needed - MatrixType getMatrixType(const ID &) override; - - /// callback to assemble a Matrix - void assembleMatrix(const ID &) override; - - /// callback to assemble a lumped Matrix - void assembleLumpedMatrix(const ID &) override; - - std::tuple - getDefaultSolverID(const AnalysisMethod & method) override; - - ModelSolverOptions - getDefaultSolverOptions(const TimeStepSolverType & type) const override; - - /// resize all fields on nodes added event - void resizeFields(); - -public: - /// get coordinates of qpoints in same order as other elemental fields - Array getQpointsCoord(); - -#ifdef AKANTU_COHESIVE_ELEMENT - /// get apperture at nodes from displacement field of cohesive model - void getApertureOnQpointsFromCohesive( - const SolidMechanicsModelCohesive & coh_model, bool first_time = false); - - /// insert fluid elements based on cohesive elements and their damage - void updateFluidElementsFromCohesive( - const SolidMechanicsModelCohesive & coh_model); -#endif - - void applyExternalFluxAtElementGroup(const Real & rate, - const ElementGroup & source_facets, - GhostType ghost_type = _not_ghost); - - void injectIntoFacetsByCoord(const Vector & position, - const Real & injection_rate); - /* ------------------------------------------------------------------------ */ - /* Methods for explicit */ - /* ------------------------------------------------------------------------ */ -public: - /// compute and get the stable time step - Real getStableTimeStep(); - - /// set the stable timestep - void setTimeStep(Real dt, const ID & solver_id = "") override; - - // temporary protection to prevent bad usage: should check for bug -protected: - /// compute the internal heat flux \todo Need code review: currently not - /// public method - void assembleInternalFlux(); - -public: - /// calculate the lumped capacity vector for heat transfer problem - void assembleCapacityLumped(); - - /* ------------------------------------------------------------------------ */ - /* Mesh Event Handler inherited members */ - /* ------------------------------------------------------------------------ */ -protected: - void onElementsAdded(const Array & element_list, - const NewElementsEvent & /*event*/) override; - /* ------------------------------------------------------------------------ */ - /* Methods for static */ - /* ------------------------------------------------------------------------ */ -public: - /// assemble the conductivity matrix - void assemblePermeabilityMatrix(); - - /// assemble the conductivity matrix - void assembleCapacity(); - - /// compute the capacity on quadrature points - void computeRho(Array & rho, ElementType type, GhostType ghost_type); - -private: - /// calculate the lumped capacity vector for heat transfer problem (w - /// ghost type) - void assembleCapacityLumped(const GhostType & ghost_type); - - /// assemble the conductivity matrix (w/ ghost type) - void assemblePermeabilityMatrix(const GhostType & ghost_type); - - /// compute the conductivity tensor for each quadrature point in an array - void computePermeabilityOnQuadPoints(const GhostType & ghost_type); - - /// compute vector k \grad T for each quadrature point - void computeKgradP(const GhostType & ghost_type); - - /// compute the thermal energy - Real computeThermalEnergyByNode(); - - /* ------------------------------------------------------------------------ */ - /* Data Accessor inherited members */ - /* ------------------------------------------------------------------------ */ -public: - inline UInt getNbData(const Array & elements, - const SynchronizationTag & tag) const override; - inline void packData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) const override; - inline void unpackData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) override; - - inline UInt getNbData(const Array & indexes, - const SynchronizationTag & tag) const override; - inline void packData(CommunicationBuffer & buffer, - const Array & indexes, - const SynchronizationTag & tag) const override; - inline void unpackData(CommunicationBuffer & buffer, - const Array & indexes, - const SynchronizationTag & tag) override; - - /* ------------------------------------------------------------------------ */ - /* Dumpable interface */ - /* ------------------------------------------------------------------------ */ -public: - std::shared_ptr - createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - std::shared_ptr - createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - std::shared_ptr - createElementalField(const std::string & field_name, - const std::string & group_name, bool padding_flag, - UInt spatial_dimension, ElementKind kind) override; - - virtual void dump(const std::string & dumper_name); - virtual void dump(const std::string & dumper_name, UInt step); - virtual void dump(const std::string & dumper_name, Real time, UInt step); - - void dump() override; - - virtual void dump(UInt step); - - virtual void dump(Real time, UInt step); - - /* ------------------------------------------------------------------------ */ - /* Accessors */ - /* ------------------------------------------------------------------------ */ -public: - AKANTU_GET_MACRO(Compressibility, compressibility, Real); - AKANTU_GET_MACRO(Pushability, pushability, Real); - /// get the dimension of the system space - AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); - /// get the current value of the time step - AKANTU_GET_MACRO(TimeStep, time_step, Real); - /// get the assembled heat flux - AKANTU_GET_MACRO(InternalFlux, *internal_flux, Array &); - /// get the boundary vector - AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); - /// get the external heat rate vector - AKANTU_GET_MACRO(ExternalFlux, *external_flux, Array &); - /// get the temperature gradient - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PressureGradient, pressure_gradient, - Real); - /// get the permeability on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PermeabilityOnQpoints, - permeability_on_qpoints, Real); - /// get the pressure on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PressureOnQpoints, pressure_on_qpoints, - Real); - /// get the pressure change on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(DeltaPressureOnQpoints, - delta_pres_on_qpoints, Real); - /// get the aperture on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ApertureOnQpoints, aperture_on_qpoints, - Real); - /// get the aperture on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ApertureOnQpoints, aperture_on_qpoints, - Real); - /// internal variables - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KgradP, k_gradp_on_qpoints, Real); - /// get the temperature - AKANTU_GET_MACRO(Pressure, *pressure, Array &); - /// get the temperature derivative - AKANTU_GET_MACRO(PressureRate, *pressure_rate, Array &); - - inline bool isModelVelocityDependent() const { return use_aperture_speed; } - // /// get the energy denominated by thermal - // Real getEnergy(const std::string & energy_id, const ElementType & type, - // UInt index); - // /// get the energy denominated by thermal - // Real getEnergy(const std::string & energy_id); - - // /// get the thermal energy for a given element - // Real getThermalEnergy(const ElementType & type, UInt index); - // /// get the thermal energy for a given element - // Real getThermalEnergy(); - -protected: - /* ------------------------------------------------------------------------ */ - FEEngine & getFEEngineBoundary(const ID & name = "") override; - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -private: - /// number of iterations - UInt n_iter; - - /// time step - Real time_step; - - /// pressure array - std::unique_ptr> pressure; - - /// pressure derivatives array - std::unique_ptr> pressure_rate; - - // /// increment array (@f$\delta \dot P@f$ or @f$\delta P@f$) - // Array * increment{nullptr}; - - /// the speed of the changing temperature - ElementTypeMapArray pressure_gradient; - - /// pressure field on quadrature points - ElementTypeMapArray pressure_on_qpoints; - - /// pressure change on quadrature points - ElementTypeMapArray delta_pres_on_qpoints; - - /// conductivity tensor on quadrature points - ElementTypeMapArray permeability_on_qpoints; - - /// aperture on quadrature points - ElementTypeMapArray aperture_on_qpoints; - - /// aperture on quadrature points - ElementTypeMapArray prev_aperture_on_qpoints; - - /// vector k \grad T on quad points - ElementTypeMapArray k_gradp_on_qpoints; - - /// external flux vector - std::unique_ptr> external_flux; - - /// residuals array - std::unique_ptr> internal_flux; - - /// boundary vector - std::unique_ptr> blocked_dofs; - - // viscosity - Real viscosity{0.}; - - /// compressibility Cf = 1/pho drho/dP - Real compressibility{0.}; - - // default value of aperture on a newly added nodesynchronizer - Real default_aperture{0.}; - - bool need_to_reassemble_capacity{true}; - bool need_to_reassemble_capacity_lumped{true}; - UInt pressure_release{0}; - UInt solution_release{0}; - UInt permeability_matrix_release{0}; - - std::unordered_map initial_permeability{{_not_ghost, true}, - {_ghost, true}}; - std::unordered_map permeability_release{{_not_ghost, 0}, - {_ghost, 0}}; - bool use_aperture_speed{false}; - - // pushability Ch = 1/h dh/dP - Real pushability{0.}; - - // damage at which cohesive element will be duplicated in fluid mesh - Real insertion_damage{0.}; -}; - -} // namespace akantu - -/* -------------------------------------------------------------------------- */ -/* inline functions */ -/* -------------------------------------------------------------------------- */ -#include "fluid_diffusion_model_inline_impl.hh" - -#endif /* __AKANTU_FLUID_DIFFUSION_MODEL_HH__ */ diff --git a/src/model/fluid_diffusion/fluid_diffusion_model_event_handler.hh b/src/model/fluid_diffusion/fluid_diffusion_model_event_handler.hh deleted file mode 100644 index 52914f7d8..000000000 --- a/src/model/fluid_diffusion/fluid_diffusion_model_event_handler.hh +++ /dev/null @@ -1,56 +0,0 @@ -/** - * @file fluid_diffusion_model_event_handler.hh - * - * @author Emil Gallyamov - * - * @date creation: Mon Aug 19 2019 - * @date last modification: - * - * @brief EventHandler implementation for FluidDiffusionEvents - * modification of the SolidMecanicsModelEventHandler - * - * @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 . - * - */ - -/* -------------------------------------------------------------------------- */ - -#ifndef __AKANTU_FLUID_DIFFUSION_MODEL_EVENT_HANDLER_HH__ -#define __AKANTU_FLUID_DIFFUSION_MODEL_EVENT_HANDLER_HH__ - -namespace akantu { - -/// akantu::FluidDiffusionModelEvent is the base event for model -namespace FluidDiffusionModelEvent {} - -/// akantu::SolidMechanicsModelEvent -class FluidDiffusionModelEventHandler { - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - virtual ~FluidDiffusionModelEventHandler() = default; - - template friend class EventHandlerManager; - -}; - -} // namespace akantu - -#endif /* __AKANTU_FLUID_DIFFUSION_MODEL_EVENT_HANDLER_HH__ */ diff --git a/src/model/fluid_diffusion/fluid_diffusion_model_inline_impl.hh b/src/model/fluid_diffusion/fluid_diffusion_model_inline_impl.hh deleted file mode 100644 index ea8a89cd5..000000000 --- a/src/model/fluid_diffusion/fluid_diffusion_model_inline_impl.hh +++ /dev/null @@ -1,42 +0,0 @@ -/** - * @file heat_transfer_model_inline_impl.cc - * - * @author Guillaume Anciaux - * @author Srinivasa Babu Ramisetti - * @author Nicolas Richart - * - * @date creation: Fri Aug 20 2010 - * @date last modification: Wed Jan 31 2018 - * - * @brief Implementation of the inline functions of the HeatTransferModel class - * - * @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 . - * - */ - -#ifndef __AKANTU_FLUID_DIFFUSION_MODEL_INLINE_IMPL_CC__ -#define __AKANTU_FLUID_DIFFUSION_MODEL_INLINE_IMPL_CC__ - -namespace akantu { - -/* -------------------------------------------------------------------------- */ - -} // namespace akantu - -#endif /* __AKANTU_FLUID_DIFFUSION_MODEL_INLINE_IMPL_CC__ */ diff --git a/src/model/model.hh b/src/model/model.hh index bdec4baeb..238f4a5a7 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,365 +1,359 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Interface of a model * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_HH_ #define AKANTU_MODEL_HH_ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Model(Mesh & mesh, const ModelType & type, UInt dim = _all_dimensions, const ID & id = "model"); ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif -#ifdef AKANTU_FLUID_DIFFUSION - case ModelType::_fluid_diffusion_model: - this->initFullImpl(FluidDiffusionModelOptions{ - use_named_args, std::forward(_pack)...}); - break; -#endif #ifdef AKANTU_PHASE_FIELD case ModelType::_phase_field_model: this->initFullImpl(PhaseFieldModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_EMBEDDED case ModelType::_embedded_model: this->initFullImpl(EmbeddedInterfaceModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel() = 0; virtual void initFEEngineBoundary(); /// function to print the containt of the class void printself(std::ostream & /*stream*/, int /*indent*/ = 0) const override{}; public: /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_TO_IMPLEMENT(); } protected: template void allocNodalField(std::unique_ptr> & array, UInt nb_component, const ID & name) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID &) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /* ------------------------------------------------------------------------ */ /* Pack and unpack hexlper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & field_id); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual std::shared_ptr createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldUInt(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createElementalField( const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind /*kind*/) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); /* ------------------------------------------------------------------------ */ virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream & /*stream*/, const Model & /*_this*/); ID id; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// parser to the pointer to use Parser & parser; /// default ElementKind for dumper ElementKind dumper_default_element_kind{_ek_regular}; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.hh" #endif /* AKANTU_MODEL_HH_ */ diff --git a/src/model/model_options.hh b/src/model/model_options.hh index 1007ca3e1..a1fcebf60 100644 --- a/src/model/model_options.hh +++ b/src/model/model_options.hh @@ -1,168 +1,154 @@ /** * @file model_options.hh * * @author Nicolas Richart * * @date creation: Mon Dec 04 2017 * @date last modification: Wed Jan 31 2018 * * @brief A Documented file. * * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_named_argument.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_OPTIONS_HH_ #define AKANTU_MODEL_OPTIONS_HH_ namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(analysis_method); } struct ModelOptions { explicit ModelOptions(AnalysisMethod analysis_method = _static) : analysis_method(analysis_method) {} template ModelOptions(use_named_args_t /*unused*/, pack &&... _pack) : ModelOptions(OPTIONAL_NAMED_ARG(analysis_method, _static)) {} virtual ~ModelOptions() = default; AnalysisMethod analysis_method; }; #ifdef AKANTU_SOLID_MECHANICS /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelOptions : public ModelOptions { explicit SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template SolidMechanicsModelOptions(use_named_args_t /*unused*/, pack &&... _pack) : SolidMechanicsModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif /* -------------------------------------------------------------------------- */ #ifdef AKANTU_COHESIVE_ELEMENT namespace { DECLARE_NAMED_ARGUMENT(is_extrinsic); } /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions { SolidMechanicsModelCohesiveOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool extrinsic = false) : SolidMechanicsModelOptions(analysis_method), is_extrinsic(extrinsic) {} template SolidMechanicsModelCohesiveOptions(use_named_args_t /*unused*/, pack &&... _pack) : SolidMechanicsModelCohesiveOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(is_extrinsic, false)) {} bool is_extrinsic{false}; }; #endif #ifdef AKANTU_HEAT_TRANSFER /* -------------------------------------------------------------------------- */ struct HeatTransferModelOptions : public ModelOptions { explicit HeatTransferModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template HeatTransferModelOptions(use_named_args_t /*unused*/, pack &&... _pack) : HeatTransferModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif -#ifdef AKANTU_FLUID_DIFFUSION -/* -------------------------------------------------------------------------- - */ -struct FluidDiffusionModelOptions : public ModelOptions { - explicit FluidDiffusionModelOptions(AnalysisMethod analysis_method = _static) - : ModelOptions(analysis_method) {} - - template - FluidDiffusionModelOptions(use_named_args_t, pack &&... _pack) - : FluidDiffusionModelOptions( - OPTIONAL_NAMED_ARG(analysis_method, _static)) {} -}; -#endif - #ifdef AKANTU_PHASE_FIELD /* -------------------------------------------------------------------------- */ struct PhaseFieldModelOptions : public ModelOptions { explicit PhaseFieldModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass) : ModelOptions(analysis_method) {} template PhaseFieldModelOptions(use_named_args_t, pack &&... _pack) : PhaseFieldModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} }; #endif /* -------------------------------------------------------------------------- */ #ifdef AKANTU_EMBEDDED namespace { DECLARE_NAMED_ARGUMENT(init_intersections); } /* -------------------------------------------------------------------------- */ struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions { /** * @brief Constructor for EmbeddedInterfaceModelOptions * @param analysis_method see SolidMechanicsModelOptions * @param init_intersections compute intersections */ EmbeddedInterfaceModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool init_intersections = true) : SolidMechanicsModelOptions(analysis_method), has_intersections(init_intersections) {} template EmbeddedInterfaceModelOptions(use_named_args_t /*unused*/, pack &&... _pack) : EmbeddedInterfaceModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(init_intersections, true)) {} /// Should consider reinforcements bool has_intersections; }; #endif } // namespace akantu #endif /* AKANTU_MODEL_OPTIONS_HH_ */