fe_map.h

Go to the documentation of this file.

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

#include "libmesh/reference_counted_object.h"
#include "libmesh/point.h"
#include "libmesh/vector_value.h"
#include "libmesh/enum_elem_type.h"
#include "libmesh/fe_type.h"
#include "libmesh/auto_ptr.h"
#include "libmesh/node.h"

namespace libMesh
{

// forward declarations
class Elem;

class FEMap
{
public:

  FEMap();
  virtual ~FEMap(){}

  static UniquePtr<FEMap> build(FEType fe_type);

  template<unsigned int Dim>
  void init_reference_to_physical_map(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<Node*>& elem_nodes);

  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);

  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);

  const std::vector<Point>& get_xyz() const
  { return xyz; }

  const std::vector<Real>& get_jacobian() const
  { return jac; }

  const std::vector<Real>& get_JxW() const
  { return JxW; }

  const std::vector<RealGradient>& get_dxyzdxi() const
  { return dxyzdxi_map; }

  const std::vector<RealGradient>& get_dxyzdeta() const
  { return dxyzdeta_map; }

  const std::vector<RealGradient>& get_dxyzdzeta() const
  { return dxyzdzeta_map; }

  const std::vector<RealGradient>& get_d2xyzdxi2() const
  { return d2xyzdxi2_map; }

  const std::vector<RealGradient>& get_d2xyzdeta2() const
  { return d2xyzdeta2_map; }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  const std::vector<RealGradient>& get_d2xyzdzeta2() const
  { return d2xyzdzeta2_map; }

#endif

  const std::vector<RealGradient>& get_d2xyzdxideta() const
  { return d2xyzdxideta_map; }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  const std::vector<RealGradient>& get_d2xyzdxidzeta() const
  { return d2xyzdxidzeta_map; }

  const std::vector<RealGradient>& get_d2xyzdetadzeta() const
  { return d2xyzdetadzeta_map; }

#endif

  const std::vector<Real>& get_dxidx() const
  { return dxidx_map; }

  const std::vector<Real>& get_dxidy() const
  { return dxidy_map; }

  const std::vector<Real>& get_dxidz() const
  { return dxidz_map; }

  const std::vector<Real>& get_detadx() const
  { return detadx_map; }

  const std::vector<Real>& get_detady() const
  { return detady_map; }

  const std::vector<Real>& get_detadz() const
  { return detadz_map; }

  const std::vector<Real>& get_dzetadx() const
  { return dzetadx_map; }

  const std::vector<Real>& get_dzetady() const
  { return dzetady_map; }

  const std::vector<Real>& get_dzetadz() const
  { return dzetadz_map; }

  const std::vector<std::vector<Real> >& get_psi() const
  { return psi_map; }

  const std::vector<std::vector<Real> >& get_phi_map() const
  { return phi_map; }

  const std::vector<std::vector<Real> >& get_dphidxi_map() const
  { return dphidxi_map; }

  const std::vector<std::vector<Real> >& get_dphideta_map() const
  { return dphideta_map; }

  const std::vector<std::vector<Real> >& get_dphidzeta_map() const
  { return dphidzeta_map; }

  const std::vector<std::vector<Point> >& get_tangents() const
  { return tangents; }

  const std::vector<Point>& get_normals() const
  { return normals; }

  const std::vector<Real>& get_curvatures() const
  { return curvatures;}

  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()
  { return dpsidxi_map; }

  std::vector<std::vector<Real> >& get_dpsideta()
  { return dpsideta_map; }

  std::vector<std::vector<Real> >& get_d2psidxi2()
  { return d2psidxi2_map; }

  std::vector<std::vector<Real> >& get_d2psidxideta()
  { return d2psidxideta_map; }

  std::vector<std::vector<Real> >& get_d2psideta2()
  { return d2psideta2_map; }

  std::vector<std::vector<Real> >& get_phi_map()
  { return phi_map; }

  std::vector<std::vector<Real> >& get_dphidxi_map()
  { return dphidxi_map; }

  std::vector<std::vector<Real> >& get_dphideta_map()
  { return dphideta_map; }

  std::vector<std::vector<Real> >& get_dphidzeta_map()
  { return dphidzeta_map; }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  std::vector<std::vector<Real> >& get_d2phidxi2_map()
  { return d2phidxi2_map; }

  std::vector<std::vector<Real> >& get_d2phidxideta_map()
  { return d2phidxideta_map; }

  std::vector<std::vector<Real> >& get_d2phidxidzeta_map()
  { return d2phidxidzeta_map; }

  std::vector<std::vector<Real> >& get_d2phideta2_map()
  { return d2phideta2_map; }

  std::vector<std::vector<Real> >& get_d2phidetadzeta_map()
  { return d2phidetadzeta_map; }

  std::vector<std::vector<Real> >& get_d2phidzeta2_map()
  { 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()
  { return JxW; }

protected:

  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;

  std::vector<RealGradient> d2xyzdxi2_map;

  std::vector<RealGradient> d2xyzdxideta_map;

  std::vector<RealGradient> d2xyzdeta2_map;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  std::vector<RealGradient> d2xyzdxidzeta_map;

  std::vector<RealGradient> d2xyzdetadzeta_map;

  std::vector<RealGradient> d2xyzdzeta2_map;

#endif

  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;

  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;

  std::vector<std::vector<Real> > d2psidxi2_map;

  std::vector<std::vector<Real> > d2psidxideta_map;

  std::vector<std::vector<Real> > d2psideta2_map;

  std::vector<std::vector<Point> > tangents;

  std::vector<Point> normals;

  std::vector<Real> curvatures;

  std::vector<Real> jac;

  std::vector<Real> JxW;
};

}

#endif // LIBMESH_FE_MAP_H