errorEstimator.h

Go to the documentation of this file.

#ifndef __ERRORESTIMATOR_H
#define __ERRORESTIMATOR_H

#include <typeinfo>
#include <cxxabi.h>

#include <Eigen/Dense>

template<typename ConfiguratorType>
class ErrorEstimator
{
public:
    typedef typename ConfiguratorType::InitType MeshType;
    typedef typename ConfiguratorType::RealType RealType;
    // Strangely, lpNorm fails if this is set to
    // "typename ConfiguratorType::VectorType". Bug?
    typedef typename Eigen::VectorXd VectorType;
    typedef typename Eigen::Matrix<RealType, 3, 1> VecType;

    ErrorEstimator(MeshType &grid, RealType gamma = 0.4)
    : _mesh(grid), _gamma(gamma)
    {}

    RealType markAndRefineTriangles(){
        RealType err = markTriangles();
        _mesh.refineMarkedTriangles();
        return err;
    }

    RealType markTriangles()
    {
        VectorType errEstimates(_mesh.getNumTriangs());

        computeLocalErrorEstimates(errEstimates);
        RealType refineAbove = errEstimates.lpNorm<Eigen::Infinity>() * _gamma;

        for (int t = 0; t < _mesh.getNumTriangs(); ++t)
            if (errEstimates[t] > refineAbove)
                _mesh.mark(t);

        return errEstimates.norm();
    }

protected:
    RealType f (const typename ConfiguratorType::Point3DType &x) const
    {
        RealType a = 0.8;
        RealType b = 0.0;
        int r = 0;

        for(int i = 0; i < 4; i++) {
            r = r || ( (x[1] > b) && (fabs(x[0]) + fabs(x[1]-b) < a) );
            b += a *= .75;
        }

        return 2.0 * static_cast<RealType>(r);
    }

    void computeLocalErrorEstimates(VectorType &errEstimates) const
    {
        for (int t = 0; t < _mesh.getNumTriangs(); ++t) {
            RealType h = _mesh.getTriang(t).getEdge(_mesh.getLongestEdgeIndex(t)).norm();
            RealType J_t = jumpResidual(t);
            errEstimates[t] = sqrt(h) * J_t;
        }
    }

    RealType jumpResidual(int t) const
    {
        RealType Jt = 0.0;

        RealType val = getMidpointValue(t);
        for (int e = 0; e < 3; ++e) {
            int neighbour = _mesh.getNeighbour(t, e);
            // If no neighbour
            if (neighbour == -1)
                continue;

            RealType nVal = getMidpointValue(neighbour);
            Jt += fabs(val - nVal);
        }

        return Jt;
    }

    RealType getMidpointValue(int t) const {
        typename ConfiguratorType::ElementType &el = _mesh.getTriang(t);

        RealType res = 0.0;
        for (int i = 0; i < 3; ++i) {
            res += f(el.getNode(i));
        }
        res /= 3.0;

        return res;
    }

    MeshType &_mesh;
    RealType _gamma;
};

#endif /* __ERRORESTIMATOR_H */