refinedquadraturerule.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef REFINED_QUADRATURE_RULE_HH
#define REFINED_QUADRATURE_RULE_HH
 
#include <dune/common/fmatrix.hh>
#include <dune/fufem/quadraturerules/compositequadraturerule.hh>
 
 
template <class ct, int dim>
class RefinedSimplexQuadratureRule:
    public CompositeQuadratureRule<ct,dim>
{
    protected:
        typedef typename Dune::QuadratureRule<ct,dim> Quad;
        typedef typename Dune::QuadraturePoint<ct,dim> QuadPoint;
        typedef typename Dune::FieldVector<ct,dim> Point;
    public:
        RefinedSimplexQuadratureRule(const Quad& quad)
        : CompositeQuadratureRule<ct,dim>(quad,1)
        {}
};
 
template <class ct>
class RefinedSimplexQuadratureRule<ct,3> :
    public Dune::QuadratureRule<ct,3>
{
    protected:
        typedef typename Dune::QuadratureRule<ct,3> Quad;
        typedef typename Dune::QuadraturePoint<ct,3> QuadPoint;
        typedef typename Dune::FieldVector<ct,3> Point;
 
    public:
 
        RefinedSimplexQuadratureRule(const Quad& quad)
        {
            Point p,pp;
            for(size_t i=0; i<quad.size(); ++i)
            {
                double w = quad[i].weight() / 8.0;
 
                // die vier abgeschnittenen Ecken
                p = quad[i].position();
                p *= 0.5;
                this->push_back(QuadPoint(p, w));
 
                p = quad[i].position();
                p *= 0.5;
                p[0] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                p = quad[i].position();
                p *= 0.5;
                p[1] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                p = quad[i].position();
                p *= 0.5;
                p[2] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                // die aus dem übrigen Oktaeder entstehenden 4 Tetraeder
                Dune::FieldMatrix<double,3,3> A(0.0);
                A[0][0] = -0.5;
                A[0][1] = -0.5;
                A[1][0] =  0.5;
                A[2][1] =  0.5;
                A[2][2] =  0.5;
                pp = quad[i].position();
                A.mv(pp,p);
                p[0] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                A = 0.0;
                A[0][0] = -0.5;
                A[1][1] = -0.5;
                A[1][2] = -0.5;
                A[2][2] =  0.5;
//                pp = quad[i].position();
                A.mv(pp,p);
                p[0] += 0.5;
                p[1] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                A = 0.0;
                A[0][1] =  0.5;
                A[1][0] = -0.5;
                A[1][1] = -0.5;
                A[2][0] =  0.5;
                A[2][1] =  0.5;
                A[2][2] =  0.5;
//                pp = quad[i].position();
                A.mv(pp,p);
                p[1] += 0.5;
                this->push_back(QuadPoint(p, w));
 
                A = 0.0;
                A[0][1] =  0.5;
                A[0][2] =  0.5;
                A[1][1] = -0.5;
                A[2][0] =  0.5;
                A[2][1] =  0.5;
//                pp = quad[i].position();
                A.mv(pp,p);
                p[1] += 0.5;
                this->push_back(QuadPoint(p, w));
            }
        }
 
};
 
#endif