inf_fe_static.c

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

// $Id: inf_fe_static.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 "inf_fe_macro.h"
#include "fe.h"
#include "fe_interface.h"
#include "fe_compute_data.h"
#include "elem.h"


// ------------------------------------------------------------
// InfFE class static member initialization
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
ElemType InfFE<Dim,T_radial,T_map>::_compute_node_indices_fast_current_elem_type = INVALID_ELEM;

#ifdef DEBUG

template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
bool InfFE<Dim,T_radial,T_map>::_warned_for_nodal_soln = false;


template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
bool InfFE<Dim,T_radial,T_map>::_warned_for_shape      = false;

#endif




// ------------------------------------------------------------
// InfFE static class members
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
unsigned int InfFE<Dim,T_radial,T_map>::n_dofs (const FEType& fet,
						const ElemType inf_elem_type)
{
  const ElemType base_et (Base::get_elem_type(inf_elem_type));
    
  if (Dim > 1)
    return FEInterface::n_dofs(Dim-1, fet, base_et) * Radial::n_dofs(fet.radial_order);
  else
    return Radial::n_dofs(fet.radial_order);
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType& fet,
							const ElemType inf_elem_type,
							const unsigned int n)
{
  const ElemType base_et (Base::get_elem_type(inf_elem_type));

  unsigned int n_base, n_radial;
  compute_node_indices(inf_elem_type, n, n_base, n_radial);
  
//   std::cout << "elem_type=" << inf_elem_type 
// 	    << ",  fet.radial_order=" << fet.radial_order
// 	    << ",  n=" << n 
// 	    << ",  n_radial=" << n_radial 
// 	    << ",  n_base=" << n_base
// 	    << std::endl;

  if (Dim > 1)
    return FEInterface::n_dofs_at_node(Dim-1, fet, base_et, n_base) 
        * Radial::n_dofs_at_node(fet.radial_order, n_radial);
  else
    return Radial::n_dofs_at_node(fet.radial_order, n_radial);
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_per_elem (const FEType& fet,
							 const ElemType inf_elem_type)
{
  const ElemType base_et (Base::get_elem_type(inf_elem_type));

  if (Dim > 1)
    return FEInterface::n_dofs_per_elem(Dim-1, fet, base_et) 
        * Radial::n_dofs_per_elem(fet.radial_order);
  else
    return Radial::n_dofs_per_elem(fet.radial_order);
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::nodal_soln(const FEType& /* fet */,
					   const Elem* /* elem */,
					   const std::vector<Number>& /* elem_soln */,
					   std::vector<Number>&       nodal_soln)
{
#ifdef DEBUG
  if (!_warned_for_nodal_soln)
    {
      std::cerr << "WARNING: nodal_soln(...) does _not_ work for infinite elements." << std::endl
		<< " Will return an empty nodal solution.  Use " << std::endl
		<< " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" << std::endl;
      _warned_for_nodal_soln = true;
    }
#endif

  /*
   * In the base the infinite element couples to
   * conventional finite elements.  To not destroy
   * the values there, clear \p nodal_soln.  This
   * indicates to the user of \p nodal_soln to
   * not use this result.
   */
  nodal_soln.clear();
  libmesh_assert (nodal_soln.begin() == nodal_soln.end());
  return;
}








template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
Real InfFE<Dim,T_radial,T_map>::shape(const FEType& fet,
				      const ElemType inf_elem_type,
				      const unsigned int i,
				      const Point& p)
{
  libmesh_assert (Dim != 0);

#ifdef DEBUG
  // this makes only sense when used for mapping
  if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
    {
      std::cerr << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
		<< " return the correct trial function!  Use " << std::endl
		<< " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" 
		<< std::endl;
      _warned_for_shape = true;
    }
#endif

  const ElemType     base_et  (Base::get_elem_type(inf_elem_type));
  const Order        o_radial (fet.radial_order);
  const Real         v        (p(Dim-1));

  unsigned int i_base, i_radial;
  compute_shape_indices(fet, inf_elem_type, i, i_base, i_radial);

  //TODO:[SP/DD]  exp(ikr) is still missing here!
  if (Dim > 1)
    return FEInterface::shape(Dim-1, fet, base_et, i_base, p)
        * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
        * InfFE<Dim,T_radial,T_map>::Radial::decay(v);
  else
    return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
        * InfFE<Dim,T_radial,T_map>::Radial::decay(v);
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
Real InfFE<Dim,T_radial,T_map>::shape(const FEType& fet,
				      const Elem* inf_elem,
				      const unsigned int i,
				      const Point& p)
{
  libmesh_assert (inf_elem != NULL);
  libmesh_assert (Dim != 0);

#ifdef DEBUG
  // this makes only sense when used for mapping
  if ((T_radial != INFINITE_MAP) && !_warned_for_shape)
    {
      std::cerr << "WARNING: InfFE<Dim,T_radial,T_map>::shape(...) does _not_" << std::endl
		<< " return the correct trial function!  Use " << std::endl
		<< " InfFE<Dim,T_radial,T_map>::compute_data(..) instead!" 
		<< std::endl;
      _warned_for_shape = true;
    }
#endif

  const Order        o_radial (fet.radial_order);
  const Real         v        (p(Dim-1));
  AutoPtr<Elem>      base_el  (inf_elem->build_side(0));

  unsigned int i_base, i_radial;
  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);

  if (Dim > 1)
    return FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p)
        * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
        * InfFE<Dim,T_radial,T_map>::Radial::decay(v);
  else
    return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
        * InfFE<Dim,T_radial,T_map>::Radial::decay(v);
}







template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::compute_data(const FEType& fet,
					     const Elem* inf_elem,
					     FEComputeData& data)
{
  libmesh_assert (inf_elem != NULL);
  libmesh_assert (Dim != 0);

  const Order        o_radial             (fet.radial_order);
  const Order        radial_mapping_order (Radial::mapping_order());    
  const Point&       p                    (data.p);
  const Real         v                    (p(Dim-1));
  AutoPtr<Elem>      base_el              (inf_elem->build_side(0));

  /*
   * compute \p interpolated_dist containing the mapping-interpolated 
   * distance of the base point to the origin.  This is the same
   * for all shape functions.  Set \p interpolated_dist to 0, it
   * is added to.
   */
  Real interpolated_dist = 0.;
  switch (Dim)
    {  
    case 1:
      {	
        libmesh_assert (inf_elem->type() == INFEDGE2);
	interpolated_dist =  Point(inf_elem->point(0) - inf_elem->point(1)).size();
	break;
      }

    case 2:
      {
	const unsigned int n_base_nodes = base_el->n_nodes();

	const Point    origin                 = inf_elem->origin();
	const Order    base_mapping_order     (base_el->default_order());
	const ElemType base_mapping_elem_type (base_el->type());

	// interpolate the base nodes' distances
	for (unsigned int n=0; n<n_base_nodes; n++)
	    interpolated_dist += Point(base_el->point(n) - origin).size()
		* FE<1,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
	break;
      }

    case 3:
      {
	const unsigned int n_base_nodes = base_el->n_nodes();

	const Point    origin                 = inf_elem->origin();
	const Order    base_mapping_order     (base_el->default_order());
	const ElemType base_mapping_elem_type (base_el->type());

	// interpolate the base nodes' distances
	for (unsigned int n=0; n<n_base_nodes; n++)
	    interpolated_dist += Point(base_el->point(n) - origin).size()
		* FE<2,LAGRANGE>::shape (base_mapping_elem_type, base_mapping_order, n, p);
	break;
      }
#ifdef DEBUG
    default:
	libmesh_error();
#endif
    }



#ifdef USE_COMPLEX_NUMBERS

  // assumption on time-harmonic behavior
  const short int sign (-1);

  // the wave number
  const Real wavenumber = 2. * libMesh::pi * data.frequency / data.speed;

  // the exponent for time-harmonic behavior
  const Real exponent = sign                                                            /* +1. or -1.                */
      * wavenumber                                                                      /* k                         */
      * interpolated_dist                                                               /* together with next line:  */
      * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1);                /* phase(s,t,v)              */

  const Number time_harmonic = Number(cos(exponent), sin(exponent));                    /* e^(sign*i*k*phase(s,t,v)) */

  /*
   * compute \p shape for all dof in the element
   */
  if (Dim > 1)
    {
      const unsigned int n_dof = n_dofs (fet, inf_elem->type());
      data.shape.resize(n_dof);

      for (unsigned int i=0; i<n_dof; i++)
        {
	  // compute base and radial shape indices
	  unsigned int i_base, i_radial;
	  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
 
	  data.shape[i] = (InfFE<Dim,T_radial,T_map>::Radial::decay(v)                  /* (1.-v)/2. in 3D          */
			   *  FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p)  /* S_n(s,t)                 */
			   * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial))    /* L_n(v)                   */
	      * time_harmonic;                                                          /* e^(sign*i*k*phase(s,t,v) */
	}
    }
  else
    {	
      std::cerr << "compute_data() for 1-dimensional InfFE not implemented." << std::endl;
      libmesh_error();
    }

#else

  const Real speed = data.speed;

  /*
   * This is quite weird: the phase is actually 
   * a measure how @e advanced the pressure is that
   * we compute.  In other words: the further away 
   * the node \p data.p is, the further we look into
   * the future... 
   */
  data.phase = interpolated_dist                                                       /* phase(s,t,v)/c  */
      * InfFE<Dim,INFINITE_MAP,T_map>::eval(v, radial_mapping_order, 1) / speed;

  if (Dim > 1)
    {
      const unsigned int n_dof = n_dofs (fet, inf_elem->type());
      data.shape.resize(n_dof);

      for (unsigned int i=0; i<n_dof; i++)
        {
	  // compute base and radial shape indices
	  unsigned int i_base, i_radial;
	  compute_shape_indices(fet, inf_elem->type(), i, i_base, i_radial);
 
	  data.shape[i] = InfFE<Dim,T_radial,T_map>::Radial::decay(v)                  /* (1.-v)/2. in 3D */
	                  *  FEInterface::shape(Dim-1, fet, base_el.get(), i_base, p)  /* S_n(s,t)        */
	                  * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial);    /* L_n(v)          */
	}
    }
  else
    {	
      std::cerr << "compute_data() for 1-dimensional InfFE not implemented." << std::endl;
      libmesh_error();
    }

#endif
}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::compute_node_indices (const ElemType inf_elem_type,
						      const unsigned int outer_node_index,
						      unsigned int& base_node,
						      unsigned int& radial_node)
{
  switch (inf_elem_type)
    {
    case INFEDGE2:
      {
	libmesh_assert (outer_node_index < 2);
	base_node   = 0;
	radial_node = outer_node_index;
	return;
      }


    // linear base approximation, easy to determine
    case INFQUAD4:
      {
	libmesh_assert (outer_node_index < 4);
	base_node   = outer_node_index % 2;
	radial_node = outer_node_index / 2;
	return;
      }

    case INFPRISM6:
      {
	libmesh_assert (outer_node_index < 6);
	base_node   = outer_node_index % 3;
	radial_node = outer_node_index / 3;
	return;
      }

    case INFHEX8:
      {
	libmesh_assert (outer_node_index < 8);
	base_node   = outer_node_index % 4;
	radial_node = outer_node_index / 4;
	return;
      }


    // higher order base approximation, more work necessary
    case INFQUAD6:
      {
	switch (outer_node_index)
	  {
	  case 0:
	  case 1:
	    {
	      radial_node = 0;
	      base_node   = outer_node_index;
	      return;
	    }

	  case 2:
	  case 3:
	    {
	      radial_node = 1;
	      base_node   = outer_node_index-2;
	      return;
	    }

	  case 4:
	    {
	      radial_node = 0;
	      base_node   = 2;
	      return;
	    }

	  case 5:
	    {
	      radial_node = 1;
	      base_node   = 2;
	      return;
	    }

	  default:
	    {
	      libmesh_error();
	      return;
	    }
	  }
      }


    case INFHEX16:
    case INFHEX18:
      {
	switch (outer_node_index)
	  {
	  case 0:
	  case 1:
	  case 2:
	  case 3:
	    {
	      radial_node = 0;
	      base_node   = outer_node_index;
	      return;
	    }

	  case 4:
	  case 5:
	  case 6:
	  case 7:
	    {
	      radial_node = 1;
	      base_node   = outer_node_index-4;