diff --git a/src/fe_engine/fe_engine.hh b/src/fe_engine/fe_engine.hh
index 24bc43b8c..ee6fa7f90 100644
--- a/src/fe_engine/fe_engine.hh
+++ b/src/fe_engine/fe_engine.hh
@@ -1,370 +1,363 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "element_type_map.hh"
#include "mesh_events.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_FE_ENGINE_HH_
#define AKANTU_FE_ENGINE_HH_
namespace akantu {
class Mesh;
class Integrator;
class ShapeFunctions;
class DOFManager;
class Element;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/**
* The generic FEEngine class derived in a FEEngineTemplate class
* containing the
* shape functions and the integration method
*/
class FEEngine : public MeshEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
FEEngine(Mesh & mesh, Int element_dimension = _all_dimensions,
const ID & id = "fem");
~FEEngine() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
virtual void initShapeFunctions(GhostType ghost_type = _not_ghost) = 0;
/// extract the nodal values and store them per element
template
static void
extractNodalToElementField(const Mesh & mesh, const Array & nodal_f,
Array & elemental_f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter);
/// filter a field
template
static void
filterElementalData(const Mesh & mesh, const Array & quad_f,
Array & filtered_f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
virtual void
integrate(const Array & f, Array & intf, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// integrate a scalar value f on all elements of type "type"
[[nodiscard]] virtual Real
integrate(const Array & f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
- /// integrate f for all integration points of type "type" but don't sum over
- /// all integration points
- virtual void integrateOnIntegrationPoints(
- const Array & f, Array & intf, Int nb_degree_of_freedom,
- ElementType type, GhostType ghost_type = _not_ghost,
- const Array & filter_elements = empty_filter) const = 0;
-
/// integrate one element scalar value on all elements of type "type"
[[nodiscard]] Real integrate(const Ref f,
const Element & element) const {
return integrate(f, element.type, element.element, element.ghost_type);
}
private:
[[nodiscard]] virtual Real
integrate(const Ref f, ElementType type, Idx index,
GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* compatibility with old FEEngine fashion */
/* ------------------------------------------------------------------------ */
public:
/// get the number of integration points
[[nodiscard]] virtual Int
getNbIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const = 0;
/// get the precomputed shapes
[[nodiscard]] const virtual Array &
getShapes(ElementType type, GhostType ghost_type = _not_ghost,
Idx id = 0) const = 0;
/// get the derivatives of shapes
[[nodiscard]] virtual const Array &
getShapesDerivatives(ElementType type, GhostType ghost_type = _not_ghost,
Idx id = 0) const = 0;
/// get integration points
[[nodiscard]] virtual const MatrixXr &
getIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
/// Compute the gradient nablauq on the integration points of an element type
/// from nodal values u
virtual void gradientOnIntegrationPoints(
const Array & u, Array & nablauq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Interpolate a nodal field u at the integration points of an element type
/// -> uq
virtual void interpolateOnIntegrationPoints(
const Array & u, Array & uq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Interpolate a nodal field u at the integration points of many element
/// types -> uq
virtual void interpolateOnIntegrationPoints(
const Array & u, ElementTypeMapArray & uq,
const ElementTypeMapArray * filter_elements = nullptr) const = 0;
/// pre multiplies a tensor by the shapes derivaties
virtual void
computeBtD(const Array & Ds, Array & BtDs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left and right multiplies a tensor by the shapes derivaties
virtual void
computeBtDB(const Array & Ds, Array & BtDBs, Int order_d,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left multiples a vector by the shape functions
virtual void
computeNtb(const Array & bs, Array & Ntbs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left and right multiplies a tensor by the shapes
virtual void
computeNtbN(const Array & bs, Array & NtbNs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Compute the interpolation point position in the global coordinates for
/// many element types
virtual void computeIntegrationPointsCoordinates(
ElementTypeMapArray & integration_points_coordinates,
const ElementTypeMapArray * filter_elements = nullptr) const = 0;
/// Compute the interpolation point position in the global coordinates for an
/// element type
virtual void computeIntegrationPointsCoordinates(
Array & integration_points_coordinates, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Build pre-computed matrices for interpolation of field form integration
/// points at other given positions (interpolation_points)
virtual void initElementalFieldInterpolationFromIntegrationPoints(
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & interpolation_points_coordinates_matrices,
ElementTypeMapArray & integration_points_coordinates_inv_matrices,
const ElementTypeMapArray * element_filter) const = 0;
/// interpolate field at given position (interpolation_points) from given
/// values of this field at integration points (field)
virtual void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & result, const GhostType ghost_type,
const ElementTypeMapArray * element_filter) const = 0;
/// Interpolate field at given position from given values of this field at
/// integration points (field)
/// using matrices precomputed with
/// initElementalFieldInterplationFromIntegrationPoints
virtual void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray &
interpolation_points_coordinates_matrices,
const ElementTypeMapArray &
integration_points_coordinates_inv_matrices,
ElementTypeMapArray & result, const GhostType ghost_type,
const ElementTypeMapArray * element_filter) const = 0;
/// interpolate on a phyiscal point inside an element
virtual void interpolate(const Ref real_coords,
const Ref nodal_values,
Ref interpolated,
const Element & element) const = 0;
/// compute the shape on a provided point
virtual void computeShapes(const Ref real_coords, Int elem,
ElementType type, Ref shapes,
GhostType ghost_type = _not_ghost) const = 0;
/// compute the shape derivatives on a provided point
virtual void
computeShapeDerivatives(const Ref real_coords, Int element,
ElementType type, Ref shape_derivatives,
GhostType ghost_type = _not_ghost) const = 0;
/// assembles the lumped version of @f[ \int N^t rho N @f]
virtual void assembleFieldLumped(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type = _not_ghost) const = 0;
/// assembles the matrix @f[ \int N^t rho N @f]
virtual void assembleFieldMatrix(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* Other methods */
/* ------------------------------------------------------------------------ */
/// pre-compute normals on integration points
virtual void
computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) = 0;
/// pre-compute normals on integration points
virtual void
computeNormalsOnIntegrationPoints(const Array & /*field*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// pre-compute normals on integration points
virtual void computeNormalsOnIntegrationPoints(
const Array & /*field*/, Array & /*normal*/,
ElementType /*type*/, GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
private:
/// initialise the class
void init();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
using ElementTypesIteratorHelper =
ElementTypeMapArray::ElementTypesIteratorHelper;
[[nodiscard]] ElementTypesIteratorHelper
elementTypes(Int dim = _all_dimensions, GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const;
/// get the dimension of the element handeled by this fe_engine object
AKANTU_GET_MACRO_AUTO(ElementDimension, element_dimension);
/// get the mesh contained in the fem object
AKANTU_GET_MACRO_AUTO(Mesh, mesh);
/// get the mesh contained in the fem object
AKANTU_GET_MACRO_NOT_CONST(Mesh, mesh, Mesh &);
/// get the in-radius of an element
template
[[nodiscard]] static inline Real
getElementInradius(const Eigen::MatrixBase & coord,
ElementType type);
[[nodiscard]] inline Real getElementInradius(const Element & element) const;
/// get the normals on integration points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalsOnIntegrationPoints,
normals_on_integration_points, Real);
/// get cohesive element type for a given facet type
[[nodiscard]] static inline constexpr ElementType
getCohesiveElementType(ElementType type_facet);
/// get igfem element type for a given regular type
[[nodiscard]] static inline Vector
getIGFEMElementTypes(ElementType type);
/// get the interpolation element associated to an element type
[[nodiscard]] static inline constexpr auto
getInterpolationType(ElementType el_type);
/// get the shape function class (probably useless: see getShapeFunction in
/// fe_engine_template.hh)
[[nodiscard]] virtual const ShapeFunctions &
getShapeFunctionsInterface() const = 0;
/// get the integrator class (probably useless: see getIntegrator in
/// fe_engine_template.hh)
[[nodiscard]] virtual const Integrator & getIntegratorInterface() const = 0;
AKANTU_GET_MACRO(ID, id, ID);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
ID id;
/// spatial dimension of the problem
Int element_dimension;
/// the mesh on which all computation are made
Mesh & mesh;
/// normals at integration points
ElementTypeMapArray normals_on_integration_points;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const FEEngine & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "fe_engine_inline_impl.hh"
#include "fe_engine_template.hh"
#endif /* AKANTU_FE_ENGINE_HH_ */
diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index b53502f77..a442e0ee7 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,481 +1,474 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "fe_engine.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_
#define AKANTU_FE_ENGINE_TEMPLATE_HH_
namespace akantu {
class Integrator;
class ShapeFunctions;
} // namespace akantu
namespace akantu {
class DOFManager;
namespace fe_engine {
namespace details {
template struct AssembleLumpedTemplateHelper;
template struct AssembleFieldMatrixHelper;
} // namespace details
} // namespace fe_engine
template struct AssembleFieldMatrixStructHelper;
struct DefaultIntegrationOrderFunctor {
template static inline constexpr int getOrder() {
return ElementClassProperty::polynomial_degree;
}
};
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind = _ek_regular,
class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor>
class FEEngineTemplate : public FEEngine {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using Integ = I;
using Shape = S;
FEEngineTemplate(Mesh & mesh, Int spatial_dimension = _all_dimensions,
const ID & id = "fem", bool do_not_precompute = false);
~FEEngineTemplate() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
void initShapeFunctions(GhostType ghost_type = _not_ghost) override;
void initShapeFunctions(const Array & nodes,
GhostType ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
void
integrate(const Array & f, Array & intf, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// integrate a scalar value on all elements of type "type"
[[nodiscard]] Real
integrate(const Array & f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// integrate one element scalar value on all elements of type "type"
[[nodiscard]] Real
integrate(const Ref f, ElementType type, Int index,
GhostType ghost_type = _not_ghost) const override;
- /// integrate partially around an integration point (@f$ intf_q = f_q * J_q *
- /// w_q @f$)
- void integrateOnIntegrationPoints(
- const Array & f, Array & intf, Int nb_degree_of_freedom,
- ElementType type, GhostType ghost_type = _not_ghost,
- const Array & filter_elements = empty_filter) const override;
-
private:
template ::value and
kind_ == _ek_regular> * = nullptr>
inline void interpolateImpl(const Eigen::MatrixBase & real_coords,
const Eigen::MatrixBase & nodal_values,
Eigen::MatrixBase & interpolated,
const Element & element) const;
template ::value and
kind_ != _ek_regular> * = nullptr>
inline void interpolateImpl(const Eigen::MatrixBase & /*real_coords*/,
const Eigen::MatrixBase & /*nodal_values*/,
Eigen::MatrixBase & /*interpolated*/,
const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
public:
/// interpolate on a phyiscal point inside an element
void interpolate(const Ref real_coords,
const Ref nodal_values,
Ref interpolated,
const Element & element) const override;
/// get the number of integration points
[[nodiscard]] Int
getNbIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const override;
/// get shapes precomputed
[[nodiscard]] const Array & getShapes(ElementType type,
GhostType ghost_type = _not_ghost,
Int id = 0) const override;
/// get the derivatives of shapes
[[nodiscard]] const Array &
getShapesDerivatives(ElementType type, GhostType ghost_type = _not_ghost,
Int id = 0) const override;
/// get integration points
[[nodiscard]] inline const Matrix &
getIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const override;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
/// compute the gradient of a nodal field on the integration points
void gradientOnIntegrationPoints(
const Array & u, Array & nablauq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate a nodal field on the integration points
void interpolateOnIntegrationPoints(
const Array & u, Array & uq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate a nodal field on the integration points given a
/// by_element_type
void interpolateOnIntegrationPoints(
const Array & u, ElementTypeMapArray & uq,
const ElementTypeMapArray * filter_elements =
nullptr) const override;
/// pre multiplies a tensor by the shapes derivaties
void
computeBtD(const Array & Ds, Array & BtDs, ElementType type,
GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// left and right multiplies a tensor by the shapes derivaties
void
computeBtDB(const Array & Ds, Array & BtDBs, Int order_d,
ElementType type, GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// left multiples a vector by the shape functions
void computeNtb(const Array & bs, Array & Ntbs, ElementType type,
GhostType ghost_type,
const Array & filter_elements) const override;
/// left and right multiplies a tensor by the shapes
void
computeNtbN(const Array & bs, Array & NtbNs, ElementType type,
GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// compute the position of integration points given by an element_type_map
/// from nodes position
inline void computeIntegrationPointsCoordinates(
ElementTypeMapArray & quadrature_points_coordinates,
const ElementTypeMapArray * filter_elements =
nullptr) const override;
/// compute the position of integration points from nodes position
inline void computeIntegrationPointsCoordinates(
Array & quadrature_points_coordinates, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate field at given position (interpolation_points) from given
/// values of this field at integration points (field)
inline void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & result, GhostType ghost_type,
const ElementTypeMapArray * element_filter) const override;
/// Interpolate field at given position from given values of this field at
/// integration points (field)
/// using matrices precomputed with
/// initElementalFieldInterplationFromIntegrationPoints
inline void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray &
interpolation_points_coordinates_matrices,
const ElementTypeMapArray & quad_points_coordinates_inv_matrices,
ElementTypeMapArray & result, GhostType ghost_type,
const ElementTypeMapArray * element_filter) const override;
/// Build pre-computed matrices for interpolation of field form integration
/// points at other given positions (interpolation_points)
inline void initElementalFieldInterpolationFromIntegrationPoints(
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & interpolation_points_coordinates_matrices,
ElementTypeMapArray & quad_points_coordinates_inv_matrices,
const ElementTypeMapArray * element_filter = nullptr) const override;
/// find natural coords from real coords provided an element
void inverseMap(const Ref real_coords, Int element,
ElementType type, Ref natural_coords,
GhostType ghost_type = _not_ghost) const;
/// return true if the coordinates provided are inside the element, false
/// otherwise
inline bool contains(const Vector & real_coords, Int element,
ElementType type,
GhostType ghost_type = _not_ghost) const;
private:
template ::value and
kind_ != _ek_cohesive> * = nullptr>
inline void computeShapesImpl(const Eigen::MatrixBase & real_coords,
Int element, ElementType type,
Eigen::MatrixBase & shapes,
GhostType ghost_type = _not_ghost) const;
template ::value and
kind_ == _ek_cohesive> * = nullptr>
inline void computeShapesImpl(const Eigen::MatrixBase & /*real_coords*/,
Int /*element*/, ElementType /*type*/,
Eigen::MatrixBase & /*shapes*/,
GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
template and kind_ != _ek_cohesive> * =
nullptr>
inline void
computeShapeDerivativesImpl(const Eigen::MatrixBase & real_coords,
Int element, ElementType type,
Eigen::MatrixBase & shape_derivatives,
GhostType ghost_type = _not_ghost) const;
template and kind_ == _ek_cohesive> * =
nullptr>
inline void
computeShapeDerivativesImpl(const Eigen::MatrixBase & /*real_coords*/,
Int /*element*/, ElementType /*type*/,
Eigen::MatrixBase & /*shape_derivatives*/,
GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
public:
/// compute the shape on a provided point
inline void computeShapes(const Ref real_coords, Int element,
ElementType type, Ref shapes,
GhostType ghost_type = _not_ghost) const override {
this->template computeShapesImpl(real_coords, element, type, shapes,
ghost_type);
}
/// compute the shape derivatives on a provided point
inline void
computeShapeDerivatives(const Ref real_coords, Int element,
ElementType type, Ref shape_derivatives,
GhostType ghost_type = _not_ghost) const override {
this->template computeShapeDerivativesImpl(
real_coords, element, type, shape_derivatives, ghost_type);
}
/* ------------------------------------------------------------------------ */
/* Other methods */
/* ------------------------------------------------------------------------ */
/// pre-compute normals on integration points
void
computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) override;
void
computeNormalsOnIntegrationPoints(const Array & field,
GhostType ghost_type = _not_ghost) override;
void computeNormalsOnIntegrationPoints(
const Array & field, Array & normal, ElementType type,
GhostType ghost_type = _not_ghost) const override;
template * = nullptr>
void computeNormalsOnIntegrationPoints(const Array & /*field*/,
Array & /*normal*/,
GhostType /*ghost_type*/) const {
AKANTU_TO_IMPLEMENT();
}
template <
ElementType type, ElementKind kind_ = kind,
std::enable_if_t * = nullptr>
void computeNormalsOnIntegrationPoints(const Array & field,
Array & normal,
GhostType ghost_type) const;
template <
ElementType type, ElementKind kind_ = kind,
std::enable_if_t * = nullptr>
void computeNormalsOnIntegrationPoints(const Array & field,
Array & normal,
GhostType ghost_type) const;
public:
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
void assembleFieldLumped(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type) const override;
private:
template * = nullptr>
void assembleFieldMatrixImpl(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type) const;
template * = nullptr>
void assembleFieldMatrixImpl(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type) const;
public:
/// assemble a field as a matrix (ex. rho to mass matrix)
void assembleFieldMatrix(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type) const override {
this->assembleFieldMatrixImpl(field_funct, matrix_id, dof_id, dof_manager,
type, ghost_type);
}
private:
friend struct fe_engine::details::AssembleLumpedTemplateHelper;
friend struct fe_engine::details::AssembleFieldMatrixHelper;
friend struct AssembleFieldMatrixStructHelper;
/// templated function to compute the scaling to assemble a lumped matrix
template
void assembleFieldLumped(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
GhostType ghost_type) const;
/// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j
/// dV = \int \rho \varphi_i dV @f$
template
void assembleLumpedRowSum(const Array & field, const ID & matrix_id,
const ID & dof_id, DOFManager & dof_manager,
GhostType ghost_type) const;
/// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
template
void assembleLumpedDiagonalScaling(const Array & field,
const ID & matrix_id, const ID & dof_id,
DOFManager & dof_manager,
GhostType ghost_type) const;
/// assemble a field as a matrix (ex. rho to mass matrix)
template
void assembleFieldMatrix(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
GhostType ghost_type) const;
#ifdef AKANTU_STRUCTURAL_MECHANICS
/// assemble a field as a matrix for structural elements (ex. rho to mass
/// matrix)
template
void assembleFieldMatrix(const Array & field_1,
Int nb_degree_of_freedom, SparseMatrix & M,
Array * n,
ElementTypeMapArray & rotation_mat,
GhostType ghost_type) const;
#endif
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler interface */
/* ------------------------------------------------------------------------ */
public:
void onElementsAdded(const Array & /*new_elements*/,
const NewElementsEvent & /*unused*/) override;
void onElementsRemoved(const Array & /*unused*/,
const ElementTypeMapArray & /*unused*/,
const RemovedElementsEvent & /*unused*/) override;
void onElementsChanged(const Array & /*unused*/,
const Array & /*unused*/,
const ElementTypeMapArray & /*unused*/,
const ChangedElementsEvent & /*unused*/) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the shape class (probably useless: see getShapeFunction)
const ShapeFunctions & getShapeFunctionsInterface() const override {
return shape_functions;
};
/// get the shape class
const Shape & getShapeFunctions() const { return shape_functions; };
/// get the integrator class (probably useless: see getIntegrator)
const Integrator & getIntegratorInterface() const override {
return integrator;
};
/// get the integrator class
const Integ & getIntegrator() const { return integrator; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
Integ integrator;
Shape shape_functions;
};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "fe_engine_template_tmpl.hh"
#include "fe_engine_template_tmpl_field.hh"
/* -------------------------------------------------------------------------- */
/* Shape Linked specialization */
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_STRUCTURAL_MECHANICS)
#include "fe_engine_template_tmpl_struct.hh"
#endif
/* -------------------------------------------------------------------------- */
/* Shape IGFEM specialization */
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_IGFEM)
#include "fe_engine_template_tmpl_igfem.hh"
#endif
#endif /* AKANTU_FE_ENGINE_TEMPLATE_HH_ */
diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh
index 925af4921..0af6d18a3 100644
--- a/src/fe_engine/fe_engine_template_tmpl.hh
+++ b/src/fe_engine/fe_engine_template_tmpl.hh
@@ -1,843 +1,806 @@
/**
* Copyright (©) 2011-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "dof_manager.hh"
#include "fe_engine_template.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
FEEngineTemplate::FEEngineTemplate(
Mesh & mesh, Int spatial_dimension, const ID & id, bool do_not_precompute)
: FEEngine(mesh, spatial_dimension, id),
integrator(mesh, spatial_dimension, id),
shape_functions(mesh, spatial_dimension, id) {
if (not do_not_precompute) {
initShapeFunctions(_not_ghost);
initShapeFunctions(_ghost);
}
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
FEEngineTemplate::~FEEngineTemplate() =
default;
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
void FEEngineTemplate::
gradientOnIntegrationPoints(const Array & u, Array & nablauq,
Int nb_degree_of_freedom, ElementType type,
GhostType ghost_type,
const Array & filter_elements) const {
AKANTU_DEBUG_IN();
auto nb_element = mesh.getNbElement(type, ghost_type);
if (filter_elements != empty_filter) {
nb_element = filter_elements.size();
}
auto nb_points =
shape_functions.getIntegrationPoints(type, ghost_type).cols();
#ifndef AKANTU_NDEBUG
auto element_dimension = mesh.getSpatialDimension(type);
AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
"The vector u(" << u.getID()
<< ") has not the good size.");
AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
"The vector u("
<< u.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(
nablauq.getNbComponent() == nb_degree_of_freedom * element_dimension,
"The vector nablauq(" << nablauq.getID()
<< ") has not the good number of component.");
#endif
nablauq.resize(nb_element * nb_points);
auto call = [&](auto && integral_type) {
constexpr ElementType type = std::decay_t::value;
if (element_dimension == ElementClass::getSpatialDimension())
shape_functions.template gradientOnIntegrationPoints(
u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
};
tuple_dispatch>(call, type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
void FEEngineTemplate::initShapeFunctions(
GhostType ghost_type) {
initShapeFunctions(mesh.getNodes(), ghost_type);
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
void FEEngineTemplate::initShapeFunctions(
const Array & nodes, GhostType ghost_type) {
AKANTU_DEBUG_IN();
for (const auto & type :
mesh.elementTypes(element_dimension, ghost_type, kind)) {
integrator.initIntegrator(nodes, type, ghost_type);
const auto & control_points = getIntegrationPoints(type, ghost_type);
shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
void FEEngineTemplate::integrate(
const Array & f, Array & intf, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type,
const Array & filter_elements) const {
auto nb_element = mesh.getNbElement(type, ghost_type);
if (filter_elements != empty_filter) {
nb_element = filter_elements.size();
}
#ifndef AKANTU_NDEBUG
auto nb_quadrature_points = getNbIntegrationPoints(type);
AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID() << " size " << f.size()
<< ") has not the good size ("
<< nb_element << ").");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
"The vector f("
<< f.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
"The vector intf("
<< intf.getID()
<< ") has not the good number of component.");
#endif
intf.resize(nb_element);
auto && call = [&](auto && enum_type) {
constexpr ElementType type = std ::decay_t::value;
integrator.template integrate(
f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
};
tuple_dispatch>(call, type);
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
Real FEEngineTemplate::integrate(
const Array & f, ElementType type, GhostType ghost_type,
const Array & filter_elements) const {
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
auto nb_element = mesh.getNbElement(type, ghost_type);
if (filter_elements != empty_filter) {
nb_element = filter_elements.size();
}
auto nb_quadrature_points = getNbIntegrationPoints(type, ghost_type);
AKANTU_DEBUG_ASSERT(
f.size() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID() << ") has not the good size. (" << f.size()
<< "!=" << nb_quadrature_points * nb_element << ")");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
"The vector f("
<< f.getID()
<< ") has not the good number of component.");
#endif
auto && call = [&](auto && enum_type) {
constexpr ElementType type = std ::decay_t::value;
return integrator.template integrate(f, ghost_type, filter_elements);
};
return tuple_dispatch>(call, type);
}
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind, class IntegrationOrderFunctor>
Real FEEngineTemplate::integrate(
const Ref f, ElementType type, Int index,
GhostType ghost_type) const {
auto && call = [&](auto && enum_type) {
constexpr ElementType type = std ::decay_t::value;
return integrator.template integrate(f, index, ghost_type);
};
return tuple_dispatch>(call, type);
}
-/* -------------------------------------------------------------------------- */
-template