inf_fe.h

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

   * the quadrature rule object of the current \p FE class.
   * From this \p QBase*, it determines the necessary data,
   * and @e builds two appropriate quadrature classes, one for radial,
   * and another for base integration, using the convenient
   * \p QBase::build() method.
   */
  virtual void attach_quadrature_rule (QBase* q);

  /**
   * @returns the number of shape functions associated with
   * this infinite element.
   */
  virtual unsigned int n_shape_functions () const
  { return _n_total_approx_sf; }
  
  /**
   * @returns the total number of quadrature points.  Call this
   * to get an upper bound for the \p for loop in your simulation
   * for matrix assembly of the current element.
   */
  virtual unsigned int n_quadrature_points () const
  { libmesh_assert (radial_qrule != NULL); return _n_total_qp; }


protected:

  //-------------------------------------------------------------
  // static members used by the "work-horses"

  /**
   * @returns the value of the \f$ i^{th} \f$ polynomial evaluated
   * at \p v.  This method provides the approximation
   * in radial direction for the overall shape functions,
   * which is defined in \p InfFE::shape().
   * This method is allowed to be static, since it is independent
   * of dimension and base_family.  It is templated, though,
   * w.r.t. to radial \p FEFamily.
   *
   * Specialized for \p T_radial=INFINITE_MAP, this function returns 
   * the value of the \f$ i^{th} \f$ @e mapping shape function
   * in radial direction evaluated at \p v.  Currently, only one specific 
   * mapping shape is used.  Namely the one by Marques JMMC, Owen DRJ:
   * Infinite elements in quasi-static materially nonlinear problems,
   * @e Computers @e and @e Structures, 1984.
   */
  static Real eval(const Real v,
		   const Order o_radial,
		   const unsigned int i);
  
  /**
   * @returns the value of the first derivative of the
   * \f$ i^{th} \f$ polynomial at coordinate \p v.
   * See \p eval for details.
   */
  static Real eval_deriv(const Real v,
			 const Order o_radial,
			 const unsigned int i);



  //-------------------------------------------------------------
  // Non-static members used by the "work-horses"
  /**
   * Updates the protected member \p base_elem to the appropriate base element
   * for the given \p inf_elem.
   */
  void update_base_elem (const Elem* inf_elem);

  /**
   * Do not use this derived member in \p InfFE<Dim,T_radial,T_map>.
   */
  virtual void init_base_shape_functions(const std::vector<Point>&,
					 const Elem*)
  { libmesh_error(); }

  /**
   * Some of the member data only depend on the radial part of the
   * infinite element.  The parts that only change when the radial
   * order changes, are initialized here.  
   */
  void init_radial_shape_functions(const Elem* inf_elem);

  /** 
   * Initialize all the data fields like \p weight, \p mode, 
   * \p phi, \p dphidxi, \p dphideta, \p dphidzeta, etc.
   * for the current element.  This method prepares the data
   * related to the base part, and some of the combined fields.
   */
  void init_shape_functions(const Elem* inf_elem);

  /** 
   * Not implemented yet.  Initialize all the data fields like \p weight, 
   * \p phi, etc for the side \p s.
   */  
  void init_face_shape_functions (const std::vector<Point>& qp,
				  const Elem* side);

  /** 
   * Combines the shape functions, which were formed in
   * \p init_shape_functions(Elem*), with geometric data.
   * Has to be called every time the geometric configuration
   * changes.  Afterwards, the fields are ready to be used
   * to compute global derivatives, the jacobian etc, see
   * \p FEBase::compute_map().
   */
  void combine_base_radial(const Elem* inf_elem);

  /** 
   * After having updated the jacobian and the transformation
   * from local to global coordinates in FEBase::compute_map(),
   * the first derivatives of the shape functions are 
   * transformed to global coordinates, giving \p dphi,
   * \p dphidx/y/z, \p dphasedx/y/z, \p dweight. This method
   * should barely be re-defined in derived classes, but
   * still should be usable for children. Therefore, keep
   * it protected.
   * Overloaded method from the \p FEBase version.
   */
  virtual void compute_shape_functions(const Elem*);



  //-------------------------------------------------------------
  // Miscellaneous static members

  /**
   * @returns the location (in physical space) of the point
   * \p p located on the reference element.
   */
  static Point map (const Elem* inf_elem,
		    const Point& reference_point);

  /**
   * Computes the indices in the base \p base_node and in radial 
   * direction \p radial_node (either 0 or 1) associated to the 
   * node \p outer_node_index of an infinite element of type
   * \p inf_elem_type.
   */
  static void compute_node_indices (const ElemType inf_elem_type,
				    const unsigned int outer_node_index,
				    unsigned int& base_node,
				    unsigned int& radial_node);

  /**
   * Does the same as \p compute_node_indices(), but stores
   * the maps for the current element type.  Provided the
   * infinite element type changes seldom, this is probably
   * faster than using \p compute_node_indices () alone.
   * This is possible since the number of @e nodes is not likely
   * to change.
   */
  static void compute_node_indices_fast (const ElemType inf_elem_type,
					 const unsigned int outer_node_index,
					 unsigned int& base_node,
					 unsigned int& radial_node);

  /**
   * Computes the indices of shape functions in the base \p base_shape and 
   * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further 
   * out) associated to the shape with global index \p i of an infinite element 
   * of type \p inf_elem_type.
   */
  static void compute_shape_indices (const FEType& fet,
				     const ElemType inf_elem_type,
				     const unsigned int i,
				     unsigned int& base_shape,
				     unsigned int& radial_shape);

  //--------------------------------------------------------------
  // protected members, which are not to be accessed from outside
  /**
   * the radial distance of the base nodes from the origin
   */
  std::vector<Real>  dist;

  /**
   * the additional radial weight \f$ 1/{r^2} \f$ in local coordinates,
   * over all quadrature points. The weight does not vary in base
   * direction.  However, for uniform access to the data fields from the 
   * outside, this data field is expanded to @e all quadrature points.
   */
  std::vector<Real>  dweightdv;

  /**
   * the radial decay \f$ 1/r \f$ in local coordinates.
   * Needed when setting up the overall shape functions.
   * Note that it is this decay which assures to satisfy
   * the Sommerfeld radiation condition in advance.
   */
  std::vector<Real>  som;
  /**
   * the first local derivative of the radial decay \f$ 1/r \f$ in local 
   * coordinates.  Needed when setting up the overall shape functions.
   */
  std::vector<Real>  dsomdv;

  /**
   * the radial approximation shapes in local coordinates
   * Needed when setting up the overall shape functions.
   */
  std::vector<std::vector<Real> >   mode;

  /**
   * the first local derivative of the radial approximation shapes.
   * Needed when setting up the overall shape functions.
   */
  std::vector<std::vector<Real> >   dmodedv;
  
  /**
   * the radial mapping shapes in local coordinates
   */
  std::vector<std::vector<Real> >   radial_map;


  /**
   * the first local derivative of the radial mapping shapes
   */
  std::vector<std::vector<Real> >   dradialdv_map;

  /**
   * the first local derivative (for 3D, the first in the base) 
   * of the phase term in local coordinates.
   * Needed in the overall weak form of infinite element formulations.
   */
  std::vector<Real>  dphasedxi;

  /**
   * the second local derivative (for 3D, the second in the base) 
   * of the phase term in local coordinates.
   * Needed in the overall weak form of infinite element formulations.
   */
  std::vector<Real>  dphasedeta;

  /**
   * the third local derivative (for 3D, the derivative in radial
   * direction) of the phase term in local coordinates.
   * Needed in the overall weak form of infinite element formulations.
   */
  std::vector<Real>  dphasedzeta;




  //--------------------------------------------------------------
  // numbering scheme maps

  /**
   * The internal structure of the \p InfFE
   * -- tensor product of base element times radial
   * nodes -- has to be determined from the node numbering
   * of the current infinite element.  This vector
   * maps the @e infinte \p Elem node number to the 
   * @e radial node (either 0 or 1).
   */
  std::vector<unsigned int> _radial_node_index;

  /**
   * The internal structure of the \p InfFE
   * -- tensor product of base element times radial
   * nodes -- has to be determined from the node numbering
   * of the current element.  This vector
   * maps the @e infinte \p Elem node number to the 
   * associated node in the @e base element.
   */
  std::vector<unsigned int> _base_node_index;

  /**
   * The internal structure of the \p InfFE
   * -- tensor product of base element shapes times radial
   * shapes -- has to be determined from the dof numbering
   * scheme of the current infinite element.  This vector
   * maps the infinite \p Elem dof index to the @e radial 
   * \p InfFE shape index (\p 0..radial_order+1 ).
   */
  std::vector<unsigned int> _radial_shape_index;

  /**
   * The internal structure of the \p InfFE
   * -- tensor product of base element shapes times radial
   * shapes -- has to be determined from the dof numbering
   * scheme of the current infinite element.  This vector
   * maps the infinite \p Elem dof index to the associated
   * dof in the @e base \p FE.
   */
  std::vector<unsigned int> _base_shape_index;




  //--------------------------------------------------------------
  // some more protected members

  /**
   * The number of total approximation shape functions for 
   * the current configuration
   */
  unsigned int _n_total_approx_sf;

  /**
   * The total number of quadrature points
   * for the current configuration
   */
  unsigned int _n_total_qp;

  /**
   * this vector contains the combined integration weights, so
   * that \p FEBase::compute_map() can still be used
   */
  std::vector<Real>  _total_qrule_weights;

  /**
   * The quadrature rule for the base element associated 
   * with the current infinite element
   */
  QBase* base_qrule;

  /**
   * The quadrature rule for the base element associated 
   * with the current infinite element
   */
  QBase* radial_qrule;

  /**
   * The base element associated with the
   * current infinite element
   */
  Elem* base_elem;

  /**
   * Have a \p FE<Dim-1,T_base> handy for base approximation.
   * Since this one is created using the \p FEBase::build() method,
   * the \p InfFE class is not required to be templated w.r.t.
   * to the base approximation shape.
   */
  FEBase* base_fe;

  /**
   * This \p FEType stores the characteristics for which
   * the data structures \p phi, \p phi_map etc are currently
   * initialized.  This avoids re-initializing the radial
   * part.  But note that currently @e only \p order may change,
   * neither the FE families nor \p base_order!
   */
  FEType current_fe_type;


private:

  /**
   * @returns \p false, currently not required.
   */
  virtual bool shapes_need_reinit() const;

  /**
   * When \p compute_node_indices_fast() is used, this static
   * variable remembers the element type for which the
   * static variables in  \p compute_node_indices_fast()
   * are currently set.  Using a class member for the
   * element type helps initializing it to a default value.
   */
  static ElemType _compute_node_indices_fast_current_elem_type;


#ifdef DEBUG

  /**
   * static members that are used to issue warning messages only once.
   */
  static bool _warned_for_nodal_soln;
  static bool _warned_for_shape;

#endif

  /**
   * Make all \p InfFE<Dim,T_radial,T_map> classes
   * friends of each other, so that the protected
   * \p eval() may be accessed.
   */
  template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
  friend class InfFE;

};




// ------------------------------------------------------------
// InfFE class inline members




// ------------------------------------------------------------
// InfFE::Radial class inline members

template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
inline
Real InfFE<Dim,T_radial,T_map>::Radial::decay(const Real v)
{
  switch (Dim)
  //TODO:[DD] What decay do i have in 2D and 1D?
  {
    case 3:
      return (1.-v)/2.;

    case 2:
      return 0.;

    case 1:
      return 0.;

    default:
      libmesh_error();
      return 0.;
  }
}



// ------------------------------------------------------------
// InfFE::Base class inline members



#endif //ifdef ENABLE_INFINITE_ELEMENTS


#endif