QuadratureFormula.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: QuadratureFormula.hpp,v 1.10 2007/06/18 23:39:21 delpinux Exp $

#ifndef QUADRATURE_FORMULA_HPP
#define QUADRATURE_FORMULA_HPP

#include <TinyVector.hpp>
#include <ThreadStaticBase.hpp>

#include <ErrorHandler.hpp>

class QuadratureFormulaP0Triangle3D
  : public ThreadStaticBase<QuadratureFormulaP0Triangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 1
  };

private:
  TinyVector<numberOfQuadraturePoints,
             TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3, real_t>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints, TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP0Triangle3D()
    : __weight(0)
  {
    __integrationVertices[0] = TinyVector<3,real_t>(1./3.,1./3.,0);
  }
};

class QuadratureFormulaP1Triangle3D
  : public ThreadStaticBase<QuadratureFormulaP1Triangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 3
  };

private:
  TinyVector<numberOfQuadraturePoints,
             TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3, real_t>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints, TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP1Triangle3D()
  {
    this->__setQuadratureVertices();
  }
private:
  void __setQuadratureVertices()
  {
    __weight = 1./(2.*numberOfQuadraturePoints);

    __integrationVertices[0][0] = 0;
    __integrationVertices[0][1] = 0.5;
    __integrationVertices[0][2] = 0;

    __integrationVertices[1][0] = 0.5;
    __integrationVertices[1][1] = 0;
    __integrationVertices[1][2] = 0;

    __integrationVertices[2][0] = 0.5;
    __integrationVertices[2][1] = 0.5;
    __integrationVertices[2][2] = 0;
  }
};

class QuadratureFormulaP2Triangle3D
  : public ThreadStaticBase<QuadratureFormulaP2Triangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 6
  };

private:
  TinyVector<numberOfQuadraturePoints,
             TinyVector<3, real_t> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,
                   TinyVector<3, real_t> >&
  vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }


  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP2Triangle3D()
  {
    this->__setQuadratureVertices();
  }
private:
  void __setQuadratureVertices()
  {
    // Picked-up from getfem++
    const real_t a = 0.445948490915965;
    const real_t b = 0.091576213509771;
    const real_t c = 0.111690794839005;
    const real_t d = 0.054975871827661;
    __integrationVertices[0] = TinyVector<3, real_t>(    a,    a,0);
    __weight[0] = c;
    __integrationVertices[1] = TinyVector<3, real_t>(1-2*a,    a,0);
    __weight[1] = c;
    __integrationVertices[2] = TinyVector<3, real_t>(    a,1-2*a,0);
    __weight[2] = c;
    __integrationVertices[3] = TinyVector<3, real_t>(    b,    b,0);
    __weight[3] = d;
    __integrationVertices[4] = TinyVector<3, real_t>(1-2*b,    b,0);
    __weight[4] = d;
    __integrationVertices[5] = TinyVector<3, real_t>(    b,1-2*b,0);
    __weight[5] = d;
  }
};

class QuadratureFormulaQ0Quadrangle3D
  : public ThreadStaticBase<QuadratureFormulaQ0Quadrangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 1
  };

private:
  TinyVector<numberOfQuadraturePoints,
             TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints, TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ0Quadrangle3D()
    : __weight(1)
  {
    __integrationVertices[0] = TinyVector<3,real_t>(0.5,0.5,0);
  }
};

class QuadratureFormulaQ1Quadrangle3D
  : public ThreadStaticBase<QuadratureFormulaQ1Quadrangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 4
  };

private:
  TinyVector<numberOfQuadraturePoints,
             TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints, TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ1Quadrangle3D()
  {
    this->__setQuadratureVertices();
  }

private:
  void __setQuadratureVertices()
  {
    __weight = 1./numberOfQuadraturePoints;

    TinyVector<2, real_t> X;
    
    X[0] = 0.5 - std::sqrt(3.)/6.;
    X[1] = 0.5 + std::sqrt(3.)/6.;

    for (size_t i=0; i<numberOfQuadraturePoints; ++i) {
      __integrationVertices[i][0] = X[i%2];
      __integrationVertices[i][1] = X[(i/2)%2];
      __integrationVertices[i][2] = 0;
    }
  }
};

class QuadratureFormulaQ2Quadrangle3D
  : public ThreadStaticBase<QuadratureFormulaQ2Quadrangle3D>
{
public:
  enum {
    numberOfQuadraturePoints = 9
  };

private:
  TinyVector<numberOfQuadraturePoints, TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints, TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ2Quadrangle3D()
  {
    this->__setQuadratureVertices();
  }

private:
  void __setQuadratureVertices()
  {
    TinyVector<3, real_t> weight1D(5./18.,8./18.,5./18.);

    TinyVector<3, real_t> X;

    X[0] = 0.5 - std::sqrt(15.)/10.;
    X[1] = 0.5;
    X[2] = 0.5 + std::sqrt(15.)/10.;

    for (size_t i=0; i<numberOfQuadraturePoints; ++i) {
      // tensorial weight
      __weight[i] = weight1D[i%3]*weight1D[(i/3)%3];

      __integrationVertices[i][0] = X[i%3];
      __integrationVertices[i][1] = X[(i/3)%3];
      __integrationVertices[i][2] = 0;
    }
  }
};

class QuadratureFormulaQ0Hexahedron
  : public ThreadStaticBase<QuadratureFormulaQ0Hexahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 1
  };

private:
  TinyVector<numberOfQuadraturePoints,TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ0Hexahedron()
    : __weight(1.)
  {
    __integrationVertices[0] = TinyVector<3,real_t>(0.5,0.5,0.5);
  }
};

class QuadratureFormulaQ1Hexahedron
  : public ThreadStaticBase<QuadratureFormulaQ1Hexahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 8
  };

private:
  TinyVector<numberOfQuadraturePoints,TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ1Hexahedron()
  {
    this->__setQuadratureVertices();
  }

private:
  void __setQuadratureVertices()
  {
    __weight = 1./numberOfQuadraturePoints;

    TinyVector<2, real_t> X;
    
    X[0] = 0.5 - std::sqrt(3.)/6.;
    X[1] = 0.5 + std::sqrt(3.)/6.;

    for (size_t i=0; i<numberOfQuadraturePoints; ++i) {
      __integrationVertices[i][0] = X[i%2];
      __integrationVertices[i][1] = X[(i/2)%2];
      __integrationVertices[i][2] = X[(i/4)%2];
    }
  }
};

class QuadratureFormulaQ2Hexahedron
  : public ThreadStaticBase<QuadratureFormulaQ2Hexahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 27
  };
private:
  TinyVector<numberOfQuadraturePoints,TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaQ2Hexahedron()
  {
    this->__setQuadratureVertices();
  }

private:
  void __setQuadratureVertices()
  {
    TinyVector<3, real_t> weight1D(5./18.,8./18.,5./18.);

    TinyVector<3, real_t> X;
    
    X[0] = 0.5 - std::sqrt(15.)/10.;
    X[1] = 0.5;
    X[2] = 0.5 + std::sqrt(15.)/10.;

    for (size_t i=0; i<numberOfQuadraturePoints; ++i) {
      // tensorial weight
      __weight[i] = weight1D[i%3]*weight1D[(i/3)%3]*weight1D[(i/9)%3];

      __integrationVertices[i][0] = X[i%3];
      __integrationVertices[i][1] = X[(i/3)%3];
      __integrationVertices[i][2] = X[(i/9)%3];
    }
  }
};


class QuadratureFormulaP0Tetrahedron
  : public ThreadStaticBase<QuadratureFormulaP0Tetrahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 1
  };

private:
  TinyVector<numberOfQuadraturePoints, TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP0Tetrahedron()
    :__weight(1./6.)
  {
    __integrationVertices[0] = TinyVector<3, real_t>(0.25, 0.25, 0.25);
  }
};

class QuadratureFormulaP1Tetrahedron
  : public ThreadStaticBase<QuadratureFormulaP1Tetrahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 4
  };

private:
  TinyVector<numberOfQuadraturePoints, TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP1Tetrahedron()
  {
    this->__setQuadratureVertices();
  }
private:
  void __setQuadratureVertices()
  {
    __weight = 1./(6.*numberOfQuadraturePoints);

    TinyVector<4, real_t> mu;
    mu[0] = 0.585410196624969;
    mu[1] = 0.138196601125011;
    mu[2] = 0.138196601125011;
    mu[3] = 0.138196601125011;

    TinyVector<4, TinyVector<3,real_t> > vertex;
    vertex[0] = TinyVector<3, real_t>(0,0,0);
    vertex[1] = TinyVector<3, real_t>(1,0,0);
    vertex[2] = TinyVector<3, real_t>(0,1,0);
    vertex[3] = TinyVector<3, real_t>(0,0,1);

    for (size_t i=0; i<4; ++i) {
      __integrationVertices[i] = 0;
      for (size_t j=0; j<4; ++j) {
        __integrationVertices[i] += mu[(i-j)%4]*vertex[j];
      }
    }
  }
};

class QuadratureFormulaP2Tetrahedron
  : public ThreadStaticBase<QuadratureFormulaP2Tetrahedron>
{
public:
  enum {
    numberOfQuadraturePoints = 15
  };

private:
  TinyVector<numberOfQuadraturePoints, TinyVector<3> > __integrationVertices;

  TinyVector<numberOfQuadraturePoints, real_t> __weight;

public:
  const TinyVector<3>& operator[](const size_t& i) const
  {
    return __integrationVertices[i];
  }

  const TinyVector<numberOfQuadraturePoints,TinyVector<3> >& vertices() const
  {
    return __integrationVertices;
  }

  size_t numberOfVertices() const
  {
    return numberOfQuadraturePoints;
  }

  real_t weight(const size_t& i) const
  {
    return __weight[i];
  }

  QuadratureFormulaP2Tetrahedron()
  {
    this->__setQuadratureVertices();
  }
private:
  void __setQuadratureVertices()
  {
    const real_t a = (7.+std::sqrt(15))/34.;
    const real_t b = (7.-std::sqrt(15))/34.;
    const real_t c = (13.-3*std::sqrt(15))/34.;
    const real_t d = (13.+3*std::sqrt(15))/34.;
    const real_t e = (5.-std::sqrt(15))/20.;
    const real_t f = (5.+std::sqrt(15))/20.;
    const real_t h = (2665-14*std::sqrt(15))/226800.;
    const real_t i = (2665+14*std::sqrt(15))/226800.;

    __integrationVertices[0] = TinyVector<3, real_t>(1./4., 1./4., 1./4.);
    __weight[ 0] = 8./405.;
    __integrationVertices[ 1] = TinyVector<3, real_t>(a,a,a);
    __weight[ 1] = h;
    __integrationVertices[ 2] = TinyVector<3, real_t>(a,a,c);
    __weight[ 2] = h;
    __integrationVertices[ 3] = TinyVector<3, real_t>(a,c,a);
    __weight[ 3] = h;
    __integrationVertices[ 4] = TinyVector<3, real_t>(c,a,a);
    __weight[ 4] = h;
    __integrationVertices[ 5] = TinyVector<3, real_t>(b,b,b);
    __weight[ 5] = i;
    __integrationVertices[ 6] = TinyVector<3, real_t>(b,b,d);
    __weight[ 6] = i;
    __integrationVertices[ 7] = TinyVector<3, real_t>(b,d,b);
    __weight[ 7] = i;
    __integrationVertices[ 8] = TinyVector<3, real_t>(d,b,b);
    __weight[ 8] = i;
    __integrationVertices[ 9] = TinyVector<3, real_t>(e,e,f);
    __weight[ 9] = 5./567.;
    __integrationVertices[10] = TinyVector<3, real_t>(e,f,e);
    __weight[10] = 5./567.;
    __integrationVertices[11] = TinyVector<3, real_t>(f,e,e);
    __weight[11] = 5./567.;
    __integrationVertices[12] = TinyVector<3, real_t>(e,f,f);
    __weight[12] = 5./567.;
    __integrationVertices[13] = TinyVector<3, real_t>(f,e,f);
    __weight[13] = 5./567.;
    __integrationVertices[14] = TinyVector<3, real_t>(f,f,e);
    __weight[14] = 5./567.;
  }
};

#endif // QUADRATURE_FORMULA_HPP

---