fe_base.h

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

  { return detadx_map; }

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

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

  /**
   * @returns the dzeta/dx entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_dzetadx() const
  { return dzetadx_map; }

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

  /**
   * @returns the dzeta/dz entry in the transformation
   * matrix from physical to local coordinates. 
   */
  const std::vector<Real>& get_dzetadz() const
  { return dzetadz_map; }
  
#ifdef ENABLE_INFINITE_ELEMENTS

  /**
   * @returns the global first derivative of the phase term 
   * which is used in infinite elements, evaluated at the 
   * quadrature points.  
   *
   * In case of the general finite element class \p FE this 
   * field is initialized to all zero, so that the variational 
   * formulation for an @e infinite element returns correct element 
   * matrices for a mesh using both finite and infinite elements.
   */
  const std::vector<RealGradient>& get_dphase() const
      { return dphase; }


  /**
   * @returns the multiplicative weight at each quadrature point.
   * This weight is used for certain infinite element weak 
   * formulations, so that @e weighted Sobolev spaces are
   * used for the trial function space.  This renders the
   * variational form easily computable. 
   *
   * In case of the general finite element class \p FE this 
   * field is initialized to all ones, so that the variational 
   * formulation for an @e infinite element returns correct element 
   * matrices for a mesh using both finite and infinite elements.
   */
  const std::vector<Real>& get_Sobolev_weight() const
      { return weight; }

  /**
   * @returns the first global derivative of the multiplicative 
   * weight at each quadrature point. See \p get_Sobolev_weight()
   * for details.  In case of \p FE initialized to all zero.
   */
  const std::vector<RealGradient>& get_Sobolev_dweight() const
      { return dweight; }

#endif

  
  /**
   * @returns the tangent vectors for face integration.
   */
  const std::vector<std::vector<Point> >& get_tangents() const
  { return tangents; }
  
  /**
   * @returns the normal vectors for face integration.
   */
  const std::vector<Point>& get_normals() const
  { return normals; }

  /**
   * @returns the curvatures for use in face integration.
   */
  const std::vector<Real>& get_curvatures() const
  { return curvatures;}
  
  /**
   * Provides the class with the quadrature rule.  Implement
   * this in derived classes.
   */
  virtual void attach_quadrature_rule (QBase* q) = 0;

  /**
   * @returns the total number of approximation shape functions
   * for the current element.  Useful during matrix assembly.  
   * Implement this in derived classes.
   */
  virtual unsigned int n_shape_functions () const = 0;

  /**
   * @returns the total number of quadrature points.  Useful
   * during matrix assembly.  Implement this in derived classes.
   */
  virtual unsigned int n_quadrature_points () const = 0;
  
  /**
   * @returns the element type that the current shape functions
   * have been calculated for.  Useful in determining when shape
   * functions must be recomputed.
   */
  ElemType get_type()  const { return elem_type; }

  /**
   * @returns the p refinement level that the current shape
   * functions have been calculated for.
   */
  unsigned int get_p_level() const { return _p_level; }

  /**
   * @returns the FE Type (approximation order and family) of the finite element.
   */
  FEType get_fe_type()  const { return fe_type; }

  /**
   * @returns the approximation order of the finite element.
   */
  Order get_order()  const { return static_cast<Order>(fe_type.order + _p_level); }

  /**
   * @returns the continuity level of the finite element.
   */
  virtual FEContinuity get_continuity() const = 0;

  /**
   * @returns true if the finite element's higher order shape functions are
   * hierarchic
   */
  virtual bool is_hierarchic() const = 0;

  /**
   * @returns the finite element family of this element.
   */
  FEFamily get_family()  const { return fe_type.family; }

  /**
   * Prints the Jacobian times the weight for each quadrature point.
   */ 
  void print_JxW(std::ostream& os) const; 
  
  /**
   * Prints the value of each shape function at each quadrature point.
   */ 
  void print_phi(std::ostream& os) const; 
  
  /**
   * Prints the value of each shape function's derivative
   * at each quadrature point.
   */ 
  void print_dphi(std::ostream& os) const;

#ifdef ENABLE_SECOND_DERIVATIVES

  /**
   * Prints the value of each shape function's second derivatives
   * at each quadrature point.
   */ 
  void print_d2phi(std::ostream& os) const;

#endif
  
  /**
   * Prints the spatial location of each quadrature point
   * (on the physical element).
   */ 
  void print_xyz(std::ostream& os) const;

  /**
   * Prints all the relevant information about the current element.
   */
  void print_info(std::ostream& os) const;

  /**
   * Same as above, but allows you to print to a stream.
   */
  friend std::ostream& operator << (std::ostream& os, const FEBase& fe);
  
  
protected:


    
#ifdef ENABLE_INFINITE_ELEMENTS

  /**
   * Initialize the data fields for the base of an
   * an infinite element.  Implement this in the derived 
   * class \p FE<Dim,T>.
   */
  virtual void init_base_shape_functions(const std::vector<Point>& qp,
					 const Elem* e) = 0;

#endif

  /**
   * Compute the jacobian and some other additional
   * data fields. Takes the integration weights
   * as input, along with a pointer to the element.
   */
  virtual void compute_map(const std::vector<Real>& qw,
		   const Elem* e);
  
  /**
   * Compute the jacobian and some other additional
   * data fields. Takes the integration weights
   * as input, along with a pointer to the element.
   * The element is assumed to have a constant Jacobian
   */
  virtual void compute_affine_map(const std::vector<Real>& qw,
		   const Elem* e);

  /**
   * Compute the jacobian and some other additional
   * data fields at the single point with index p.
   */
  void compute_single_point_map(const std::vector<Real>& qw,
		                const Elem* e,
				unsigned int p);
  
  /**
   * A utility function for use by compute_*_map
   */
  void resize_map_vectors(unsigned int n_qp);

  /** 
   * Same as compute_map, but for a side.  Useful for boundary integration.
   */  
  void compute_face_map(const std::vector<Real>& qw,
			const Elem* side);

  /** 
   * Same as before, but for an edge.  Useful for some projections.
   */  
  void compute_edge_map(const std::vector<Real>& qw,
			const Elem* side);

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

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the x value of the pth entry of the dxzydxi_map.
   */
  Real dxdxi_map(const unsigned int p) const   { return dxyzdxi_map[p](0); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the y value of the pth entry of the dxzydxi_map.
   */
  Real dydxi_map(const unsigned int p) const   { return dxyzdxi_map[p](1); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the z value of the pth entry of the dxzydxi_map.
   */
  Real dzdxi_map(const unsigned int p) const   { return dxyzdxi_map[p](2); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the x value of the pth entry of the dxzydeta_map.
   */
  Real dxdeta_map(const unsigned int p) const  { return dxyzdeta_map[p](0); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the y value of the pth entry of the dxzydeta_map.
   */
  Real dydeta_map(const unsigned int p) const  { return dxyzdeta_map[p](1); } 

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the z value of the pth entry of the dxzydeta_map.
   */
  Real dzdeta_map(const unsigned int p) const  { return dxyzdeta_map[p](2); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the x value of the pth entry of the dxzydzeta_map.
   */
  Real dxdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](0); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the y value of the pth entry of the dxzydzeta_map.
   */
  Real dydzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](1); }

  /**
   * Used in \p FEBase::compute_map(), which should be
   * be usable in derived classes, and therefore protected.
   * Returns the z value of the pth entry of the dxzydzeta_map.
   */
  Real dzdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](2); }




  
  /**
   * The dimensionality of the object
   */
  const unsigned int dim;

  /**
   * The spatial locations of the quadrature points
   */
  std::vector<Point> xyz;


  
  /**
   * Vector of parital derivatives: 
   * d(x)/d(xi), d(y)/d(xi), d(z)/d(xi) 
   */
  std::vector<RealGradient> dxyzdxi_map;

  /**
   * Vector of parital derivatives: 
   * d(x)/d(eta), d(y)/d(eta), d(z)/d(eta)
   */
  std::vector<RealGradient> dxyzdeta_map;

  /**
   * Vector of parital derivatives: 
   * d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta)
   */
  std::vector<RealGradient> dxyzdzeta_map;
  
  /**
   * Vector of second partial derivatives in xi:
   * d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2
   */
  std::vector<RealGradient> d2xyzdxi2_map;

  /**
   * Vector of mixed second partial derivatives in xi-eta:
   * d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(xi)d(eta)
   */
  std::vector<RealGradient> d2xyzdxideta_map;

  /**
   * Vector of second partial derivatives in eta:
   * d^2(x)/d(eta)^2 
   */
  std::vector<RealGradient> d2xyzdeta2_map;

#ifdef ENABLE_SECOND_DERIVATIVES

  /**
   * Vector of second partial derivatives in xi-zeta:
   * d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta), d^2(z)/d(xi)d(zeta)
   */
  std::vector<RealGradient> d2xyzdxidzeta_map;

  /**
   * Vector of mixed second partial derivatives in eta-zeta:
   * d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2(z)/d(eta)d(zeta)
   */
  std::vector<RealGradient> d2xyzdetadzeta_map;

  /**
   * Vector of second partial derivatives in zeta:
   * d^2(x)/d(zeta)^2 
   */
  std::vector<RealGradient> d2xyzdzeta2_map;

#endif
  
  /**
   * Map for partial derivatives:
   * d(xi)/d(x). Needed for the Jacobian.
   */
  std::vector<Real>  dxidx_map;

  /**
   * Map for partial derivatives:
   * d(xi)/d(y). Needed for the Jacobian.
   */
  std::vector<Real>  dxidy_map;

  /**
   * Map for partial derivatives: