VolumeElementForcesAndSourcesCoreOnSide.hpp

Go to the documentation of this file.

/** \file VolumeElementForcesAndSourcesCoreOnSide.hpp
  \brief Side volume element for skeleton integration.
 
  This header defines VolumeElementForcesAndSourcesCoreOnSide, which extends
  VolumeElementForcesAndSourcesCore to evaluate volume finite elements on a
  selected face (the "skeleton"). It provides connectivity maps, face
  orientation information, and access to face geometry (normals, tangents)
  for implementing mortar methods, flux terms, and interface conditions.
 
  The nested UserDataOperator adapts the standard volume operator interface
  for side evaluations. Users inherit from UserDataOperator and override
  methods like doWork() to assemble residuals and tangent contributions using
  face normals at Gauss points, face coordinates, and the skeleton sense.
 
  Typical usage:
  1. Instantiate VolumeElementForcesAndSourcesCoreOnSide for each side
  2. Attach a derived UserDataOperator to handle side-specific assembly
  3. Loop over face elements on the skeleton
  4. For each face, find adjacent volume elements (sides) and loop over both
  5. Access getFaceSideNumber(), getSkeletonSense(), getNormal(), etc. in your
     operator to retrieve face geometry and orientation for each side
 
  \author Anonymous author contributing to MoFEM under MiT License
  \date June 2024
*/
 
#ifndef __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
#define __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
 
using namespace boost::numeric;
 
namespace MoFEM {
 
/**
 * \brief Base volume element used to integrate on skeleton
 * \ingroup mofem_forces_and_sources_volume_element
 */
struct VolumeElementForcesAndSourcesCoreOnSide
    : public VolumeElementForcesAndSourcesCore {
 
  using VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore;
 
  /**
   * @brief Get the face nodes mapped on the volume element
   *
   * \todo This is not general; e.g., for quads the number of nodes is 4.
   *
   * @return const std::array<int, 4>&
   */
  inline const std::array<int, 4> &getFaceConnMap() const;
 
  /**
   * @brief Get face nodes mapped on volume
   *
   * \todo This is not general; e.g., for prisms or hexes the size of the fixed
   * array is incorrect.
   *
   * @return const std::array<int, 8>&
   */
  inline const std::array<int, 8> &getTetConnMap() const;
 
  /**
   * @brief Get node on volume opposite to volume element
   *
   * \todo This is not general; for prisms or hexes the opposite node is
   * not unique.
   *
   * @return int
   */
  inline int getOppositeNode() const;
 
  /**
   * @brief Sense face on volume
   * @deprecated use getSkeletonSense()
   *
   * @return int
   */
  DEPRECATED inline int getFaceSense() const;
 
  /**
   * @brief Sense face on volume
   *
   * @return int
   */
  inline int getSkeletonSense() const;
 
  /**
   * @brief Face number on the volume
   *
   * @return int
   */
  inline int getFaceSideNumber() const;
 
  /** \brief default operator for TET element
   * \ingroup mofem_forces_and_sources_volume_element
   */
  struct UserDataOperator;
 
  /**
   * \brief Compute quadrature rule index for a given polynomial order
   * 
   * Note: This function for volume side by default return -1, since the integartion
   * rule is set by 
   * 
   *
   * @param order Requested polynomial order
   * @return Integration rule index compatible with the requested order
   */
  int getRule(int order);
 
  /**
   * \brief Configure Gauss points for the element side
   *
   * @param order Requested polynomial order for integration
   * @return Error code (0 on success)
   */
  MoFEMErrorCode setGaussPts(int order);
 
  /**
   * \brief Execute user operator   
   * @return Error code (0 on success)
   */
  MoFEMErrorCode operator()();
 
private:
 
  /**
   * \brief Implementation detail for operator()
   *
   * @return Error code (0 on success)
   */
  MoFEMErrorCode operatorImpl();
 
  int faceSense;      ///< Sense of face, could be 1 or -1
  int faceSideNumber; ///< Face side number
  int nbNodesOnFace;  ///< Number of nodes on face
  std::array<int, 4> faceConnMap;
  std::array<int, 8> tetConnMap;
  int oppositeNode;
  bool operatorImplCalled = false; ///< flag to check if operatorImpl was called
};
 
/** \brief default operator for TET element
 * \ingroup mofem_forces_and_sources_volume_element
 */
struct VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator
    : public VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
 
  /**
   * \brief Access the owning volume finite element
   *
   * Volume element is side element, to face element.
   *
   * @return Pointer to the parent volume finite element
   */
  inline VolumeElementForcesAndSourcesCoreOnSide *getVolumeFE() const;
 
  /**
   * \brief Access the face finite element
   * Face element, execute side element, that is volume, on which the
   * volume side operators are executed.
   *
   * @return Pointer to the face finite element attached to this side
   */
  inline FaceElementForcesAndSourcesCore *getFaceFE() const;
 
  /**
   * \brief Get face sense with respect to the volume
   * \deprecated Use getSkeletonSense()
   * @return Face sense (+1 or -1)
   */
  DEPRECATED inline int getFaceSense() const;
 
  /**
   * @brief Get the skeleton sense
   *
   * Calls getFaceSense().
   *
   * @return Skeleton sense (+1 or -1)
   */
  inline int getSkeletonSense() const;
 
  /**
   * \brief Get face side number with respect to the volume
   * @return Face side number
   */
  inline int getFaceSideNumber() const;
 
  /**
   * \brief Query if an edge belongs to the current face
   * @param ee Local edge index (0-5 for tetrahedra)
   * @return True when the edge lies on the active face
   */
  inline bool getEdgeFace(const int ee) const;
 
  /**
   * \brief Get the face normal for this side element
   * @return Face normal vector
   */
  inline VectorDouble &getNormal();
 
  /** \brief Get first triangle tangent
   */
  inline VectorDouble &getTangent1();
 
  /** \brief Get second triangle tangent
   */
  inline VectorDouble &getTangent2();
 
  /** \brief Get normal as tensor
   */
  inline auto getFTensor1Normal();
 
  /** \brief Get tangent one as tensor
   */
  inline auto getFTensor1Tangent1();
 
  /** \brief Get tangent two as tensor
   */
  inline auto getFTensor1Tangent2();
 
  /** \brief If higher-order geometry is available, return normals at Gauss points.
   *
   * \note The returned matrix has zero rows and columns when no higher-order
   * approximation of geometry is available.
   */
  inline MatrixDouble &getNormalsAtGaussPts();
 
  /** \brief If higher-order geometry is available, return normals at Gauss points.
   *
   * \param gg Gauss point number
   */
  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
 
  /** \brief Get normal at integration points
   *
   * Example:
   * \code
   * double nrm2;
   * FTensor::Index<'i',3> i;
   * auto t_normal = getFTensor1NormalsAtGaussPts();
   * for(int gg = gg!=data.getN().size1();gg++) {
   *   nrm2 = sqrt(t_normal(i)*t_normal(i));
   *   ++t_normal;
   * }
   * \endcode
   */
  inline auto getFTensor1NormalsAtGaussPts();
 
  /** \brief Get face coordinates at Gauss points.
   *
   * \note Coordinates should match the values returned by getCoordsAtGaussPts
   * on tetrahedra. If they differ, it indicates an error or an unusual setup.
   */
  inline MatrixDouble &getFaceCoordsAtGaussPts();
 
protected:
  /**
   * \brief Assign base class pointer for the operator
   * @param ptr Pointer to the forces and sources core instance
   * @return Error code (0 on success)
   */
  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
};
 
const std::array<int, 4> &
VolumeElementForcesAndSourcesCoreOnSide::getFaceConnMap() const {
  return faceConnMap;
}
 
const std::array<int, 8> &
VolumeElementForcesAndSourcesCoreOnSide::getTetConnMap() const {
  return tetConnMap;
}
 
int VolumeElementForcesAndSourcesCoreOnSide::getOppositeNode() const {
  return oppositeNode;
}
 
int VolumeElementForcesAndSourcesCoreOnSide::getFaceSense() const {
  return getSkeletonSense();
}
 
int VolumeElementForcesAndSourcesCoreOnSide::getSkeletonSense() const {
  return faceSense;
}
 
int VolumeElementForcesAndSourcesCoreOnSide::getFaceSideNumber() const {
  return faceSideNumber;
}
 
FaceElementForcesAndSourcesCore *
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getFaceFE() const {
  return static_cast<FaceElementForcesAndSourcesCore *>(
      getVolumeFE()->sidePtrFE);
}
 
VectorDouble &
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getNormal() {
  return getFaceFE()->nOrmal;
}
 
VectorDouble &
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getTangent1() {
  return getFaceFE()->tangentOne;
}
 
VectorDouble &
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getTangent2() {
  return getFaceFE()->tangentTwo;
}
 
VolumeElementForcesAndSourcesCoreOnSide *
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getVolumeFE() const {
  return static_cast<VolumeElementForcesAndSourcesCoreOnSide *>(ptrFE);
}
 
int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getFaceSense()
    const {
  return getSkeletonSense();
}
 
int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getSkeletonSense() const {
  return getVolumeFE()->faceSense;
}
 
auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFTensor1Normal() {
  double *ptr = &*getNormal().data().begin();
  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
}
 
auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFTensor1Tangent1() {
  double *ptr = &*getTangent1().data().begin();
  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
}
 
auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFTensor1Tangent2() {
  double *ptr = &*getTangent2().data().begin();
  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
}
 
inline auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFTensor1NormalsAtGaussPts() {
  double *ptr = &*getNormalsAtGaussPts().data().begin();
  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
                                                            &ptr[2]);
}
 
int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFaceSideNumber() const {
  return getVolumeFE()->faceSideNumber;
}
 
MatrixDouble &VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getNormalsAtGaussPts() {
  return getFaceFE()->normalsAtGaussPts;
}
 
MatrixDouble &VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
    getFaceCoordsAtGaussPts() {
  return getFaceFE()->coordsAtGaussPts;
}
 
bool VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getEdgeFace(
    const int ee) const {
  constexpr bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
                                         {false, true, false, true}, // e1
                                         {false, false, true, true}, // e2
                                         {true, false, true, false}, // e3
                                         {true, true, false, false}, // e4
                                         {false, true, true, false}};
  return edges_on_faces[ee][getFaceSideNumber()];
}
 
ublas::matrix_row<MatrixDouble>
VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getNormalsAtGaussPts(
    const int gg) {
  return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
}
 
/**
 * @deprecated Use VolumeElementForcesAndSourcesCore.
 */
DEPRECATED typedef VolumeElementForcesAndSourcesCoreOnSide
    VolumeElementForcesAndSourcesCoreOnSideBase;
 
/**
 * @deprecated Do not use this template; it is obsolete.
 */
template <int SWITCH>
struct VolumeElementForcesAndSourcesCoreOnSideSwitch
    : public VolumeElementForcesAndSourcesCoreOnSide {
  using VolumeElementForcesAndSourcesCoreOnSide::
      VolumeElementForcesAndSourcesCoreOnSide;
  using UserDataOperator =
      VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;
};
 
} // namespace MoFEM
 
#endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__