inf_fe.h

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

// $Id: inf_fe.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 __inf_fe_h__
#define __inf_fe_h__

// C++ includes

// Local includes
#include "libmesh_config.h"
#ifdef ENABLE_INFINITE_ELEMENTS

#include "fe_base.h"


// forward declarations
class Elem;
class FEComputeData;



/**
 * A specific instatiation of the \p FEBase class. This
 * class is templated, and specific template instantiations
 * will result in different Infinite Element families, similar
 * to the \p FE class.  \p InfFE builds a \p FE<Dim-1,T_base>,
 * and most of the requests related to the base are handed over
 * to this object.  All methods related to the radial part 
 * are collected in the nested class \p Radial.  Similarly,
 * most of the static methods concerning base approximation
 * are contained in \p Base.
 * 
 * Having different shape approximation families in radial direction 
 * introduces the requirement for an additional \p Order in this
 * class. Therefore, the \p FEType internals change when infinite
 * elements are enabled. 
 * When the specific infinite element type is not known at compile
 * time, use the \p FEBase::build() member to create abstract 
 * (but still optimized) infinite elements at run time.
 *
 * The @e node numbering scheme is the one from the current
 * infinite element.  Each node in the base holds exactly
 * the same number of dofs as an adjacent conventional \p FE
 * would contain.  The nodes further out hold the additional
 * dof necessary for radial approximation.  The order of the outer nodes'
 * components is such that the radial shapes have highest 
 * priority, followed by the base shapes.
 *
 * \author Daniel Dreyer
 * \date 2003
 * \version $Revision: 2789 $
 */

//-------------------------------------------------------------
// InfFE class definition
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
class InfFE : public FEBase
{

/*
 * Protect the nested class
 */
protected:

  /**
   * Infinite elements are in some sense directional, compared
   * to conventional finite elements.  All methods related 
   * to the radial part, which extends perpendicular from the base,
   * are collected in this nested class.  This class offers
   * static methods, which are only available to \p InfFE members.
   *
   * \author Daniel Dreyer
   * \date 2003
   * \version $Revision: 2789 $
   */
  //-------------------------------------------------------------
  // InfFE::Radial class definition
  class Radial
  {
  private:

    /**
     * Never use an object of this type.
     */
    Radial() {}

  public:

    /**
     * @returns the decay in radial direction of
     * the \p Dim dimensional infinite element.
     */
    static Real decay (const Real v);

    /**
     * @returns the first (local) derivative of the
     * decay in radial direction of the infinite element.
     */
    static Real decay_deriv (const Real) { return -.5; }

    /**
     * @returns the radial weight D, used as an additional weight
     * for the test function, evaluated at local radial coordinate \p v.
     */
    static Real D (const Real v) { return (1.-v)*(1.-v)/4.; }

    /**
     * @returns the first (local) radial derivative of the radial weight D.
     */
    static Real D_deriv (const Real v) { return (v-1.)/2.; }

   /**
     * @returns the Order of the mapping functions
     * in radial direction. Currently, this is @e always \p FIRST.
     */
    static Order mapping_order() { return FIRST; }

    /**
     * @returns the number of shape functions in radial direction
     * associated with this infinite element.
     * Either way, if the modes are stored as nodal dofs (\p n_dofs_at_node) 
     * or as element dofs (\p n_dofs_per_elem), in each case we have the
     * same number of modes in radial direction. Note that for the case of 1D
     * infinite elements, in the base the dof-per-node scheme is used.
     * 
     * From the formulation of the infinite elements, we have
     * 1 mode, when \p o_radial=CONST.
     * Therefore, we have a total of \p o_radial+1 modes in radial direction.
     */
    static unsigned int n_dofs (const Order o_radial)
	{ return static_cast<unsigned int>(o_radial)+1; }

    /**
     * @returns the number of dofs in radial direction on "onion slice" 
     * \p n (either 0 or 1) for an infinite element of type \p inf_elem_type and 
     * radial order \p o_radial.
     *
     * Currently, the first radial mode is associated with the nodes in
     * the base.  All higher radial modes are associated with
     * the physically existing nodes further out.
     */
    static unsigned int n_dofs_at_node (const Order o_radial,
					const unsigned int n_onion);

    /**
     * @returns the number of modes in radial direction interior to the element,
     * not associated with any interior nodes.
     * Note that these modes are a discontinuous approximation, therefore
     * we have no special formulation for coupling in the base, like in the 
     * case of associating (possibly) multiple dofs per (outer) node.
     */
    static unsigned int n_dofs_per_elem (const Order o_radial)
	{ return static_cast<unsigned int>(o_radial)+1; }

  };



  /**
   * This nested class contains most of the static methods related 
   * to the base part of an infinite element.  Only static members
   * are provided, and these should only be accessible from within \p InfFE.
   *
   * \author Daniel Dreyer
   * \date 2003
   * \version $Revision: 2789 $
   */
  //-------------------------------------------------------------
  // InfFE::Base class definition
  class Base
  {
  private:

    /**
     * Never use an object of this type.
     */
    Base() {}

  public:

    /**
     * Build the base element of an infinite element.  Be careful,
     * this method allocates memory!  So be sure to delete the
     * new element afterwards.
     */
    static Elem* build_elem (const Elem* inf_elem);

    /**
     * @returns the base element associated to
     * \p type.  This is, for example, \p TRI3 for
     * \p INFPRISM6.
     */
    static ElemType get_elem_type (const ElemType type);

    /**
     * @returns the number of shape functions used in the
     * mapping in the @e base element of type \p base_elem_type
     * mapped with order \p base_mapping_order
     */
    static unsigned int n_base_mapping_sf (const ElemType base_elem_type,
					   const Order base_mapping_order);

  };







public:

  //-------------------------------------------------------------
  // InfFE continued

  /**
   * Constructor.
   * Initializes some data structures.  Builds a \p FE<Dim-1,T_base>
   * object to handle  approximation in the base, so that
   * there is no need to template \p InfFE<Dim,T_radial,T_map> also with 
   * respect to the base approximation \p T_base.
   * 
   * The same remarks concerning compile-time optimization for 
   * \p FE also hold for \p InfFE.  Use the 
   * \p FEBase::build_InfFE(const unsigned int, const FEType&) 
   * method to build specific instantiations of \p InfFE at
   * run time.
   */
  InfFE(const FEType& fet);

  /**
   * Desctructor.  Clean up.
   */
  ~InfFE();





  //-------------------------------------------------------------
  // The static public members for access from FEInterface etc
  /**
   * @returns the value of the \f$ i^{th} \f$ shape function at
   * point \p p.  This method lets you specify the relevant
   * data directly, and is therefore allowed to be static.
   * Note that this class member is by far not as efficient as 
   * its counterpart in \p FE<Dim,T>, and is @e not employed
   * in the \p reinit() cycle.
   * Also note that this method does @e not return physically
   * correct shapes, instead use \p compute_data().  The \p shape()
   * methods should only be used for mapping.
   */
  static Real shape(const FEType& fet,
		    const ElemType t,
		    const unsigned int i,
		    const Point& p);

  /**
   * @returns the value of the \f$ i^{th} \f$ shape function at
   * point \p p.  This method lets you specify the relevant
   * data directly, and is therefore allowed to be static.
   * Note that this class member is not as efficient as its 
   * counterpart in \p FE<Dim,T>, and is @e not employed 
   * in the \p reinit() cycle.
   * Also note that this method does @e not return physically
   * correct shapes, instead use \p compute_data().  The \p shape()
   * methods should only be used for mapping.
   */
  static Real shape(const FEType& fet,
		    const Elem* elem,
		    const unsigned int i,
		    const Point& p);

  /**
   * Generalized version of \p shape(), takes an \p Elem*.  The \p data
   * contains both input and output parameters.  For frequency domain
   * simulations, the complex-valued shape is returned.  In time domain
   * both the computed shape, @e and the phase is returned.  Note that
   * the phase (proportional to the distance of the \p Point \p data.p
   * from the envelope) is actually a measure how @e far into the @e future
   * the results are.  Pretty weird, hm!?
   */
  static void compute_data(const FEType& fe_t,
			   const Elem* inf_elem,
			   FEComputeData& data);
  
  /**
   * @returns the number of shape functions associated with
   * a finite element of type \p t and approximation order \p o.
   */
  static unsigned int n_shape_functions (const FEType& fet,
					 const ElemType t)
    { return n_dofs(fet, t); }

  /**
   * @returns the number of shape functions associated with this
   * infinite element.  Currently, we have \p o_radial+1 modes in 
   * radial direction, and \p FE<Dim-1,T>::n_dofs(...) in the base.
   */
  static unsigned int n_dofs(const FEType& fet,
			     const ElemType inf_elem_type);

  /**
   * @returns the number of dofs at infinite element @e node \p n 
   * (not dof!) for an element of type \p t and order \p o.  
   */
  static unsigned int n_dofs_at_node(const FEType& fet,
				     const ElemType inf_elem_type,
				     const unsigned int n);

  /**
   * @returns the number of dofs interior to the element,
   * not associated with any interior nodes.
   */
  static unsigned int n_dofs_per_elem(const FEType& fet,
				      const ElemType inf_elem_type);

  /**
   * @returns the continuity of the element.
   */
  virtual FEContinuity get_continuity() const
    { return C_ZERO; }  // FIXME - is this true??

  /**
   * @returns true if the element's higher order shape functions are
   * hierarchic
   */
  virtual bool is_hierarchic() const
    { return false; }  // FIXME - Inf FEs don't handle p elevation yet

  /**
   * Usually, this method would build the nodal soln from the 
   * element soln.  But infinite elements require additional 
   * simulation-specific data to compute physically correct 
   * results.  Use \p compute_data() to compute results.  For
   * compatibility an empty vector is returned.
   */
  static void nodal_soln(const FEType& fet,
			 const Elem* elem, 
			 const std::vector<Number>& elem_soln,
			 std::vector<Number>& nodal_soln);

  /**
   * @returns the location (on the reference element) of the
   * point \p p located in physical space.  First, the location
   * in the base face is computed. This requires inverting the 
   * (possibly nonlinear) transformation map in the base, so it is 
   * not trivial. The optional parameter \p tolerance defines
   * how close is "good enough."  The map inversion iteration
   * computes the sequence \f$ \{ p_n \} \f$, and the iteration is
   * terminated when \f$ \|p - p_n\| < \mbox{\texttt{tolerance}} \f$.
   * Once the base face point is determined, the radial local
   * coordinate is directly evaluated.
   * If \p interpolated is true, the interpolated distance from the
   * base element to the infinite element origin is used for the map
   * in radial direction.
   */
  static Point inverse_map (const Elem* elem,
			    const Point& p,
			    const Real tolerance = TOLERANCE,
			    const bool secure = true,
			    const bool interpolated = true);


  /**
   * Takes a number points in physical space (in the \p physical_points
   * vector) and finds their location on the reference element for the
   * input element \p elem.  The values on the reference element are
   * returned in the vector \p reference_points
   */
  static void inverse_map (const Elem* elem,
			   const std::vector<Point>& physical_points,
			   std::vector<Point>&       reference_points,
			   const Real tolerance = TOLERANCE,
			   const bool secure = true);


  //-------------------------------------------------------------
  // The work-horses of InfFE. These are often used during matrix assembly
  /**
   * This is at the core of this class. Use this for each
   * new element in the mesh.  Reinitializes all the physical 
   * element-dependent data based on the current element 
   * \p elem.
   */
  virtual void reinit (const Elem* elem,
		       const std::vector<Point>* const pts);
    
  /**
   * Not implemented yet.  Reinitializes all the physical 
   * element-dependent data based on the \p side of an infinite 
   * element.
   */
  virtual void reinit (const Elem* elem,
		       const unsigned int side,
		       const Real tolerance = TOLERANCE);

  /**
   * Not implemented yet.  Reinitializes all the physical 
   * element-dependent data based on the \p edge of an infinite 
   * element.
   */
  virtual void edge_reinit (const Elem* elem,
		            const unsigned int edge,
			    const Real tolerance = TOLERANCE);

  /**
   * The use of quadrature rules with the \p InfFE class is somewhat
   * different from the approach of the \p FE class.  While the
   * \p FE class requires an appropriately initialized quadrature
   * rule object, and simply uses it, the \p InfFE class requires only