quadrature_monomial.c

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

// $Id: quadrature_monomial.C 2932 2008-07-14 22:29:43Z 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


#include "quadrature_monomial.h"


// See the files:

// quadrature_monomial_2D.C
// quadrature_monomial_3D.C

// for implementation of specific element types.

// Constructor
QMonomial::QMonomial(const unsigned int d,
		     const Order o) : QBase(d,o)
{
}


// Destructor
QMonomial::~QMonomial()
{
}

void QMonomial::wissmann_rule(const Real rule_data[][3],
			      const unsigned int n_pts)
{
  for (unsigned int i=0, c=0; i<n_pts; ++i)
    {
      _points[c]  = Point( rule_data[i][0], rule_data[i][1] );
      _weights[c++] = rule_data[i][2];

      // This may be an (x1,x2) -> (-x1,x2) point, in which case
      // we will also generate the mirror point using the same weight.
      if (rule_data[i][0] != static_cast<Real>(0.0))
	{
	  _points[c]  = Point( -rule_data[i][0], rule_data[i][1] );
	  _weights[c++] = rule_data[i][2];
	}
    }
}



void QMonomial::stroud_rule(const Real rule_data[][3],
			    const unsigned int* rule_symmetry,
			    const unsigned int n_pts)
{
  for (unsigned int i=0, c=0; i<n_pts; ++i)
    {
      const Real
	x=rule_data[i][0],
	y=rule_data[i][1],
	w=rule_data[i][2];

      switch(rule_symmetry[i])
	{
	case 0: // Single point (no symmetry)
	  {
	    _points[c]  = Point( x, y);
	    _weights[c++] = w;
	    
	    break;
	  }
	case 1: // Fully-symmetric (x,y)
	  {
	    _points[c]    = Point( x, y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x, y);
	    _weights[c++] = w;

	    _points[c]    = Point( x,-y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x,-y);
	    _weights[c++] = w;
	    
	    _points[c]    = Point( y, x);
	    _weights[c++] = w;

	    _points[c]    = Point(-y, x);
	    _weights[c++] = w;

	    _points[c]    = Point( y,-x);
	    _weights[c++] = w;

	    _points[c]    = Point(-y,-x);
	    _weights[c++] = w;
	    
	    break;
	  }
	case 2: // Fully-symmetric (x,x)
	  {
	    _points[c]    = Point( x, x);
	    _weights[c++] = w;

	    _points[c]    = Point(-x, x);
	    _weights[c++] = w;

	    _points[c]    = Point( x,-x);
	    _weights[c++] = w;

	    _points[c]    = Point(-x,-x);
	    _weights[c++] = w;

	    break;
	  }
	case 3: // Fully-symmetric (x,0)
	  {
	    libmesh_assert(y==0.0);

	    _points[c]    = Point( x,0.);
	    _weights[c++] = w;

	    _points[c]    = Point(-x,0.);
	    _weights[c++] = w;

	    _points[c]    = Point(0., x);
	    _weights[c++] = w;

	    _points[c]    = Point(0.,-x);
	    _weights[c++] = w;

	    break;
	  }
	case 4: // Rotational invariant
	  {
	    _points[c]    = Point( x, y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x,-y);
	    _weights[c++] = w;

	    _points[c]    = Point(-y, x);
	    _weights[c++] = w;

	    _points[c]    = Point( y,-x);
	    _weights[c++] = w;

	    break;
	  }
	case 5: // Partial symmetry (Wissman's rules)
	  {
	    libmesh_assert (x != 0.0);
	    
	    _points[c]    = Point( x, y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x, y);
	    _weights[c++] = w;

	    break;
	  }
	case 6: // Rectangular symmetry
	  {
	    _points[c]    = Point( x, y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x, y);
	    _weights[c++] = w;

	    _points[c]    = Point(-x,-y);
	    _weights[c++] = w;

	    _points[c]    = Point( x,-y);
	    _weights[c++] = w;

	    break;
	  }
	case 7: // Central symmetry
	  {
	    libmesh_assert (x == 0.0);
	    libmesh_assert (y != 0.0);
	    
	    _points[c]    = Point(0., y);
	    _weights[c++] = w;

	    _points[c]    = Point(0.,-y);
	    _weights[c++] = w;

	    break;
	  }
	default:
	  {
	    std::cerr << "Unknown symmetry!" << std::endl;
	    libmesh_error();
	  }
	} // end switch(rule_symmetry[i])
    }
}