fe_xyz.c

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

// $Id: fe_xyz.C 2856 2008-06-05 15:41:49Z 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 "libmesh_logging.h"
#include "fe.h"
#include "fe_macro.h"
#include "elem.h"




// ------------------------------------------------------------
// XYZ-specific implementations
template <unsigned int Dim>
void FEXYZ<Dim>::init_shape_functions(const std::vector<Point>& qp,
				      const Elem* elem)
{
  libmesh_assert (elem  != NULL);
  this->calculations_started = true;

  // If the user forgot to request anything, we'll be safe and
  // calculate everything:
#ifdef ENABLE_SECOND_DERIVATIVES
  if (!this->calculate_phi && !this->calculate_dphi && !this->calculate_d2phi)
    this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = true;
#else
  if (!this->calculate_phi && !this->calculate_dphi)
    this->calculate_phi = this->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
  {
    // (note: GCC 3.4.0 requires the use of this-> here)
    if (this->calculate_phi)
      this->phi.resize     (n_approx_shape_functions);
    if (this->calculate_dphi)
      {
        this->dphi.resize    (n_approx_shape_functions);
        this->dphidx.resize  (n_approx_shape_functions);
        this->dphidy.resize  (n_approx_shape_functions);
        this->dphidz.resize  (n_approx_shape_functions);
      }
          
      if (Dim > 1)
        this->dphideta_map.resize  (n_mapping_shape_functions);
    
      if (Dim == 3)
        this->dphidzeta_map.resize (n_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
    if (this->calculate_d2phi)
      {
        this->d2phi.resize     (n_approx_shape_functions);
        this->d2phidx2.resize  (n_approx_shape_functions);
        this->d2phidxdy.resize (n_approx_shape_functions);
        this->d2phidxdz.resize (n_approx_shape_functions);
        this->d2phidy2.resize  (n_approx_shape_functions);
        this->d2phidydz.resize (n_approx_shape_functions);
        this->d2phidz2.resize  (n_approx_shape_functions);
        this->d2phidxi2.resize (n_approx_shape_functions);
      }
      if (Dim > 1)
        {
          this->d2phidxideta_map.resize (n_mapping_shape_functions);
          this->d2phideta2_map.resize   (n_mapping_shape_functions);
        }
      if (Dim > 2)
        {
          this->d2phidxidzeta_map.resize  (n_mapping_shape_functions);
          this->d2phidetadzeta_map.resize (n_mapping_shape_functions);
          this->d2phidzeta2_map.resize    (n_mapping_shape_functions);
        }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    this->phi_map.resize         (n_mapping_shape_functions);
    this->dphidxi_map.resize     (n_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
    this->d2phidxi2_map.resize   (n_mapping_shape_functions);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    for (unsigned int i=0; i<n_approx_shape_functions; i++)
      {
        if (this->calculate_phi)
	  this->phi[i].resize           (n_qp);
        if (this->calculate_dphi)
          {
	    this->dphi[i].resize        (n_qp);
	    this->dphidx[i].resize      (n_qp);
	    this->dphidy[i].resize      (n_qp);
	    this->dphidz[i].resize      (n_qp);
          }
#ifdef ENABLE_SECOND_DERIVATIVES
        if (this->calculate_d2phi)
          {
	    this->d2phi[i].resize       (n_qp);
	    this->d2phidx2[i].resize    (n_qp);
	    this->d2phidxdy[i].resize   (n_qp);
	    this->d2phidy2[i].resize    (n_qp);
	    this->d2phidydz[i].resize   (n_qp);
	    this->d2phidz2[i].resize    (n_qp);
          }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
      }
       
    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
      {
	this->phi_map[i].resize         (n_qp);
	this->dphidxi_map[i].resize     (n_qp);
#ifdef ENABLE_SECOND_DERIVATIVES
	this->d2phidxi2_map[i].resize   (n_qp);
        if (Dim > 1)
          {
	    this->d2phidxideta_map[i].resize   (n_qp);
	    this->d2phideta2_map[i].resize     (n_qp);
          }
	if (Dim > 2)
          {
	    this->d2phidxidzeta_map[i].resize  (n_qp);
	    this->d2phidetadzeta_map[i].resize (n_qp);
	    this->d2phidzeta2_map[i].resize    (n_qp);
          }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	   
	if (Dim > 1)
	  this->dphideta_map[i].resize  (n_qp);
	   
	if (Dim == 3)
	  this->dphidzeta_map[i].resize (n_qp);
      }
  }


      
#ifdef ENABLE_INFINITE_ELEMENTS
  //------------------------------------------------------------
  // Initialize the data fields, which should only be used for infinite 
  // elements, to some sensible values, so that using a FE with the
  // variational formulation of an InfFE, correct element matrices are
  // returned

 {
    this->weight.resize  (n_qp);
    this->dweight.resize (n_qp);
    this->dphase.resize  (n_qp);
    
    for (unsigned int p=0; p<n_qp; p++)
      {
        this->weight[p] = 1.;
	this->dweight[p].zero();
	this->dphase[p].zero();
      }

 }
#endif // ifdef ENABLE_INFINITE_ELEMENTS

  // Optimize for the affine elements case:
  bool has_affine_map = elem->has_affine_map();
  
  switch (Dim)
    {

      //------------------------------------------------------------
      // 1D
    case 1:
      {
        // Compute the value of the mapping shape function i at quadrature point p
        // (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
            for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
                this->phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
                this->dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                for (unsigned int p=1; p<n_qp; p++)
                  {
                    this->phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
                    this->dphidxi_map[i][p]  = this->dphidxi_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
          for (unsigned int i=0; i<n_mapping_shape_functions; i++)
            for (unsigned int p=0; p<n_qp; p++)
              {
                this->phi_map[i][p]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
                this->dphidxi_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
              }

        break;

      }


      
      //------------------------------------------------------------
      // 2D
    case 2:
      {
 	// Compute the value of the mapping shape function i at quadrature point p
	// (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
	    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
	        this->phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
	        this->dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
	        this->dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
                this->d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
                this->d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	        for (unsigned int p=1; p<n_qp; p++)
                  {
	            this->phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
	            this->dphidxi_map[i][p]  = this->dphidxi_map[i][0];
	            this->dphideta_map[i][p] = this->dphideta_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
                    this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
                    this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
	  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        this->phi_map[i][p]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
	        this->dphidxi_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
	        this->dphideta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
                this->d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
                this->d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	      }
			
       	break;
      }


      
      //------------------------------------------------------------
      // 3D
    case 3:
      {
	// Compute the value of the mapping shape function i at quadrature point p
	// (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
	    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
	        this->phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
	        this->dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
	        this->dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
	        this->dphidzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
                this->d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
                this->d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
                this->d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[0]);
                this->d2phidetadzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[0]);
                this->d2phidzeta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	        for (unsigned int p=1; p<n_qp; p++)
                  {
	            this->phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
	            this->dphidxi_map[i][p]  = this->dphidxi_map[i][0];
	            this->dphideta_map[i][p] = this->dphideta_map[i][0];
	            this->dphidzeta_map[i][p] = this->dphidzeta_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
                    this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
                    this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
                    this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
                    this->d2phidetadzeta_map[i][p] = this->d2phidetadzeta_map[i][0];
                    this->d2phidzeta2_map[i][p] = this->d2phidzeta2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
	  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        this->phi_map[i][p]       = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
	        this->dphidxi_map[i][p]   = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
	        this->dphideta_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
	        this->dphidzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
                this->d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
                this->d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
                this->d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[p]);
                this->d2phidetadzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[p]);
                this->d2phidzeta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	      }
			
	break;
 
      }


    default:
      libmesh_error();
    }
  
  // Stop logging the shape function initialization
  STOP_LOG("init_shape_functions()", "FE");
}




template <unsigned int Dim>
void FEXYZ<Dim>::compute_shape_functions (const Elem* elem)
{
  libmesh_assert (elem != NULL);