fe_map.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library 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 2.1 of the License, or (at your option) any later version.
 
// This library 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 this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_FE_MAP_H
#define LIBMESH_FE_MAP_H
 
// libMesh includes
#include "libmesh/reference_counted_object.h"
#include "libmesh/point.h"
#include "libmesh/vector_value.h"
#include "libmesh/fe_type.h"
 
// C++ includes
#include <memory>
 
namespace libMesh
{
 
// forward declarations
class Elem;
class Node;
 
class FEMap
{
public:
 
  FEMap(Real jtol = 0);
  virtual ~FEMap() = default;
 
  static std::unique_ptr<FEMap> build(FEType fe_type);
 
  static FEFamily map_fe_type(const Elem & elem);
 
  template<unsigned int Dim>
  void init_reference_to_physical_map(const std::vector<Point> & qp,
                                      const Elem * elem);
 
  // Non-templated version for runtime selection
  void init_reference_to_physical_map(unsigned int dim,
                                      const std::vector<Point> & qp,
                                      const Elem * elem);
 
  void compute_single_point_map(const unsigned int dim,
                                const std::vector<Real> & qw,
                                const Elem * elem,
                                unsigned int p,
                                const std::vector<const Node *> & elem_nodes,
                                bool compute_second_derivatives);
 
  virtual void compute_affine_map(const unsigned int dim,
                                  const std::vector<Real> & qw,
                                  const Elem * elem);
 
  virtual void compute_null_map(const unsigned int dim,
                                const std::vector<Real> & qw);
 
 
  virtual void compute_map(const unsigned int dim,
                           const std::vector<Real> & qw,
                           const Elem * elem,
                           bool calculate_d2phi);
 
  virtual void compute_face_map(int dim,
                                const std::vector<Real> & qw,
                                const Elem * side);
 
  void compute_edge_map(int dim,
                        const std::vector<Real> & qw,
                        const Elem * side);
 
  template<unsigned int Dim>
  void init_face_shape_functions(const std::vector<Point> & qp,
                                 const Elem * side);
 
  template<unsigned int Dim>
  void init_edge_shape_functions(const std::vector<Point> & qp,
                                 const Elem * edge);
 
  static Point map (const unsigned int dim,
                    const Elem * elem,
                    const Point & reference_point);
 
  static Point map_deriv (const unsigned int dim,
                          const Elem * elem,
                          const unsigned int j,
                          const Point & reference_point);
 
  static Point inverse_map (const unsigned int dim,
                            const Elem * elem,
                            const Point & p,
                            const Real tolerance = TOLERANCE,
                            const bool secure = true,
                            const bool extra_checks = true);
 
  static void inverse_map (unsigned int dim,
                           const Elem * elem,
                           const std::vector<Point> & physical_points,
                           std::vector<Point> &       reference_points,
                           const Real tolerance = TOLERANCE,
                           const bool secure = true,
                           const bool extra_checks = true);
 
  const std::vector<Point> & get_xyz() const
  { libmesh_assert(!calculations_started || calculate_xyz);
    calculate_xyz = true; return xyz; }
 
  const std::vector<Real> & get_jacobian() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return jac; }
 
  const std::vector<Real> & get_JxW() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return JxW; }
 
  const std::vector<RealGradient> & get_dxyzdxi() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxyzdxi_map; }
 
  const std::vector<RealGradient> & get_dxyzdeta() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxyzdeta_map; }
 
  const std::vector<RealGradient> & get_dxyzdzeta() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxyzdzeta_map; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  const std::vector<RealGradient> & get_d2xyzdxi2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdxi2_map; }
 
  const std::vector<RealGradient> & get_d2xyzdeta2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdeta2_map; }
 
  const std::vector<RealGradient> & get_d2xyzdzeta2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdzeta2_map; }
 
  const std::vector<RealGradient> & get_d2xyzdxideta() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdxideta_map; }
 
  const std::vector<RealGradient> & get_d2xyzdxidzeta() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdxidzeta_map; }
 
  const std::vector<RealGradient> & get_d2xyzdetadzeta() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xyzdetadzeta_map; }
 
#endif
 
  const std::vector<Real> & get_dxidx() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxidx_map; }
 
  const std::vector<Real> & get_dxidy() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxidy_map; }
 
  const std::vector<Real> & get_dxidz() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dxidz_map; }
 
  const std::vector<Real> & get_detadx() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return detadx_map; }
 
  const std::vector<Real> & get_detady() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return detady_map; }
 
  const std::vector<Real> & get_detadz() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return detadz_map; }
 
  const std::vector<Real> & get_dzetadx() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dzetadx_map; }
 
  const std::vector<Real> & get_dzetady() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dzetady_map; }
 
  const std::vector<Real> & get_dzetadz() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dzetadz_map; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  const std::vector<std::vector<Real>> & get_d2xidxyz2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2xidxyz2_map; }
 
  const std::vector<std::vector<Real>> & get_d2etadxyz2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2etadxyz2_map; }
 
  const std::vector<std::vector<Real>> & get_d2zetadxyz2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2zetadxyz2_map; }
#endif
 
  const std::vector<std::vector<Real>> & get_psi() const
  { return psi_map; }
 
  const std::vector<std::vector<Real>> & get_phi_map() const
  { libmesh_assert(!calculations_started || calculate_xyz);
    calculate_xyz = true; return phi_map; }
 
  const std::vector<std::vector<Real>> & get_dphidxi_map() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphidxi_map; }
 
  const std::vector<std::vector<Real>> & get_dphideta_map() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphideta_map; }
 
  const std::vector<std::vector<Real>> & get_dphidzeta_map() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphidzeta_map; }
 
  const std::vector<std::vector<Point>> & get_tangents() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return tangents; }
 
  const std::vector<Point> & get_normals() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return normals; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  const std::vector<Real> & get_curvatures() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return curvatures;}
 
#endif
 
  void print_JxW(std::ostream & os) const;
 
  void print_xyz(std::ostream & os) const;
 
  /* FIXME: PB: The public functions below break encapsulation! I'm not happy
     about it, but I don't have time to redo the infinite element routines.
     These are needed in inf_fe_boundary.C and inf_fe.C. A proper implementation
     would subclass this class and hide all the radial function business. */
  std::vector<std::vector<Real>> & get_psi()
  { return psi_map; }
 
  std::vector<std::vector<Real>> & get_dpsidxi()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dpsidxi_map; }
 
  const std::vector<std::vector<Real>> & get_dpsidxi() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dpsidxi_map; }
 
  std::vector<std::vector<Real>> & get_dpsideta()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dpsideta_map; }
 
  const std::vector<std::vector<Real>> & get_dpsideta() const
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dpsideta_map; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<std::vector<Real>> & get_d2psidxi2()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psidxi2_map; }
 
  const std::vector<std::vector<Real>> & get_d2psidxi2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psidxi2_map; }
 
  std::vector<std::vector<Real>> & get_d2psidxideta()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psidxideta_map; }
 
  const std::vector<std::vector<Real>> & get_d2psidxideta() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psidxideta_map; }
 
  std::vector<std::vector<Real>> & get_d2psideta2()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psideta2_map; }
 
 
  const std::vector<std::vector<Real>> & get_d2psideta2() const
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2psideta2_map; }
 
#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<std::vector<Real>> & get_phi_map()
  { libmesh_assert(!calculations_started || calculate_xyz);
    calculate_xyz = true; return phi_map; }
 
  std::vector<std::vector<Real>> & get_dphidxi_map()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphidxi_map; }
 
  std::vector<std::vector<Real>> & get_dphideta_map()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphideta_map; }
 
  std::vector<std::vector<Real>> & get_dphidzeta_map()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return dphidzeta_map; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  std::vector<std::vector<Real>> & get_d2phidxi2_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phidxi2_map; }
 
  std::vector<std::vector<Real>> & get_d2phidxideta_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phidxideta_map; }
 
  std::vector<std::vector<Real>> & get_d2phidxidzeta_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phidxidzeta_map; }
 
  std::vector<std::vector<Real>> & get_d2phideta2_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phideta2_map; }
 
  std::vector<std::vector<Real>> & get_d2phidetadzeta_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phidetadzeta_map; }
 
  std::vector<std::vector<Real>> & get_d2phidzeta2_map()
  { libmesh_assert(!calculations_started || calculate_d2xyz);
    calculate_d2xyz = true; return d2phidzeta2_map; }
#endif
 
  /* FIXME: PB: This function breaks encapsulation! Needed in FE<>::reinit and
     InfFE<>::reinit. Not sure yet if the algorithm can be redone to avoid this. */
  std::vector<Real> & get_JxW()
  { libmesh_assert(!calculations_started || calculate_dxyz);
    calculate_dxyz = true; return JxW; }
 
  void add_calculations();
 
  void set_jacobian_tolerance(Real tol) { jacobian_tolerance = tol; }
 
protected:
 
  void determine_calculations()
  {
    calculations_started = true;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
    // Second derivative calculations currently have first derivative
    // calculations as a prerequisite
    if (calculate_d2xyz)
      calculate_dxyz = true;
#endif
  }
 
  void resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp);
 
  Real dxdxi_map(const unsigned int p) const   { return dxyzdxi_map[p](0); }
 
  Real dydxi_map(const unsigned int p) const   { return dxyzdxi_map[p](1); }
 
  Real dzdxi_map(const unsigned int p) const   { return dxyzdxi_map[p](2); }
 
  Real dxdeta_map(const unsigned int p) const  { return dxyzdeta_map[p](0); }
 
  Real dydeta_map(const unsigned int p) const  { return dxyzdeta_map[p](1); }
 
  Real dzdeta_map(const unsigned int p) const  { return dxyzdeta_map[p](2); }
 
  Real dxdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](0); }
 
  Real dydzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](1); }
 
  Real dzdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](2); }
 
  std::vector<Point> xyz;
 
  std::vector<RealGradient> dxyzdxi_map;
 
  std::vector<RealGradient> dxyzdeta_map;
 
  std::vector<RealGradient> dxyzdzeta_map;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<RealGradient> d2xyzdxi2_map;
 
  std::vector<RealGradient> d2xyzdxideta_map;
 
  std::vector<RealGradient> d2xyzdeta2_map;
 
  std::vector<RealGradient> d2xyzdxidzeta_map;
 
  std::vector<RealGradient> d2xyzdetadzeta_map;
 
  std::vector<RealGradient> d2xyzdzeta2_map;
 
#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<Real> dxidx_map;
 
  std::vector<Real> dxidy_map;
 
  std::vector<Real> dxidz_map;
 
 
  std::vector<Real> detadx_map;
 
  std::vector<Real> detady_map;
 
  std::vector<Real> detadz_map;
 
 
  std::vector<Real> dzetadx_map;
 
  std::vector<Real> dzetady_map;
 
  std::vector<Real> dzetadz_map;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  std::vector<std::vector<Real>> d2xidxyz2_map;
 
  std::vector<std::vector<Real>> d2etadxyz2_map;
 
  std::vector<std::vector<Real>> d2zetadxyz2_map;
#endif
 
  std::vector<std::vector<Real>> phi_map;
 
  std::vector<std::vector<Real>> dphidxi_map;
 
  std::vector<std::vector<Real>> dphideta_map;
 
  std::vector<std::vector<Real>> dphidzeta_map;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<std::vector<Real>> d2phidxi2_map;
 
  std::vector<std::vector<Real>> d2phidxideta_map;
 
  std::vector<std::vector<Real>> d2phidxidzeta_map;
 
  std::vector<std::vector<Real>> d2phideta2_map;
 
  std::vector<std::vector<Real>> d2phidetadzeta_map;
 
  std::vector<std::vector<Real>> d2phidzeta2_map;
 
#endif
 
  std::vector<std::vector<Real>> psi_map;
 
  std::vector<std::vector<Real>> dpsidxi_map;
 
  std::vector<std::vector<Real>> dpsideta_map;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<std::vector<Real>> d2psidxi2_map;
 
  std::vector<std::vector<Real>> d2psidxideta_map;
 
  std::vector<std::vector<Real>> d2psideta2_map;
 
#endif
 
  std::vector<std::vector<Point>> tangents;
 
  std::vector<Point> normals;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  std::vector<Real> curvatures;
 
#endif
 
  std::vector<Real> jac;
 
  std::vector<Real> JxW;
 
  mutable bool calculations_started;
 
  mutable bool calculate_xyz;
 
  mutable bool calculate_dxyz;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  mutable bool calculate_d2xyz;
 
#endif
 
  Real jacobian_tolerance;
 
private:
  void compute_inverse_map_second_derivs(unsigned p);
 
  std::vector<const Node *> _elem_nodes;
 
  static Threads::spin_mutex _point_inv_err_mutex;
};
 
 
 
inline void
FEMap::init_reference_to_physical_map(unsigned int dim,
                                      const std::vector<Point> & qp,
                                      const Elem * elem)
{
  switch (dim)
  {
  case 0:
    this->init_reference_to_physical_map<0>(qp, elem);
    break;
  case 1:
    this->init_reference_to_physical_map<1>(qp, elem);
    break;
  case 2:
    this->init_reference_to_physical_map<2>(qp, elem);
    break;
  case 3:
    this->init_reference_to_physical_map<3>(qp, elem);
    break;
  default:
    libmesh_error();
  }
}
 
}
 
#endif // LIBMESH_FE_MAP_H