diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
index 68dd33e7f..855158ebd 100644
--- a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
+++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
@@ -1,199 +1,195 @@
/**
* @file material_FE2.cc
*
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
*
* @brief Material for multi-scale simulations. It stores an
* underlying RVE on each integration point of the material.
*
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
*/
/* -------------------------------------------------------------------------- */
#include "material_FE2.hh"
#include "communicator.hh"
#include "solid_mechanics_model_RVE.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
MaterialFE2<spatial_dimension>::MaterialFE2(SolidMechanicsModel & model,
const ID & id)
: Parent(model, id), C("material_stiffness", *this) {
AKANTU_DEBUG_IN();
this->C.initialize(voigt_h::size * voigt_h::size);
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
MaterialFE2<spatial_dimension>::~MaterialFE2() = default;
/* -------------------------------------------------------------------------- */
template <UInt dim> void MaterialFE2<dim>::initialize() {
this->registerParam("element_type", el_type, _triangle_3,
_pat_parsable | _pat_modifiable,
"element type in RVE mesh");
this->registerParam("mesh_file", mesh_file, _pat_parsable | _pat_modifiable,
"the mesh file for the RVE");
this->registerParam("nb_gel_pockets", nb_gel_pockets,
_pat_parsable | _pat_modifiable,
"the number of gel pockets in each RVE");
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialFE2<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
Parent::initMaterial();
- /// create a Mesh and SolidMechanicsModel on each integration point of the
- /// material
- auto C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size);
-
- for (auto && data :
- enumerate(make_view(C(this->el_type), voigt_h::size, voigt_h::size))) {
- auto q = std::get<0>(data);
- auto & C = std::get<1>(data);
+ // create a Mesh and SolidMechanicsModel on each integration point of the
+ // material
+ auto & mesh = this->model.getMesh();
+ auto & g_ids = mesh.template getData<UInt>("global_ids", this->el_type);
+ auto const & element_filter = this->getElementFilter()(this->el_type);
+ for (auto && data : zip(element_filter)) {
+ UInt proc_el_id = std::get<0>(data);
+ UInt gl_el_id = g_ids(proc_el_id);
meshes.emplace_back(std::make_unique<Mesh>(
- spatial_dimension, "RVE_mesh_" + std::to_string(q), q + 1));
+ spatial_dimension, "RVE_mesh_" + std::to_string(gl_el_id)));
auto & mesh = *meshes.back();
mesh.read(mesh_file);
RVEs.emplace_back(std::make_unique<SolidMechanicsModelRVE>(
mesh, true, this->nb_gel_pockets, _all_dimensions,
- "SMM_RVE_" + std::to_string(q), q + 1));
+ "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);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialFE2<spatial_dimension>::computeStress(ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
// Compute thermal stresses first
Parent::computeStress(el_type, ghost_type);
Array<Real>::const_scalar_iterator sigma_th_it =
this->sigma_th(el_type, ghost_type).begin();
// Wikipedia convention:
// 2*eps_ij (i!=j) = voigt_eps_I
// http://en.wikipedia.org/wiki/Voigt_notation
Array<Real>::const_matrix_iterator C_it =
this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size);
// create vectors to store stress and strain in Voigt notation
// for efficient computation of stress
Vector<Real> voigt_strain(voigt_h::size);
Vector<Real> voigt_stress(voigt_h::size);
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
const Matrix<Real> & C_mat = *C_it;
const Real & sigma_th = *sigma_th_it;
/// copy strains in Voigt notation
for (UInt I = 0; I < voigt_h::size; ++I) {
/// copy stress in
Real voigt_factor = voigt_h::factors[I];
UInt i = voigt_h::vec[I][0];
UInt j = voigt_h::vec[I][1];
voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.;
}
// compute stresses in Voigt notation
voigt_stress.mul<false>(C_mat, voigt_strain);
/// copy stresses back in full vectorised notation
for (UInt I = 0; I < voigt_h::size; ++I) {
UInt i = voigt_h::vec[I][0];
UInt j = voigt_h::vec[I][1];
sigma(i, j) = sigma(j, i) = voigt_stress(I) + (i == j) * sigma_th;
}
++C_it;
++sigma_th_it;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialFE2<spatial_dimension>::computeTangentModuli(
- ElementType el_type, Array<Real> & tangent_matrix,
- GhostType ghost_type) {
+ ElementType el_type, Array<Real> & tangent_matrix, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<Real>::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 <UInt spatial_dimension>
void MaterialFE2<spatial_dimension>::advanceASR(
const Matrix<Real> & prestrain) {
AKANTU_DEBUG_IN();
for (auto && data :
zip(RVEs,
make_view(this->gradu(this->el_type), spatial_dimension,
spatial_dimension),
make_view(this->eigengradu(this->el_type), spatial_dimension,
spatial_dimension),
make_view(this->C(this->el_type), voigt_h::size, voigt_h::size),
this->delta_T(this->el_type))) {
auto & RVE = *(std::get<0>(data));
/// apply boundary conditions based on the current macroscopic displ.
/// gradient
RVE.applyBoundaryConditions(std::get<1>(data));
/// apply homogeneous temperature field to each RVE to obtain thermoelastic
/// effect
RVE.applyHomogeneousTemperature(std::get<4>(data));
/// advance the ASR in every RVE
RVE.advanceASR(prestrain);
/// compute the average eigen_grad_u
RVE.homogenizeEigenGradU(std::get<2>(data));
/// compute the new effective stiffness of the RVE
RVE.homogenizeStiffness(std::get<3>(data));
}
AKANTU_DEBUG_OUT();
}
INSTANTIATE_MATERIAL(material_FE2, MaterialFE2);
} // namespace akantu
diff --git a/packages/fluid_diffusion.cmake b/packages/fluid_diffusion.cmake
new file mode 100644
index 000000000..b64028200
--- /dev/null
+++ b/packages/fluid_diffusion.cmake
@@ -0,0 +1,41 @@
+#===============================================================================
+# @file heat_transfer.cmake
+#
+# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
+# @author Nicolas Richart <nicolas.richart@epfl.ch>
+#
+# @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 <http://www.gnu.org/licenses/>.
+#
+#===============================================================================
+
+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 ebd779e80..bfb83ae2a 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,677 +1,685 @@
/**
* @file aka_common.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <boost/preprocessor.hpp>
#include <limits>
#include <list>
#include <memory>
#include <string>
#include <type_traits>
#include <unordered_map>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Constants */
/* -------------------------------------------------------------------------- */
namespace {
[[gnu::unused]] constexpr UInt _all_dimensions{
std::numeric_limits<UInt>::max()};
#ifdef AKANTU_NDEBUG
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.};
#else
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{
std::numeric_limits<Real>::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) \
(htm_phi) \
(htm_gradient_phi) \
+ (fdm_pressure) \
+ (fdm_gradient_pressure) \
(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)
// 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 ---
+ _fdm_pressure, ///< synchronization of the nodal pressure
+ _fdm_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
+ _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
-
+ _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
};
#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<NodeFlag>;
return NodeFlag(under(a) & under(b));
}
inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) {
using under = std::underlying_type_t<NodeFlag>;
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<NodeFlag>;
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<GhostType_def>;
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<type> & 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 <typename T> std::string printMemorySize(UInt size);
/* -------------------------------------------------------------------------- */
struct TensorTrait {};
struct TensorProxyTrait {};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* Type traits */
/* -------------------------------------------------------------------------- */
namespace aka {
/* ------------------------------------------------------------------------ */
template <typename T> using is_tensor = std::is_base_of<akantu::TensorTrait, T>;
template <typename T>
using is_tensor_proxy = std::is_base_of<akantu::TensorProxyTrait, T>;
/* ------------------------------------------------------------------------ */
template <typename T> using is_scalar = std::is_arithmetic<T>;
/* ------------------------------------------------------------------------ */
template <typename R, typename T,
std::enable_if_t<std::is_reference<T>::value> * = nullptr>
bool is_of_type(T && t) {
return (
dynamic_cast<std::add_pointer_t<
std::conditional_t<std::is_const<std::remove_reference_t<T>>::value,
std::add_const_t<R>, R>>>(&t) != nullptr);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T> bool is_of_type(std::unique_ptr<T> & t) {
return (
dynamic_cast<std::add_pointer_t<
std::conditional_t<std::is_const<T>::value, std::add_const_t<R>, R>>>(
t.get()) != nullptr);
}
/* ------------------------------------------------------------------------ */
template <typename R, typename T,
std::enable_if_t<std::is_reference<T>::value> * = nullptr>
decltype(auto) as_type(T && t) {
static_assert(
disjunction<
std::is_base_of<std::decay_t<T>, std::decay_t<R>>, // down-cast
std::is_base_of<std::decay_t<R>, std::decay_t<T>> // up-cast
>::value,
"Type T and R are not valid for a as_type conversion");
return dynamic_cast<std::add_lvalue_reference_t<
std::conditional_t<std::is_const<std::remove_reference_t<T>>::value,
std::add_const_t<R>, R>>>(t);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T,
std::enable_if_t<std::is_pointer<T>::value> * = nullptr>
decltype(auto) as_type(T && t) {
return &as_type<R>(*t);
}
/* -------------------------------------------------------------------------- */
template <typename R, typename T>
decltype(auto) as_type(const std::shared_ptr<T> & t) {
return std::dynamic_pointer_cast<R>(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 <typename a, typename b> struct hash<std::pair<a, b>> {
hash() = default;
size_t operator()(const std::pair<a, b> & p) const {
size_t seed = ah(p.first);
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
const hash<a> ah{};
const hash<b> bh{};
};
} // namespace std
#endif // AKANTU_COMMON_HH_
diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in
index 275e536c9..8e55d4d55 100644
--- a/src/common/aka_config.hh.in
+++ b/src/common/aka_config.hh.in
@@ -1,93 +1,103 @@
/**
* @file aka_config.hh.in
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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)
+ * 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 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.
+ * 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 <http://www.gnu.org/licenses/>.
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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@;
+#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 23cd96a4b..fc74c828f 100644
--- a/src/fe_engine/shape_lagrange.hh
+++ b/src/fe_engine/shape_lagrange.hh
@@ -1,176 +1,179 @@
/**
* @file shape_lagrange.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "shape_lagrange_base.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SHAPE_LAGRANGE_HH_
#define AKANTU_SHAPE_LAGRANGE_HH_
namespace akantu {
/* -------------------------------------------------------------------------- */
template <class Shape> class ShapeCohesive;
class ShapeIGFEM;
template <ElementKind kind> 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<Real> & nodes,
const Matrix<Real> & integration_points,
- ElementType type,
- GhostType ghost_type);
+ ElementType type, GhostType ghost_type);
/// computes the shape functions derivatives for given interpolation points
template <ElementType type>
void computeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
Array<Real> & shape_derivatives, GhostType ghost_type,
const Array<UInt> & filter_elements = empty_filter) const;
void computeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
- Array<Real> & shape_derivatives, ElementType type,
- GhostType ghost_type,
+ Array<Real> & shape_derivatives, ElementType type, GhostType ghost_type,
const Array<UInt> & filter_elements) const override;
+ /// computes the shape functions variations for spatial dim != natural dim
+ template <ElementType type>
+ void computeShapeDerivativesOnIntegrationPoints1DIn2D(
+ const Array<Real> & nodes, const Matrix<Real> & integration_points,
+ Array<Real> & shape_derivatives, const GhostType & ghost_type,
+ const Array<UInt> & filter_elements = empty_filter) const;
+
/// pre compute all shapes on the element integration points from natural
/// coordinates
template <ElementType type>
void precomputeShapesOnIntegrationPoints(const Array<Real> & nodes,
GhostType ghost_type);
/// pre compute all shape derivatives on the element integration points from
/// natural coordinates
template <ElementType type>
- void
- precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes,
- GhostType ghost_type);
+ void precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes,
+ GhostType ghost_type);
/// interpolate nodal values on the integration points
template <ElementType type>
void interpolateOnIntegrationPoints(
const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
GhostType ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
template <ElementType type>
void interpolateOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
const Array<Real> & shapes, GhostType ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// interpolate on physical point
template <ElementType type>
void interpolate(const Vector<Real> & real_coords, UInt elem,
const Matrix<Real> & nodal_values,
- Vector<Real> & interpolated,
- GhostType ghost_type) const;
+ Vector<Real> & interpolated, GhostType ghost_type) const;
/// compute the gradient of u on the integration points
template <ElementType type>
void gradientOnIntegrationPoints(
const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
GhostType ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
template <ElementType type>
void computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
GhostType ghost_type,
const Array<UInt> & filter_elements) const;
template <ElementType type>
void computeBtDB(const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
GhostType ghost_type,
const Array<UInt> & filter_elements) const;
/// multiply a field by shape functions @f$ fts_{ij} = f_i * \varphi_j @f$
template <ElementType type>
void computeNtb(const Array<Real> & bs, Array<Real> & Ntbs,
GhostType ghost_type,
const Array<UInt> & filter_elements = empty_filter) const;
template <ElementType type>
void computeNtbN(const Array<Real> & bs, Array<Real> & NtbNs,
GhostType ghost_type,
const Array<UInt> & filter_elements) const;
/// find natural coords from real coords provided an element
template <ElementType type>
void inverseMap(const Vector<Real> & real_coords, UInt element,
Vector<Real> & natural_coords,
GhostType ghost_type = _not_ghost) const;
/// return true if the coordinates provided are inside the element, false
/// otherwise
template <ElementType type>
bool contains(const Vector<Real> & real_coords, UInt elem,
GhostType ghost_type) const;
/// compute the shape on a provided point
template <ElementType type>
void computeShapes(const Vector<Real> & real_coords, UInt elem,
Vector<Real> & shapes, GhostType ghost_type) const;
/// compute the shape derivatives on a provided point
template <ElementType type>
void computeShapeDerivatives(const Matrix<Real> & real_coords, UInt elem,
Tensor3<Real> & shapes,
GhostType ghost_type) const;
protected:
/// compute the shape derivatives on integration points for a given element
template <ElementType type>
inline void
computeShapeDerivativesOnCPointsByElement(const Matrix<Real> & node_coords,
const Matrix<Real> & natural_coords,
Tensor3<Real> & 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 0cd7ee9c4..6fa81e2fd 100644
--- a/src/fe_engine/shape_lagrange_inline_impl.hh
+++ b/src/fe_engine/shape_lagrange_inline_impl.hh
@@ -1,579 +1,651 @@
/**
* @file shape_lagrange_inline_impl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<type>(integration_points, ghost_type); \
precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type); \
if (ElementClass<type>::getNaturalSpaceDimension() == \
mesh.getSpatialDimension() || \
kind != _ek_regular) \
precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
template <ElementKind kind>
inline void ShapeLagrange<kind>::initShapeFunctions(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
ElementType type, GhostType ghost_type) {
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
}
#undef INIT_SHAPE_FUNCTIONS
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
inline void ShapeLagrange<kind>::computeShapeDerivativesOnCPointsByElement(
const Matrix<Real> & node_coords, const Matrix<Real> & natural_coords,
Tensor3<Real> & shapesd) const {
AKANTU_DEBUG_IN();
// compute dnds
Tensor3<Real> dnds(node_coords.rows(), node_coords.cols(),
natural_coords.cols());
ElementClass<type>::computeDNDS(natural_coords, dnds);
// compute jacobian
Tensor3<Real> J(node_coords.rows(), natural_coords.rows(),
natural_coords.cols());
ElementClass<type>::computeJMat(dnds, node_coords, J);
// compute dndx
ElementClass<type>::computeShapeDerivatives(J, dnds, shapesd);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::inverseMap(const Vector<Real> & real_coords,
UInt elem, Vector<Real> & natural_coords,
GhostType ghost_type) const {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt nb_nodes_per_element =
ElementClass<type>::getNbNodesPerInterpolationElement();
UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
Matrix<Real> 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<type>::inverseMap(real_coords, nodes_coord, natural_coords);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
bool ShapeLagrange<kind>::contains(const Vector<Real> & real_coords, UInt elem,
GhostType ghost_type) const {
UInt spatial_dimension = mesh.getSpatialDimension();
Vector<Real> natural_coords(spatial_dimension);
inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
return ElementClass<type>::contains(natural_coords);
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::interpolate(const Vector<Real> & real_coords,
UInt elem,
const Matrix<Real> & nodal_values,
Vector<Real> & interpolated,
GhostType ghost_type) const {
UInt nb_shapes = ElementClass<type>::getShapeSize();
Vector<Real> shapes(nb_shapes);
computeShapes<type>(real_coords, elem, shapes, ghost_type);
ElementClass<type>::interpolate(nodal_values, shapes, interpolated);
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeShapes(const Vector<Real> & real_coords,
UInt elem, Vector<Real> & shapes,
GhostType ghost_type) const {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
Vector<Real> natural_coords(spatial_dimension);
inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
ElementClass<type>::computeShapes(natural_coords, shapes);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeShapeDerivatives(
const Matrix<Real> & real_coords, UInt elem, Tensor3<Real> & 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<type>::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<Real> natural_coords(spatial_dimension, nb_points);
// Creates the matrix of natural coordinates
for (UInt i = 0; i < nb_points; i++) {
Vector<Real> real_point = real_coords(i);
Vector<Real> natural_point = natural_coords(i);
inverseMap<type>(real_point, elem, natural_point, ghost_type);
}
UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
Matrix<Real> 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<type>(nodes_coord, natural_coords,
shapesd);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
ShapeLagrange<kind>::ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
const ID & id)
: ShapeLagrangeBase(mesh, spatial_dimension, kind, id) {}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
Array<Real> & shape_derivatives, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt nb_nodes_per_element =
ElementClass<type>::getNbNodesPerInterpolationElement();
UInt nb_points = integration_points.cols();
UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
UInt size_of_shapesd = ElementClass<type>::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<Real> 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<Real>::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<Real> & X = *x_it;
Tensor3<Real> B(shapesd_val, spatial_dimension, nb_nodes_per_element,
nb_points);
computeShapeDerivativesOnCPointsByElement<type>(X, integration_points, B);
if (filter_elements == empty_filter) {
shapesd_val += size_of_shapesd * nb_points;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
Array<Real> & shape_derivatives, ElementType type, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
#define AKANTU_COMPUTE_SHAPES(type) \
computeShapeDerivativesOnIntegrationPoints<type>( \
nodes, integration_points, shape_derivatives, ghost_type, \
filter_elements);
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES);
#undef AKANTU_COMPUTE_SHAPES
}
+/* -------------------------------------------------------------------------- */
+template <ElementKind kind>
+template <ElementType type>
+void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints1DIn2D(
+ const Array<Real> & nodes, const Matrix<Real> & integration_points,
+ Array<Real> & shape_derivatives, const GhostType & ghost_type,
+ const Array<UInt> & filter_elements) const {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = mesh.getSpatialDimension();
+ UInt element_dimension = ElementClass<type>::getNaturalSpaceDimension();
+
+ UInt nb_nodes_per_element =
+ ElementClass<type>::getNbNodesPerInterpolationElement();
+
+ UInt nb_points = integration_points.cols();
+ UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
+
+ UInt size_of_shapesd = ElementClass<type>::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<Real> 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<Real>::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<Real> 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<Real> & x = std::get<0>(data);
+ auto & x_eq = std::get<1>(data);
+ for (auto node_nb : arange(1, nb_nodes_per_element)) {
+ Vector<Real> delta_x = Vector<Real>(x(node_nb)) - Vector<Real>(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<Real> B(shapesd_val, element_dimension, nb_nodes_per_element,
+ nb_points);
+ computeShapeDerivativesOnCPointsByElement<type>(X, integration_points, B);
+
+ if (filter_elements == empty_filter)
+ shapesd_val += size_of_shapesd * nb_points;
+ }
+
+ AKANTU_DEBUG_OUT();
+}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::precomputeShapesOnIntegrationPoints(
const Array<Real> & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
Matrix<Real> & natural_coords = integration_points(type, ghost_type);
UInt size_of_shapes = ElementClass<type>::getShapeSize();
Array<Real> & shapes_tmp =
shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
this->computeShapesOnIntegrationPoints<type>(nodes, natural_coords,
shapes_tmp, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::precomputeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
Matrix<Real> & natural_coords = integration_points(type, ghost_type);
UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
+ UInt spatial_dimension = mesh.getSpatialDimension();
+ UInt natural_dimension = ElementClass<type>::getNaturalSpaceDimension();
Array<Real> & shapes_derivatives_tmp =
shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
- this->computeShapeDerivativesOnIntegrationPoints<type>(
- nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
+ if (spatial_dimension != natural_dimension) {
+ this->computeShapeDerivativesOnIntegrationPoints1DIn2D<type>(
+ nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
+ } else {
+ this->computeShapeDerivativesOnIntegrationPoints<type>(
+ nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
+ }
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
const Array<Real> & shapes, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
UInt nb_nodes_per_element =
ElementClass<type>::getNbNodesPerInterpolationElement();
Array<Real> 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<type>(
u_el, out_uq, ghost_type, shapes, filter_elements);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
GhostType ghost_type, const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type),
"No shapes for the type "
<< shapes.printType(itp_type, ghost_type));
this->interpolateOnIntegrationPoints<type>(in_u, out_uq, nb_degree_of_freedom,
shapes(itp_type, ghost_type),
ghost_type, filter_elements);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::gradientOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_nablauq,
UInt nb_degree_of_freedom, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
InterpolationType itp_type = ElementClassProperty<type>::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<type>::getNbNodesPerInterpolationElement();
Array<Real> 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<type>(
u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type),
filter_elements);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeBtD(
const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
auto itp_type = ElementClassProperty<type>::interpolation_type;
const auto & shapes_derivatives =
this->shapes_derivatives(itp_type, ghost_type);
- auto spatial_dimension = mesh.getSpatialDimension();
+ auto spatial_dimension = ElementClass<type>::getNaturalSpaceDimension();
auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Array<Real> 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<false, false>(D, B);
}
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeBtDB(
const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
GhostType ghost_type, const Array<UInt> & filter_elements) const {
auto itp_type = ElementClassProperty<type>::interpolation_type;
const auto & shapes_derivatives =
this->shapes_derivatives(itp_type, ghost_type);
constexpr auto dim = ElementClass<type>::getSpatialDimension();
auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Array<Real> 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<dim>::size;
Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> 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<dim>::transferBMatrixToSymVoigtBMatrix(Bfull, B,
nb_nodes_per_element);
Bt_D.template mul<true, false>(B, D);
Bt_D_B.template mul<false, false>(Bt_D, B);
}
} else if (order_d == 2) {
Matrix<Real> 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<true, false>(B, D);
Bt_D_B.template mul<false, false>(Bt_D, B);
}
}
}
template <>
template <>
inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>(
const Array<Real> & /*Ds*/, Array<Real> & /*BtDBs*/, UInt /*order_d*/,
GhostType /*ghost_type*/, const Array<UInt> & /*filter_elements*/) const {
AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeNtbN(
const Array<Real> & bs, Array<Real> & NtbNs, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
auto itp_type = ElementClassProperty<type>::interpolation_type;
auto size_of_shapes = ElementClass<type>::getShapeSize();
auto nb_degree_of_freedom = bs.getNbComponent();
auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
Array<Real> 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<Real> 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<true, false>(N, b);
Nt_b_N.template mul<false, false>(Nt_b, N);
}
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeLagrange<kind>::computeNtb(
const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
Ntbs.resize(bs.size());
auto size_of_shapes = ElementClass<type>::getShapeSize();
auto itp_type = ElementClassProperty<type>::interpolation_type;
auto nb_degree_of_freedom = bs.getNbComponent();
Array<Real> 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<false, false>(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 1e8a45c53..35faa4efd 100644
--- a/src/mesh/element_type_map.cc
+++ b/src/mesh/element_type_map.cc
@@ -1,73 +1,69 @@
/**
* @file element_type_map.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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.getMesh().getSpatialDimension()
- : spatial_dimension,
- ghost_type, element_kind, true, false),
+ : MeshElementTypeMapArrayInitializer(fe_engine.getMesh(), nb_component,
+ spatial_dimension == UInt(-2)
+ ? fe_engine.getElementDimension()
+ : 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)
- : MeshElementTypeMapArrayInitializer(
- fe_engine.getMesh(), nb_component,
- spatial_dimension == UInt(-2)
- ? fe_engine.getMesh().getSpatialDimension()
- : spatial_dimension,
- ghost_type, element_kind, true, false),
+ UInt spatial_dimension, GhostType ghost_type, ElementKind element_kind)
+ : MeshElementTypeMapArrayInitializer(fe_engine.getMesh(), nb_component,
+ spatial_dimension == UInt(-2)
+ ? fe_engine.getElementDimension()
+ : 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 9fea35044..5a68bf6aa 100644
--- a/src/model/common/integration_scheme/integration_scheme.cc
+++ b/src/model/common/integration_scheme/integration_scheme.cc
@@ -1,92 +1,96 @@
/**
* @file integration_scheme.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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) {}
+ 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<Array<Real>>(
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 c6706df0a..f9c7eecd8 100644
--- a/src/model/common/integration_scheme/integration_scheme.hh
+++ b/src/model/common/integration_scheme/integration_scheme.hh
@@ -1,123 +1,125 @@
/**
* @file integration_scheme.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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<std::string> 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<std::unique_ptr<Array<Real>>> 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
new file mode 100644
index 000000000..c876c1f7c
--- /dev/null
+++ b/src/model/fluid_diffusion/fluid_diffusion_model.cc
@@ -0,0 +1,1273 @@
+/**
+ * @file fluid_diffusion_model.hh
+ *
+ * @author Emil Gallyamov <emil.gallyamov@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#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<Real> & 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<FEEngineType>(id + ":fem", mesh, spatial_dimension);
+
+#ifdef AKANTU_USE_IOHELPER
+ this->mesh.registerDumper<DumperParaview>("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<FEEngine>(getFEEngineClassBoundary<FEEngineType>(name));
+}
+/* -------------------------------------------------------------------------- */
+FluidDiffusionModel::~FluidDiffusionModel() = default;
+
+/* -------------------------------------------------------------------------- */
+void FluidDiffusionModel::assembleCapacityLumped(const GhostType & ghost_type) {
+ AKANTU_DEBUG_IN();
+
+ auto & fem = getFEEngineClass<FEEngineType>();
+ 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<Real> 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<ID, TimeStepSolverType>
+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<Array<Real>>(
+ 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<Array<Real>>(
+ 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<Real> 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<Real> 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<Real> 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<Real> 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<Real>::max();
+ Real max_aper_on_qpoint = std::numeric_limits<Real>::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<Real> 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<Real> & el_coord = std::get<0>(data);
+ Vector<Real> & 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<L_inf>());
+ }
+
+ 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<FEEngineType>();
+
+ 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<Real> & 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> & 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<akantu::Element>("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<const MaterialCohesive &>(material);
+ // const auto & normal_open_norm = mat.getNormalOpeningNorm(typecoh, gt);
+ const auto normal_open_norm_it =
+ make_view(material.getArray<Real>("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<Element>("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<Real>("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<UInt> & 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<UInt> 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<Real> 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<Element>("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<UInt> 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<UInt> 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<Real> 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<Real> 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<dumpers::Field> FluidDiffusionModel::createNodalFieldBool(
+ const std::string & field_name, const std::string & group_name,
+ __attribute__((unused)) bool padding_flag) {
+
+ std::map<std::string, Array<bool> *> 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<dumpers::Field> 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<std::string, Array<Real> *> 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<dumpers::Field> field =
+ mesh.createNodalField(real_nodal_fields[field_name], group_name);
+
+ return field;
+}
+
+/* -------------------------------------------------------------------------- */
+std::shared_ptr<dumpers::Field> FluidDiffusionModel::createElementalField(
+ const std::string & field_name, const std::string & group_name,
+ bool /*padding_flag*/, UInt /*spatial_dimension*/,
+ ElementKind element_kind) {
+
+ std::shared_ptr<dumpers::Field> field;
+
+ if (field_name == "partitions") {
+ field = mesh.createElementalField<UInt, dumpers::ElementPartitionField>(
+ mesh.getConnectivities(), group_name, this->spatial_dimension,
+ element_kind);
+ } else if (field_name == "pressure_gradient") {
+ ElementTypeMap<UInt> nb_data_per_elem =
+ this->mesh.getNbDataPerElem(pressure_gradient);
+
+ field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
+ pressure_gradient, group_name, this->spatial_dimension, element_kind,
+ nb_data_per_elem);
+ } else if (field_name == "permeability") {
+ ElementTypeMap<UInt> nb_data_per_elem =
+ this->mesh.getNbDataPerElem(permeability_on_qpoints);
+
+ field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
+ permeability_on_qpoints, group_name, this->spatial_dimension,
+ element_kind, nb_data_per_elem);
+ } else if (field_name == "aperture") {
+ ElementTypeMap<UInt> nb_data_per_elem =
+ this->mesh.getNbDataPerElem(aperture_on_qpoints);
+
+ field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
+ aperture_on_qpoints, group_name, this->spatial_dimension, element_kind,
+ nb_data_per_elem);
+ } else if (field_name == "pressure_on_qpoints") {
+ ElementTypeMap<UInt> nb_data_per_elem =
+ this->mesh.getNbDataPerElem(pressure_on_qpoints);
+
+ field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
+ pressure_on_qpoints, group_name, this->spatial_dimension, element_kind,
+ nb_data_per_elem);
+ }
+
+ return field;
+}
+
+/* -------------------------------------------------------------------------- */
+#else
+/* -------------------------------------------------------------------------- */
+std::shared_ptr<dumpers::Field> 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<dumpers::Field> 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<dumpers::Field> 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<UInt> & 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<UInt> & 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<UInt> & 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<Element> & elements,
+ const SynchronizationTag & tag) const {
+ AKANTU_DEBUG_IN();
+
+ UInt size = 0;
+ UInt nb_nodes_per_element = 0;
+ Array<Element>::const_iterator<Element> it = elements.begin();
+ Array<Element>::const_iterator<Element> 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<Element> & 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<Element> & 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<Real> & 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<Real> bary_facet(dim);
+
+ Real min_dist = std::numeric_limits<Real>::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
new file mode 100644
index 000000000..bcf576dbb
--- /dev/null
+++ b/src/model/fluid_diffusion/fluid_diffusion_model.hh
@@ -0,0 +1,381 @@
+/**
+ * @file fluid_diffusion_model.hh
+ *
+ * @author Emil Gallyamov <emil.gallyamov@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#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 <array>
+/* -------------------------------------------------------------------------- */
+
+#ifndef __AKANTU_FLUID_DIFFUSION_MODEL_HH__
+#define __AKANTU_FLUID_DIFFUSION_MODEL_HH__
+
+namespace akantu {
+template <ElementKind kind, class IntegrationOrderFunctor>
+class IntegratorGauss;
+template <ElementKind kind> class ShapeLagrange;
+} // namespace akantu
+
+namespace akantu {
+
+class FluidDiffusionModel
+ : public Model,
+ public DataAccessor<Element>,
+ public DataAccessor<UInt>,
+ public EventHandlerManager<FluidDiffusionModelEventHandler> {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+ using FEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
+
+ 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<ID, TimeStepSolverType>
+ 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<Real> 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<Real> & 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> & 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<Real> & 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<Element> & elements,
+ const SynchronizationTag & tag) const override;
+ inline void packData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ const SynchronizationTag & tag) const override;
+ inline void unpackData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ const SynchronizationTag & tag) override;
+
+ inline UInt getNbData(const Array<UInt> & indexes,
+ const SynchronizationTag & tag) const override;
+ inline void packData(CommunicationBuffer & buffer,
+ const Array<UInt> & indexes,
+ const SynchronizationTag & tag) const override;
+ inline void unpackData(CommunicationBuffer & buffer,
+ const Array<UInt> & indexes,
+ const SynchronizationTag & tag) override;
+
+ /* ------------------------------------------------------------------------ */
+ /* Dumpable interface */
+ /* ------------------------------------------------------------------------ */
+public:
+ std::shared_ptr<dumpers::Field>
+ createNodalFieldReal(const std::string & field_name,
+ const std::string & group_name,
+ bool padding_flag) override;
+
+ std::shared_ptr<dumpers::Field>
+ createNodalFieldBool(const std::string & field_name,
+ const std::string & group_name,
+ bool padding_flag) override;
+
+ std::shared_ptr<dumpers::Field>
+ 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<Real> &);
+ /// get the boundary vector
+ AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &);
+ /// get the external heat rate vector
+ AKANTU_GET_MACRO(ExternalFlux, *external_flux, Array<Real> &);
+ /// 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<Real> &);
+ /// get the temperature derivative
+ AKANTU_GET_MACRO(PressureRate, *pressure_rate, Array<Real> &);
+
+ 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<Array<Real>> pressure;
+
+ /// pressure derivatives array
+ std::unique_ptr<Array<Real>> pressure_rate;
+
+ // /// increment array (@f$\delta \dot P@f$ or @f$\delta P@f$)
+ // Array<Real> * increment{nullptr};
+
+ /// the speed of the changing temperature
+ ElementTypeMapArray<Real> pressure_gradient;
+
+ /// pressure field on quadrature points
+ ElementTypeMapArray<Real> pressure_on_qpoints;
+
+ /// pressure change on quadrature points
+ ElementTypeMapArray<Real> delta_pres_on_qpoints;
+
+ /// conductivity tensor on quadrature points
+ ElementTypeMapArray<Real> permeability_on_qpoints;
+
+ /// aperture on quadrature points
+ ElementTypeMapArray<Real> aperture_on_qpoints;
+
+ /// aperture on quadrature points
+ ElementTypeMapArray<Real> prev_aperture_on_qpoints;
+
+ /// vector k \grad T on quad points
+ ElementTypeMapArray<Real> k_gradp_on_qpoints;
+
+ /// external flux vector
+ std::unique_ptr<Array<Real>> external_flux;
+
+ /// residuals array
+ std::unique_ptr<Array<Real>> internal_flux;
+
+ /// boundary vector
+ std::unique_ptr<Array<bool>> 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<GhostType, bool> initial_permeability{{_not_ghost, true},
+ {_ghost, true}};
+ std::unordered_map<GhostType, UInt> 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
new file mode 100644
index 000000000..52914f7d8
--- /dev/null
+++ b/src/model/fluid_diffusion/fluid_diffusion_model_event_handler.hh
@@ -0,0 +1,56 @@
+/**
+ * @file fluid_diffusion_model_event_handler.hh
+ *
+ * @author Emil Gallyamov <emil.gallyamov@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+
+#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 <class EventHandler> 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
new file mode 100644
index 000000000..ea8a89cd5
--- /dev/null
+++ b/src/model/fluid_diffusion/fluid_diffusion_model_inline_impl.hh
@@ -0,0 +1,42 @@
+/**
+ * @file heat_transfer_model_inline_impl.cc
+ *
+ * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
+ * @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#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 ffcf1d6af..bdec4baeb 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,363 +1,365 @@
/**
* @file model.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <typeindex>
/* -------------------------------------------------------------------------- */
#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<std::string, std::unique_ptr<FEEngine>>;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
virtual void initFullImpl(const ModelOptions & options);
public:
template <typename... pack>
std::enable_if_t<are_named_argument<pack...>::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<decltype(_pack)>(_pack)...});
break;
#endif
#ifdef AKANTU_COHESIVE_ELEMENT
case ModelType::_solid_mechanics_model_cohesive:
this->initFullImpl(SolidMechanicsModelCohesiveOptions{
use_named_args, std::forward<decltype(_pack)>(_pack)...});
break;
#endif
#ifdef AKANTU_HEAT_TRANSFER
case ModelType::_heat_transfer_model:
this->initFullImpl(HeatTransferModelOptions{
use_named_args, std::forward<decltype(_pack)>(_pack)...});
break;
#endif
+#ifdef AKANTU_FLUID_DIFFUSION
+ case ModelType::_fluid_diffusion_model:
+ this->initFullImpl(FluidDiffusionModelOptions{
+ use_named_args, std::forward<decltype(_pack)>(_pack)...});
+ break;
+#endif
#ifdef AKANTU_PHASE_FIELD
case ModelType::_phase_field_model:
this->initFullImpl(PhaseFieldModelOptions{
use_named_args, std::forward<decltype(_pack)>(_pack)...});
break;
-#endif
+#endif
#ifdef AKANTU_EMBEDDED
case ModelType::_embedded_model:
this->initFullImpl(EmbeddedInterfaceModelOptions{
use_named_args, std::forward<decltype(_pack)>(_pack)...});
break;
#endif
default:
this->initFullImpl(ModelOptions{use_named_args,
std::forward<decltype(_pack)>(_pack)...});
}
}
template <typename... pack>
std::enable_if_t<not are_named_argument<pack...>::value>
initFull(pack &&... _pack) {
this->initFullImpl(std::forward<decltype(_pack)>(_pack)...);
}
/// initialize a new solver if needed
void initNewSolver(const AnalysisMethod & method);
protected:
/// get some default values for derived classes
virtual std::tuple<ID, TimeStepSolverType>
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 <typename T>
void allocNodalField(std::unique_ptr<Array<T>> & 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 <typename FEEngineClass>
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 <typename FEEngineClass>
inline FEEngineClass & getFEEngineClass(std::string name = "") const;
/// return the fem boundary object associated with a provided name
template <typename FEEngineClass>
inline FEEngineClass & getFEEngineClassBoundary(std::string name = "");
/// Get the type of analysis method used
AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod);
/* ------------------------------------------------------------------------ */
- /* Pack and unpack hexlper functions */
+ /* Pack and unpack hexlper functions */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbIntegrationPoints(const Array<Element> & 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<dumpers::Field> 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<dumpers::Field>
createNodalFieldReal(const std::string & /*field_name*/,
const std::string & /*group_name*/,
bool /*padding_flag*/) {
return nullptr;
}
virtual std::shared_ptr<dumpers::Field>
createNodalFieldUInt(const std::string & /*field_name*/,
const std::string & /*group_name*/,
bool /*padding_flag*/) {
return nullptr;
}
virtual std::shared_ptr<dumpers::Field>
createNodalFieldBool(const std::string & /*field_name*/,
const std::string & /*group_name*/,
bool /*padding_flag*/) {
return nullptr;
}
- virtual std::shared_ptr<dumpers::Field>
- createElementalField(const std::string & /*field_name*/,
- const std::string & /*group_name*/,
- bool /*padding_flag*/,
- UInt /*spatial_dimension*/,
- ElementKind /*kind*/) {
+ virtual std::shared_ptr<dumpers::Field> 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 a1fcebf60..1007ca3e1 100644
--- a/src/model/model_options.hh
+++ b/src/model/model_options.hh
@@ -1,154 +1,168 @@
/**
* @file model_options.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <typename... pack>
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 <typename... pack>
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 <typename... pack>
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 <typename... pack>
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 <typename... pack>
+ 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 <typename... pack>
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 <typename... pack>
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_ */