fe_map.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2024 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;
};



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