inf_fe.c

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

// $Id: inf_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 "libmesh_config.h"
#ifdef ENABLE_INFINITE_ELEMENTS
#include "inf_fe.h"
#include "quadrature_gauss.h"
#include "elem.h"
#include "libmesh_logging.h"



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



// Constructor
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
InfFE<Dim,T_radial,T_map>::InfFE (const FEType& fet) :
  FEBase       (Dim, fet),

  _n_total_approx_sf (0),
  _n_total_qp        (0),

  base_qrule   (NULL),
  radial_qrule (NULL),
  base_elem    (NULL),
  base_fe      (NULL),
  
  // initialize the current_fe_type to all the same
  // values as \p fet (since the FE families and coordinate
  // map type should @e not change), but use an invalid order
  // for the radial part (since this is the only order
  // that may change!).
  // the data structures like \p phi etc are not initialized 
  // through the constructor, but throught reinit()
  current_fe_type ( FEType(fet.order, 
			   fet.family, 
			   INVALID_ORDER, 
			   fet.radial_family,      
			   fet.inf_map) )

{
  // Sanity checks
  libmesh_assert (T_radial == fe_type.radial_family);
  libmesh_assert (T_map    == fe_type.inf_map);

  // build the base_fe object, handle the AutoPtr
  if (Dim != 1)
    {
      AutoPtr<FEBase> ap_fb(FEBase::build(Dim-1, fet));
      base_fe = ap_fb.release();
    }
}




// Destructor
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
InfFE<Dim,T_radial,T_map>::~InfFE ()
{
  // delete pointers, if necessary
  if (base_qrule != NULL)
    {
      delete base_qrule;
      base_qrule = NULL;
    }
  
  if (radial_qrule != NULL)
    {
      delete radial_qrule;
      radial_qrule = NULL;
    }

  if (base_elem != NULL)
    {
      delete base_elem;
      base_elem = NULL;
    }
  
  if (base_fe != NULL)
    {
      delete base_fe;
      base_fe = NULL;
    }
}





template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>:: attach_quadrature_rule (QBase* q)
{
  libmesh_assert (q       != NULL);
  libmesh_assert (base_fe != NULL);

  const Order base_int_order   = q->get_order();
  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order) + 1) +2);
  const unsigned int qrule_dim = q->get_dim();

  if (Dim != 1)
    {
      // build a Dim-1 quadrature rule of the type that we received
      AutoPtr<QBase> apq( QBase::build(q->type(), qrule_dim-1, base_int_order) );
      base_qrule = apq.release();
      base_fe->attach_quadrature_rule(base_qrule);
    }
  
  // in radial direction, always use Gauss quadrature
  radial_qrule = new QGauss(1, radial_int_order);

  // currently not used. But maybe helpful to store the QBase*
  // with which we initialized our own quadrature rules
  qrule = q;
}





template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
void InfFE<Dim,T_radial,T_base>::update_base_elem (const Elem* inf_elem)
{
  if (base_elem != NULL)
    delete base_elem;
  base_elem = Base::build_elem(inf_elem);
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::reinit(const Elem* inf_elem,
				       const std::vector<Point>* const pts)
{
  libmesh_assert (base_fe        != NULL);
  libmesh_assert (base_fe->qrule != NULL);
  libmesh_assert (base_fe->qrule == base_qrule);
  libmesh_assert (radial_qrule   != NULL);
  libmesh_assert (inf_elem       != NULL);

  if (pts == NULL)
    {
      bool init_shape_functions_required = false;
      
      // -----------------------------------------------------------------
      // init the radial data fields only when the radial order changes
      if (current_fe_type.radial_order != fe_type.radial_order)
	{
	  current_fe_type.radial_order = fe_type.radial_order;
	  
	  // Watch out: this call to QBase->init() only works for 
	  // current_fe_type = const!   To allow variable Order, 
	  // the init() of QBase has to be modified...  
	  radial_qrule->init(EDGE2);
	  
	  // initialize the radial shape functions
	  this->init_radial_shape_functions(inf_elem);
	  
	  init_shape_functions_required=true;      
	}
      
      
      bool update_base_elem_required=true;
      
      // -----------------------------------------------------------------
      // 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
      // depend on the particular element and need a reinit
      if (  ( Dim != 1) &&
	    (  (this->get_type() != inf_elem->type())  ||  
	       (base_fe->shapes_need_reinit())  )  )
	{
	  // store the new element type, update base_elem
	  // here.  Through \p update_base_elem_required,
	  // remember whether it has to be updated (see below).
	  elem_type = inf_elem->type();
	  this->update_base_elem(inf_elem);
	  update_base_elem_required=false;
	  
	  // initialize the base quadrature rule for the new element
	  base_qrule->init(base_elem->type());
	  
	  // initialize the shape functions in the base
	  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
					     base_elem);
	  
	  init_shape_functions_required=true;      
	}
      
      
      // when either the radial or base part change, 
      // we have to init the whole fields
      if (init_shape_functions_required)
	this->init_shape_functions (inf_elem);
      
      // computing the distance only works when we have the current
      // base_elem stored.  This happens when fe_type is const,
      // the inf_elem->type remains the same.  Then we have to
      // update the base elem _here_.
      if (update_base_elem_required)
	this->update_base_elem(inf_elem);
      
      // compute dist (depends on geometry, therefore has to be updated for
      // each and every new element), throw radial and base part together
      this->combine_base_radial (inf_elem);
      
      this->compute_map (_total_qrule_weights, inf_elem);

      // Compute the shape functions and the derivatives 
      // at all quadrature points.
      this->compute_shape_functions (inf_elem);
    }
  
  else // if pts != NULL
    {
      // update the elem_type
      elem_type = inf_elem->type();

      // init radial shapes
      this->init_radial_shape_functions(inf_elem);

      // update the base
      this->update_base_elem(inf_elem);

      // the finite element on the ifem base
      {
	AutoPtr<FEBase> ap_fb(FEBase::build(Dim-1, this->fe_type));
	if (base_fe != NULL)
	  delete base_fe;
	base_fe = ap_fb.release();
      }

      // inite base shapes
      base_fe->init_base_shape_functions(*pts,
					 base_elem);

      this->init_shape_functions (inf_elem);

      // combine the base and radial shapes
      this->combine_base_radial (inf_elem);

      // dummy weights
      std::vector<Real> dummy_weights (pts->size(), 1.);

      this->compute_map (dummy_weights, inf_elem);

      // finally compute the ifem shapes     
      this->compute_shape_functions (inf_elem);
    }
  
}





template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::init_radial_shape_functions(const Elem* inf_elem)
{
  libmesh_assert (radial_qrule != NULL);
  libmesh_assert (inf_elem     != NULL);

 
  /**
   * Start logging the radial shape function initialization
   */
  START_LOG("init_radial_shape_functions()", "InfFE");


  // -----------------------------------------------------------------
  // initialize most of the things related to mapping

  // The order to use in the radial map (currently independent of the element type)
  const Order        radial_mapping_order             (Radial::mapping_order());    
  const unsigned int n_radial_mapping_shape_functions (Radial::n_dofs(radial_mapping_order));



  // -----------------------------------------------------------------
  // initialize most of the things related to physical approximation
  
  const Order        radial_approx_order             (fe_type.radial_order);
  const unsigned int n_radial_approx_shape_functions (Radial::n_dofs(radial_approx_order));

  const unsigned int        n_radial_qp = radial_qrule->n_points();
  const std::vector<Point>&   radial_qp = radial_qrule->get_points();



  // -----------------------------------------------------------------
  // resize the radial data fields
  
  mode.resize      (n_radial_approx_shape_functions);       // the radial polynomials (eval)
  dmodedv.resize   (n_radial_approx_shape_functions);

  som.resize       (n_radial_qp);                           // the (1-v)/2 weight
  dsomdv.resize    (n_radial_qp);

  radial_map.resize    (n_radial_mapping_shape_functions);  // the radial map
  dradialdv_map.resize (n_radial_mapping_shape_functions);

  
  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
    {
      radial_map[i].resize    (n_radial_qp);
      dradialdv_map[i].resize (n_radial_qp);
    }


  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
    {
      mode[i].resize    (n_radial_qp);
      dmodedv[i].resize (n_radial_qp);
    }


  // compute scalar values at radial quadrature points
  for (unsigned int p=0; p<n_radial_qp; p++)
    {
      som[p]       = Radial::decay       (radial_qp[p](0)); 
      dsomdv[p]    = Radial::decay_deriv (radial_qp[p](0)); 
    }

  
  // evaluate the mode shapes in radial direction at radial quadrature points
  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
    for (unsigned int p=0; p<n_radial_qp; p++)
      {
        mode[i][p]    = InfFE<Dim,T_radial,T_map>::eval       (radial_qp[p](0), radial_approx_order, i);
        dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
      }


  // evaluate the mapping functions in radial direction at radial quadrature points
  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
    for (unsigned int p=0; p<n_radial_qp; p++)
      {
	radial_map[i][p]    = InfFE<Dim,INFINITE_MAP,T_map>::eval       (radial_qp[p](0), radial_mapping_order, i);
	dradialdv_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (radial_qp[p](0), radial_mapping_order, i);
      }

  /**
   * Stop logging the radial shape function initialization
   */
  STOP_LOG("init_radial_shape_functions()", "InfFE");

}





template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::init_shape_functions(const Elem* inf_elem)
{
  libmesh_assert (inf_elem     != NULL);
 
  
  // Start logging the radial shape function initialization
  START_LOG("init_shape_functions()", "InfFE");

  
  // -----------------------------------------------------------------
  // fast access to some const int's for the radial data
  const unsigned int n_radial_mapping_sf = radial_map.size();
  const unsigned int n_radial_approx_sf  = mode.size();
  const unsigned int n_radial_qp         = som.size();


  // -----------------------------------------------------------------
  // initialize most of the things related to mapping
  
  // The element type and order to use in the base map  
  const Order    base_mapping_order     ( base_elem->default_order() );
  const ElemType base_mapping_elem_type ( base_elem->type()          );

  // the number of base shape functions used to construct the map
  // (Lagrange shape functions are used for mapping in the base)
  unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
									base_mapping_order);

  const unsigned int n_total_mapping_shape_functions = 
    n_radial_mapping_sf * n_base_mapping_shape_functions;



  // -----------------------------------------------------------------
  // initialize most of the things related to physical approximation
  
  unsigned int n_base_approx_shape_functions;
  if (Dim > 1)
    n_base_approx_shape_functions = base_fe->n_shape_functions();
  else
    n_base_approx_shape_functions = 1;


  const unsigned int n_total_approx_shape_functions = 
      n_radial_approx_sf * n_base_approx_shape_functions;

  // update class member field
  _n_total_approx_sf = n_total_approx_shape_functions;


  // The number of the base quadrature points.
  const unsigned int        n_base_qp =  base_qrule->n_points();

  // The total number of quadrature points.
  const unsigned int        n_total_qp =  n_radial_qp * n_base_qp;

  
  // update class member field
  _n_total_qp = n_total_qp;



  // -----------------------------------------------------------------
  // initialize the node and shape numbering maps
  {
    // these vectors work as follows: the i-th entry stores
    // the associated base/radial node number
    _radial_node_index.resize    (n_total_mapping_shape_functions);
    _base_node_index.resize      (n_total_mapping_shape_functions);

    // similar for the shapes: the i-th entry stores
    // the associated base/radial shape number
    _radial_shape_index.resize   (n_total_approx_shape_functions);
    _base_shape_index.resize     (n_total_approx_shape_functions);

    const ElemType inf_elem_type (inf_elem->type());

    // fill the node index map
    for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
      {
        compute_node_indices (inf_elem_type, 
			      n,
			      _base_node_index[n], 
			      _radial_node_index[n]);
	libmesh_assert (_base_node_index[n]   < n_base_mapping_shape_functions);
	libmesh_assert (_radial_node_index[n] < n_radial_mapping_sf);
      }

    // fill the shape index map
    for (unsigned int n=0; n<n_total_approx_shape_functions; n++)
      {
	compute_shape_indices (this->fe_type,
			       inf_elem_type, 
			       n,
			       _base_shape_index[n], 
			       _radial_shape_index[n]);
	libmesh_assert (_base_shape_index[n]   < n_base_approx_shape_functions);
	libmesh_assert (_radial_shape_index[n] < n_radial_approx_sf);
      }
  }
 




  // -----------------------------------------------------------------
  // resize the base data fields
  dist.resize(n_base_mapping_shape_functions);



  // -----------------------------------------------------------------
  // resize the total data fields
  
  // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
  //