fe_boundary.c

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

// $Id: fe_boundary.C 2951 2008-08-01 16:32:30Z jwpeterson $

// 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



// C++ includes
#include <cmath> // for std::sqrt


// Local includes
#include "libmesh_common.h"
#include "fe.h"
#include "quadrature.h"
#include "elem.h"
#include "fe_macro.h"
#include "libmesh_logging.h"

//-------------------------------------------------------
// Full specializations for useless methods in 0D, 1D
#define REINIT_ERROR(_dim, _type, _func)       \
template <>                                    \
void FE<_dim,_type>::_func(const Elem*,        \
			   const unsigned int, \
                           const Real)         \
{                                              \
  std::cerr << "ERROR: This method makes no sense for low-D elements!" \
	    << std::endl;                      \
  libmesh_error();                                     \
}

REINIT_ERROR(0, CLOUGH, reinit)
REINIT_ERROR(0, HERMITE, reinit)
REINIT_ERROR(0, HIERARCHIC, reinit)
REINIT_ERROR(0, LAGRANGE, reinit)
REINIT_ERROR(0, MONOMIAL, reinit)
REINIT_ERROR(0, XYZ, reinit)
#ifdef ENABLE_HIGHER_ORDER_SHAPES
REINIT_ERROR(0, BERNSTEIN, reinit)
REINIT_ERROR(0, SZABAB, reinit)
#endif

REINIT_ERROR(1, CLOUGH, edge_reinit)
REINIT_ERROR(1, HERMITE, edge_reinit)
REINIT_ERROR(1, HIERARCHIC, edge_reinit)
REINIT_ERROR(1, LAGRANGE, edge_reinit)
REINIT_ERROR(1, XYZ, edge_reinit)
REINIT_ERROR(1, MONOMIAL, edge_reinit)
#ifdef ENABLE_HIGHER_ORDER_SHAPES
REINIT_ERROR(1, BERNSTEIN, edge_reinit)
REINIT_ERROR(1, SZABAB, edge_reinit)
#endif


//-------------------------------------------------------
// Methods for 2D, 3D
template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::reinit(const Elem* elem,
		       const unsigned int s,
		       const Real tolerance)
{
  libmesh_assert (elem  != NULL);
  libmesh_assert (qrule != NULL);
  // We now do this for 1D elements!
  // libmesh_assert (Dim != 1);

  // Build the side of interest 
  const AutoPtr<Elem> side(elem->build_side(s));


  // Find the max p_level to select 
  // the right quadrature rule for side integration
  unsigned int side_p_level = elem->p_level();
  if (elem->neighbor(s) != NULL)
     side_p_level = std::max(side_p_level, elem->neighbor(s)->p_level());
     
  // initialize quadrature rule
  qrule->init(side->type(), side_p_level);

  // FIXME - could this break if the same FE object was used
  // for both volume and face integrals? - RHS
  // We might not need to reinitialize the shape functions
  if ((this->get_type() != elem->type())    ||
      (s != last_side)                      || 
      (this->get_p_level() != side_p_level) ||
      this->shapes_need_reinit())
    {
      // Set the element type
      elem_type = elem->type();

      // Set the last_side
      last_side = s;
      
      // Set the last p level
      _p_level = side_p_level;
      
      // Initialize the face shape functions
      this->init_face_shape_functions (qrule->get_points(),  side.get());
    }
  
  // Compute the Jacobian*Weight on the face for integration
  this->compute_face_map (qrule->get_weights(), side.get());

  // make a copy of the Jacobian for integration
  const std::vector<Real> JxW_int(JxW);

  // Find where the integration points are located on the
  // full element.
  std::vector<Point> qp; this->inverse_map (elem, xyz, qp, tolerance);
  
  // compute the shape function and derivative values
  // at the points qp
  this->reinit  (elem, &qp);
      
  // copy back old data
  JxW = JxW_int;
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::edge_reinit(const Elem* elem,
		            const unsigned int e,
			    const Real tolerance)
{
  libmesh_assert (elem  != NULL);
  libmesh_assert (qrule != NULL);
  // We don't do this for 1D elements!
  libmesh_assert (Dim != 1);

  // Build the side of interest 
  const AutoPtr<Elem> edge(elem->build_edge(e));

  // initialize quadrature rule
  qrule->init(edge->type(), elem->p_level());

  // We might not need to reinitialize the shape functions
  if ((this->get_type() != elem->type()) ||
      (e != last_edge) ||
      this->shapes_need_reinit())
    {
      // Set the element type
      elem_type = elem->type();

      // Set the last_edge
      last_edge = e;
      
      // Initialize the edge shape functions
      this->init_edge_shape_functions (qrule->get_points(), edge.get());
    }
  
  // Compute the Jacobian*Weight on the face for integration
  this->compute_edge_map (qrule->get_weights(), edge.get());

  // make a copy of the Jacobian for integration
  const std::vector<Real> JxW_int(JxW);

  // Find where the integration points are located on the
  // full element.
  std::vector<Point> qp; this->inverse_map (elem, xyz, qp, tolerance);
  
  // compute the shape function and derivative values
  // at the points qp
  this->reinit  (elem, &qp);
      
  // copy back old data
  JxW = JxW_int;
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::init_face_shape_functions(const std::vector<Point>& qp,
					  const Elem* side)
{
  libmesh_assert (side  != NULL);
  
  /**
   * Start logging the shape function initialization
   */
  START_LOG("init_face_shape_functions()", "FE");

  // The element type and order to use in
  // the map
  const Order    mapping_order     (side->default_order()); 
  const ElemType mapping_elem_type (side->type());

  // The number of quadrature points.
  const unsigned int n_qp = qp.size();
	
  const unsigned int n_mapping_shape_functions =
    FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
					 mapping_order);
  
  // resize the vectors to hold current data
  // Psi are the shape functions used for the FE mapping
  psi_map.resize        (n_mapping_shape_functions);

  if (Dim > 1)
    {
      dpsidxi_map.resize    (n_mapping_shape_functions);
      d2psidxi2_map.resize  (n_mapping_shape_functions);
    }
  
  if (Dim == 3)
    {
      dpsideta_map.resize     (n_mapping_shape_functions);
      d2psidxideta_map.resize (n_mapping_shape_functions);
      d2psideta2_map.resize   (n_mapping_shape_functions);
    }
  
  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
    {
      // Allocate space to store the values of the shape functions
      // and their first and second derivatives at the quadrature points.
      psi_map[i].resize        (n_qp);
      if (Dim > 1)
        {
          dpsidxi_map[i].resize    (n_qp);
          d2psidxi2_map[i].resize  (n_qp);
        }
      if (Dim == 3)
	{
	  dpsideta_map[i].resize     (n_qp);
	  d2psidxideta_map[i].resize (n_qp);
	  d2psideta2_map[i].resize   (n_qp);
	}
  
      // Compute the value of shape function i, and its first and
      // second derivatives at quadrature point p
      // (Lagrange shape functions are used for the mapping)
      for (unsigned int p=0; p<n_qp; p++)
	{
	  psi_map[i][p]        = FE<Dim-1,LAGRANGE>::shape             (mapping_elem_type, mapping_order, i,    qp[p]);
          if (Dim > 1)
	    {
	      dpsidxi_map[i][p]    = FE<Dim-1,LAGRANGE>::shape_deriv       (mapping_elem_type, mapping_order, i, 0, qp[p]);
	      d2psidxi2_map[i][p]  = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
	    }
	  // std::cout << "d2psidxi2_map["<<i<<"][p]=" << d2psidxi2_map[i][p] << std::endl;

	  // If we are in 3D, then our sides are 2D faces.
	  // For the second derivatives, we must also compute the cross
	  // derivative d^2() / dxi deta
	  if (Dim == 3)
	    {
	      dpsideta_map[i][p]     = FE<Dim-1,LAGRANGE>::shape_deriv       (mapping_elem_type, mapping_order, i, 1, qp[p]);
	      d2psidxideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 1, qp[p]); 
	      d2psideta2_map[i][p]   = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 2, qp[p]);
	    }
	}
    }

  
  /**
   * Stop logging the shape function initialization
   */
  STOP_LOG("init_face_shape_functions()", "FE");
}


  
template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::init_edge_shape_functions(const std::vector<Point>& qp,
					  const Elem* edge)
{
  libmesh_assert (edge != NULL);
  
  /**
   * Start logging the shape function initialization
   */
  START_LOG("init_edge_shape_functions()", "FE");

  // The element type and order to use in
  // the map
  const Order    mapping_order     (edge->default_order()); 
  const ElemType mapping_elem_type (edge->type());

  // The number of quadrature points.
  const unsigned int n_qp = qp.size();
	
  const unsigned int n_mapping_shape_functions =
    FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
					 mapping_order);
  
  // resize the vectors to hold current data
  // Psi are the shape functions used for the FE mapping
  psi_map.resize        (n_mapping_shape_functions);
  dpsidxi_map.resize    (n_mapping_shape_functions);
  d2psidxi2_map.resize  (n_mapping_shape_functions);
  
  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
    {
      // Allocate space to store the values of the shape functions
      // and their first and second derivatives at the quadrature points.
      psi_map[i].resize        (n_qp);
      dpsidxi_map[i].resize    (n_qp);
      d2psidxi2_map[i].resize  (n_qp);
  
      // Compute the value of shape function i, and its first and
      // second derivatives at quadrature point p
      // (Lagrange shape functions are used for the mapping)
      for (unsigned int p=0; p<n_qp; p++)
	{
	  psi_map[i][p]        = FE<1,LAGRANGE>::shape             (mapping_elem_type, mapping_order, i,    qp[p]);
	  dpsidxi_map[i][p]    = FE<1,LAGRANGE>::shape_deriv       (mapping_elem_type, mapping_order, i, 0, qp[p]);
	  d2psidxi2_map[i][p]  = FE<1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
	}
    }
  
  /**
   * Stop logging the shape function initialization
   */
  STOP_LOG("init_edge_shape_functions()", "FE");
}

  

void FEBase::compute_face_map(const std::vector<Real>& qw,
			      const Elem* side)
{
  libmesh_assert (side  != NULL);

  START_LOG("compute_face_map()", "FE");

  // The number of quadrature points.
  const unsigned int n_qp = qw.size();
  
  
  switch (dim)
    {
    case 1:
      {
	// A 1D finite element, currently assumed to be in 1D space
	// This means the boundary is a "0D finite element", a
        // NODEELEM.

	// Resize the vectors to hold data at the quadrature points
	{  
	  xyz.resize(n_qp);
	  normals.resize(n_qp);

	  JxW.resize(n_qp);
        }

        // If we have no quadrature points, there's nothing else to do
        if (!n_qp)
          break;

        // We need to look back at the full edge to figure out the normal
        // vector
        const Elem *elem = side->parent();
        libmesh_assert (elem);
        if (side->node(0) == elem->node(0))
          normals[0] = Point(-1.);
        else
          {
            libmesh_assert (side->node(0) == elem->node(1));
            normals[0] = Point(1.);
          }