fe.c

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

// $Id: fe.C 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



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



// ------------------------------------------------------------
// FE class members
template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_shape_functions () const
{
  return FE<Dim,T>::n_dofs (elem_type,
           static_cast<Order>(fe_type.order + _p_level));
}


template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::attach_quadrature_rule (QBase* q)
{
  libmesh_assert (q != NULL); 
  qrule = q;
  // make sure we don't cache results from a previous quadrature rule
  elem_type = INVALID_ELEM;
  return;
}


template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_quadrature_points () const
{ 
  libmesh_assert (qrule != NULL);  
  return qrule->n_points(); 
}


template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::dofs_on_side(const Elem* const elem,
			     const Order o,
			     unsigned int s,
			     std::vector<unsigned int>& di)
{
  libmesh_assert(elem != NULL);
  libmesh_assert(s < elem->n_sides());

  di.clear();
  unsigned int nodenum = 0;
  const unsigned int n_nodes = elem->n_nodes();
  for (unsigned int n = 0; n != n_nodes; ++n)
    {
      const unsigned int n_dofs = n_dofs_at_node(elem->type(),
						 static_cast<Order>(o + elem->p_level()), n);
      if (elem->is_node_on_side(n, s))
	for (unsigned int i = 0; i != n_dofs; ++i)
	  di.push_back(nodenum++);
      else
	nodenum += n_dofs;
    }
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::dofs_on_edge(const Elem* const elem,
			     const Order o,
			     unsigned int e,
			     std::vector<unsigned int>& di)
{
  libmesh_assert(elem != NULL);
  libmesh_assert(e < elem->n_edges());

  di.clear();
  unsigned int nodenum = 0;
  const unsigned int n_nodes = elem->n_nodes();
  for (unsigned int n = 0; n != n_nodes; ++n)
    {
      const unsigned int n_dofs = n_dofs_at_node(elem->type(),
						 static_cast<Order>(o + elem->p_level()), n);
      if (elem->is_node_on_edge(n, e))
	for (unsigned int i = 0; i != n_dofs; ++i)
	  di.push_back(nodenum++);
      else
	nodenum += n_dofs;
    }
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::reinit(const Elem* elem,
		       const std::vector<Point>* const pts)
{
  libmesh_assert (elem != NULL);

  // Try to avoid calling init_shape_functions
  // even when shapes_need_reinit
  bool cached_nodes_still_fit = false;

  // Initialize the shape functions at the user-specified
  // points
  if (pts != NULL)
    {
      // Set the type and p level for this element
      elem_type = elem->type();
      _p_level = elem->p_level();

      // Initialize the shape functions
      this->init_shape_functions (*pts, elem);

      // The shape functions do not correspond to the qrule
      shapes_on_quadrature = false;
    }
  
  // If there are no user specified points, we use the
  // quadrature rule
  
  // update the type in accordance to the current cell
  // and reinit if the cell type has changed or (as in
  // the case of the hierarchics) the shape functions need
  // reinit, since they depend on the particular element shape
  else 
    {
      libmesh_assert (qrule   != NULL);
      qrule->init(elem->type(), elem->p_level());

      if (elem_type != elem->type() ||
          _p_level != elem->p_level() ||
          !shapes_on_quadrature)
        {
          // Set the type and p level for this element
          elem_type = elem->type();
          _p_level = elem->p_level();
          // Initialize the shape functions
          this->init_shape_functions (qrule->get_points(), elem);

          if (this->shapes_need_reinit())
            {
              cached_nodes.resize(elem->n_nodes());
              for (unsigned int n = 0; n != elem->n_nodes(); ++n)
                {
                  cached_nodes[n] = elem->point(n);
                }
            }
        }
      else
        {
          libmesh_assert(elem->n_nodes() > 1);

          cached_nodes_still_fit = true;
          if (cached_nodes.size() != elem->n_nodes())
            cached_nodes_still_fit = false;
          else
            for (unsigned int n = 1; n < elem->n_nodes(); ++n)
              {
                if ((elem->point(n) - elem->point(0)) !=
                    (cached_nodes[n] - cached_nodes[0]))
                  {
                    cached_nodes_still_fit = false;
                    break;
                  }
              }

          if (this->shapes_need_reinit() && !cached_nodes_still_fit)
            {
              this->init_shape_functions (qrule->get_points(), elem);
              cached_nodes.resize(elem->n_nodes());
              for (unsigned int n = 0; n != elem->n_nodes(); ++n)
                cached_nodes[n] = elem->point(n);
            }
        }
      
      // The shape functions correspond to the qrule
      shapes_on_quadrature = true;
    }
  
  // Compute the map for this element.  In the future we can specify
  // different types of maps
  if (pts != NULL)
    {
      std::vector<Real> dummy_weights (pts->size(), 1.);
      
      this->compute_map (dummy_weights, elem);
    }
  else
    {
      this->compute_map (qrule->get_weights(), elem);
    }

  // Compute the shape functions and the derivatives at all of the
  // quadrature points.
  if (!cached_nodes_still_fit)
    this->compute_shape_functions (elem);
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::init_shape_functions(const std::vector<Point>& qp,
				     const Elem* elem)
{
  libmesh_assert (elem  != NULL);
  calculations_started = true;

  // If the user forgot to request anything, we'll be safe and
  // calculate everything:
#ifdef ENABLE_SECOND_DERIVATIVES
  if (!calculate_phi && !calculate_dphi && !calculate_d2phi)
    calculate_phi = calculate_dphi = calculate_d2phi = true;
#else
  if (!calculate_phi && !calculate_dphi)
    calculate_phi = calculate_dphi = true;
#endif // ENABLE_SECOND_DERIVATIVES
  
  // Start logging the shape function initialization
  START_LOG("init_shape_functions()", "FE");

  
  // The number of quadrature points.
  const unsigned int n_qp = qp.size();

  // The element type and order to use in
  // the map
  const Order    mapping_order     (elem->default_order());
  const ElemType mapping_elem_type (elem->type());
    
  // Number of shape functions in the finite element approximation
  // space.
  const unsigned int n_approx_shape_functions =
    this->n_shape_functions(this->get_type(),
			    this->get_order());

  // Number of shape functions used to construt the map
  // (Lagrange shape functions are used for mapping)
  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
  // Phi are the shape functions used for the FE approximation
  // Phi_map are the shape functions used for the FE mapping
  {
    if (calculate_phi)
      phi.resize     (n_approx_shape_functions);
    if (calculate_dphi)
      {
        dphi.resize    (n_approx_shape_functions);
        dphidx.resize  (n_approx_shape_functions);
        dphidy.resize  (n_approx_shape_functions);
        dphidz.resize  (n_approx_shape_functions);
        dphidxi.resize (n_approx_shape_functions);

        if (Dim > 1)
          dphideta.resize      (n_approx_shape_functions);
    
        if (Dim == 3)
          dphidzeta.resize     (n_approx_shape_functions);
      }
#ifdef ENABLE_SECOND_DERIVATIVES
    if (calculate_d2phi)
      {
        d2phi.resize     (n_approx_shape_functions);
        d2phidx2.resize  (n_approx_shape_functions);
        d2phidxdy.resize (n_approx_shape_functions);
        d2phidxdz.resize (n_approx_shape_functions);
        d2phidy2.resize  (n_approx_shape_functions);
        d2phidydz.resize (n_approx_shape_functions);
        d2phidz2.resize  (n_approx_shape_functions);
        d2phidxi2.resize (n_approx_shape_functions);
        if (Dim > 1)
          {
            d2phidxideta.resize (n_approx_shape_functions);
            d2phideta2.resize   (n_approx_shape_functions);
          }
        if (Dim > 2)
          {
            d2phidxidzeta.resize  (n_approx_shape_functions);
            d2phidetadzeta.resize (n_approx_shape_functions);
            d2phidzeta2.resize    (n_approx_shape_functions);
          }
      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    phi_map.resize         (n_mapping_shape_functions);
    dphidxi_map.resize     (n_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
    d2phidxi2_map.resize   (n_mapping_shape_functions);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    if (Dim > 1)
      {
        dphideta_map.resize  (n_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
        d2phidxideta_map.resize   (n_mapping_shape_functions);
        d2phideta2_map.resize     (n_mapping_shape_functions);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
      }
    
    if (Dim == 3)
      {
        dphidzeta_map.resize (n_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
        d2phidxidzeta_map.resize  (n_mapping_shape_functions);
        d2phidetadzeta_map.resize (n_mapping_shape_functions);
        d2phidzeta2_map.resize    (n_mapping_shape_functions);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
      }
    
    for (unsigned int i=0; i<n_approx_shape_functions; i++)
      {
	if (calculate_phi)
	  phi[i].resize         (n_qp);
	if (calculate_dphi)
	  {
	    dphi[i].resize        (n_qp);
	    dphidx[i].resize      (n_qp);
	    dphidy[i].resize      (n_qp);
	    dphidz[i].resize      (n_qp);
	    dphidxi[i].resize     (n_qp);

	    if (Dim > 1)
	      dphideta[i].resize  (n_qp);