fe_base.h

一个用来实现偏微分方程中网格的计算库

// $Id: fe_base.h 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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 __fe_base_h__
#define __fe_base_h__

// C++ includes
#include <vector>

// Local includes
#include "reference_counted_object.h"
#include "point.h"
#include "vector_value.h"
#include "enum_elem_type.h"
#include "fe_type.h"
#include "auto_ptr.h"

#ifdef ENABLE_SECOND_DERIVATIVES
#include "tensor_value.h"
#endif


// forward declarations
template <typename T> class DenseMatrix;
template <typename T> class DenseVector;
class BoundaryInfo;
class DofConstraints;
class DofMap;
class Elem;
class MeshBase;
class PeriodicBoundaries;
template <typename T> class NumericVector;
class QBase;

#ifdef ENABLE_INFINITE_ELEMENTS

template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
class InfFE;

#endif



/**
 * This class forms the foundation from which generic finite
 * elements may be derived.  In the current implementation the
 * templated derived class \p FE offers a wide variety of commonly 
 * used finite element concepts.  Check there for details.
 * Use the \p FEBase::build() method to create an object of any of 
 * the derived classes.
 * Note that the amount of virtual members is kept to a minimum,
 * and the sophisticated template scheme of \p FE is quite
 * likely to offer acceptably fast code.
 *
 * All calls to static members of the \p FE classes should be
 * requested through the \p FEInterface.  This interface class
 * offers sort-of runtime polymorphism for the templated finite
 * element classes.  Even internal library classes, like \p DofMap,
 * request the number of dof's through this interface class.
 * Note that this also enables the co-existence of various
 * element-based schemes.  
 * This class is well 'at the heart' of the library, so 
 * things in here should better remain unchanged. 
 *
 * @author Benjamin S. Kirk, 2002
 */

// ------------------------------------------------------------
// FEBase class definition
class FEBase : public ReferenceCountedObject<FEBase>
{
protected:

  /**
   * Constructor.  Optionally initializes required data
   * structures.  Protected so that this base class
   * cannot be explicitly instantiated.
   */
  FEBase (const unsigned int dim,
	  const FEType& fet);
  
public:
  
  /**
   * Destructor.
   */
  virtual ~FEBase();

  /**
   * Builds a specific finite element type.  A \p AutoPtr<FEBase> is
   * returned to prevent a memory leak. This way the user need not
   * remember to delete the object.
   */
  static AutoPtr<FEBase> build (const unsigned int dim,
				const FEType& type); 
  
#ifdef ENABLE_INFINITE_ELEMENTS

  /**
   * Builds a specific infinite element type.  A \p AutoPtr<FEBase> is
   * returned to prevent a memory leak. This way the user need not
   * remember to delete the object.
   */
  static AutoPtr<FEBase> build_InfFE (const unsigned int dim,
				      const FEType& type); 

#endif

  /**
   * This is at the core of this class. Use this for each
   * new element in the mesh.  Reinitializes the requested physical 
   * element-dependent data based on the current element 
   * \p elem. By default the element data computed at the quadrature
   * points specified by the quadrature rule \p qrule, but any set
   * of points on the reference element may be specified in the optional
   * argument \p pts.
   *
   * Note that the FE classes decide which data to initialize based on
   * which accessor functions such as \p get_phi() or \p get_d2phi() have
   * been called, so all such accessors should be called before the first
   * \p reinit().
   */
  virtual void reinit (const Elem* elem,
		       const std::vector<Point>* const pts = NULL) = 0;
    
  /**
   * Reinitializes all the physical element-dependent data based on
   * the \p side of the element \p elem.  The \p tolerance paremeter
   * is passed to the involved call to \p inverse_map().
   */
  virtual void reinit (const Elem* elem,
		       const unsigned int side,
		       const Real tolerance = TOLERANCE) = 0;
  
  /**
   * Reinitializes all the physical element-dependent data based on
   * the \p edge of the element \p elem.  The \p tolerance paremeter
   * is passed to the involved call to \p inverse_map().
   */
  virtual void edge_reinit (const Elem* elem,
		            const unsigned int edge,
			    const Real tolerance = TOLERANCE) = 0;
  
  /**
   * @returns true if the point p is located on the reference element
   * for element type t, false otherwise.  Since we are doing floating
   * point comparisons here the parameter \p eps can be specified to
   * indicate a tolerance.  For example, \f$ x \le 1 \f$  becomes
   * \f$ x \le 1 + \epsilon \f$. 
   */
  static bool on_reference_element(const Point& p,
				   const ElemType t,
				   const Real eps = TOLERANCE);

#ifdef ENABLE_AMR

  /**
   * Computes the constraint matrix contributions (for
   * non-conforming adapted meshes) corresponding to
   * variable number \p var_number, using generic
   * projections.
   */
  static void compute_proj_constraints (DofConstraints &constraints,
                                        DofMap &dof_map,
                                        const unsigned int variable_number,
                                        const Elem* elem);

  /**
   * Creates a local projection on \p coarse_elem, based on the
   * DoF values in \p global_vector for it's children.
   */

  static void coarsened_dof_values(const NumericVector<Number> &global_vector,
			           const DofMap &dof_map,
                                   const Elem *coarse_elem,
			           DenseVector<Number> &coarse_dofs,
			           const unsigned int var,
			           const bool use_old_dof_indices = false);

#endif // #ifdef ENABLE_AMR

#ifdef ENABLE_PERIODIC

  /**
   * Computes the constraint matrix contributions (for
   * meshes with periodic boundary conditions) corresponding to
   * variable number \p var_number, using generic projections.
   */
  static void compute_periodic_constraints (DofConstraints &constraints,
                                            DofMap &dof_map,
                                            PeriodicBoundaries &boundaries,
                                            const MeshBase& mesh,
                                            const unsigned int variable_number,
                                            const Elem* elem);

#endif // ENABLE_PERIODIC

  /**
   * @returns the \p xyz spatial locations of the quadrature
   * points on the element.
   */    
  const std::vector<Point>& get_xyz() const
  { return xyz; }
  
  /**
   * @returns the shape function values at the quadrature points
   * on the element.
   */    
  const std::vector<std::vector<Real> >& get_phi() const
  { libmesh_assert(!calculations_started || calculate_phi);
    calculate_phi = true; return phi; }
  
  /**
   * @returns the element Jacobian times the quadrature weight for
   * each quadrature point.
   */    
  const std::vector<Real>& get_JxW() const
  { return JxW; }

  /**
   * @returns the shape function derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<RealGradient> >& get_dphi() const
  { libmesh_assert(!calculations_started || calculate_dphi); 
    calculate_dphi = true; return dphi; }
  
  /**
   * @returns the shape function x-derivative at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_dphidx() const
  { libmesh_assert(!calculations_started || calculate_dphi); 
    calculate_dphi = true; return dphidx; }
  
  /**
   * @returns the shape function y-derivative at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_dphidy() const
  { libmesh_assert(!calculations_started || calculate_dphi); 
    calculate_dphi = true; return dphidy; }
  
  /**
   * @returns the shape function z-derivative at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_dphidz() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = true; return dphidz; }

#ifdef ENABLE_SECOND_DERIVATIVES

  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<RealTensor> >& get_d2phi() const
  { libmesh_assert(!calculations_started || calculate_d2phi); 
    calculate_d2phi = true; return d2phi; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidx2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidx2; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidxdy() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidxdy; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidxdz() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidxdz; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidy2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidy2; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidydz() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidydz; }
  
  /**
   * @returns the shape function second derivatives at the quadrature
   * points.
   */
  const std::vector<std::vector<Real> >& get_d2phidz2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = true; return d2phidz2; }

#endif
  
  /**
   * @returns the element tangents in xi-direction at the quadrature
   * points.
   */
  const std::vector<RealGradient>& get_dxyzdxi() const
  { return dxyzdxi_map; }

  /**
   * @returns the element tangents in eta-direction at the quadrature
   * points.
   */
  const std::vector<RealGradient>& get_dxyzdeta() const
  { return dxyzdeta_map; }

  /**
   * @returns the element tangents in zeta-direction at the quadrature
   * points.
   */
  const std::vector<RealGradient>& get_dxyzdzeta() const
  { return dxyzdzeta_map; }
  
  /**
   * @returns the second partial derivatives in xi.
   */
  const std::vector<RealGradient>& get_d2xyzdxi2() const
  { return d2xyzdxi2_map; }

  /**
   * @returns the second partial derivatives in eta.
   */
  const std::vector<RealGradient>& get_d2xyzdeta2() const
  { return d2xyzdeta2_map; }

#ifdef ENABLE_SECOND_DERIVATIVES
  
  /**
   * @returns the second partial derivatives in zeta.
   */
  const std::vector<RealGradient>& get_d2xyzdzeta2() const
  { return d2xyzdzeta2_map; }
  
#endif
  
  /**
   * @returns the second partial derivatives in xi-eta.
   */
  const std::vector<RealGradient>& get_d2xyzdxideta() const
  { return d2xyzdxideta_map; }

#ifdef ENABLE_SECOND_DERIVATIVES
  
  /**
   * @returns the second partial derivatives in xi-zeta.
   */
  const std::vector<RealGradient>& get_d2xyzdxidzeta() const
  { return d2xyzdxidzeta_map; }

  /**
   * @returns the second partial derivatives in eta-zeta.
   */
  const std::vector<RealGradient>& get_d2xyzdetadzeta() const
  { return d2xyzdetadzeta_map; }

#endif
  
  /**
   * @returns the dxi/dx entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_dxidx() const
  { return dxidx_map; }

  /**
   * @returns the dxi/dy entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_dxidy() const
  { return dxidy_map; }

  /**
   * @returns the dxi/dz entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_dxidz() const
  { return dxidz_map; }

  /**
   * @returns the deta/dx entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_detadx() const