vtkfunction.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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 3 of the License, 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, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef VTK_FUNCTION_HH
#define VTK_FUNCTION_HH
 
#include <string>
#include <memory>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/io/file/vtk/function.hh>
 
#include <dumux/common/typetraits/typetraits.hh>
 
namespace Dumux {
namespace Vtk {
 
template <typename GridView, typename Mapper, typename F>
struct VectorP0VTKFunction : Dune::VTKFunction<GridView>
{
    enum { dim = GridView::dimension };
    using ctype = typename GridView::ctype;
    using Element = typename GridView::template Codim<0>::Entity;
 
public:
 
    virtual int ncomps() const { return nComps_; }
 
    virtual std::string name() const { return name_; }
 
    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>&) const
    { return accessChooser_(mycomp, mapper_.index(e), IsIndexable<decltype(field_[0])>()); }
 
    VectorP0VTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
    : field_(field), name_(name), nComps_(nComps), mapper_(mapper)
    {
        if (field.size()!=(unsigned int)(mapper.size()))
            DUNE_THROW(Dune::IOError, "VectorP0VTKFunction: size mismatch between field "
                                       << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")");
    }
 
private:
 
    double accessChooser_(int mycomp, int i, std::true_type) const
    { return vectorFieldAccess_(mycomp, i, IsIndexable<decltype(field_[0][0])>()); }
 
    double accessChooser_(int mycomp, int i, std::false_type) const { return field_[i]; }
 
    double vectorFieldAccess_(int mycomp, int i, std::false_type) const { return field_[i][mycomp]; }
 
    double vectorFieldAccess_(int mycomp, int i, std::true_type) const { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
 
    const F& field_;
    const std::string name_;
    int nComps_;
    const Mapper& mapper_;
};
 
template <typename GridView, typename Mapper, typename F>
struct VectorP1VTKFunction : Dune::VTKFunction<GridView>
{
    enum { dim = GridView::dimension };
    using ctype = typename GridView::ctype;
    using Element = typename GridView::template Codim<0>::Entity;
 
public:
 
    virtual int ncomps() const { return nComps_; }
 
    virtual std::string name() const { return name_; }
 
    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const
    {
        const unsigned int dim = Element::mydimension;
        const unsigned int nVertices = e.subEntities(dim);
 
        std::vector<Dune::FieldVector<ctype, 1>> cornerValues(nVertices);
        for (unsigned i = 0; i < nVertices; ++i)
            cornerValues[i] = accessChooser_(mycomp, mapper_.subIndex(e, i, dim), IsIndexable<decltype(field_[0])>());
 
        // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
        const Dune::MultiLinearGeometry<ctype, dim, 1> interpolation(e.type(), std::move(cornerValues));
        return interpolation.global(xi);
    }
 
    VectorP1VTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
    : field_(field), name_(name), nComps_(nComps), mapper_(mapper)
    {
        if (field.size()!=(unsigned int)( mapper.size() ))
            DUNE_THROW(Dune::IOError, "VectorP1VTKFunction: size mismatch between field "
                                      << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")");
    }
private:
 
    double accessChooser_(int mycomp, int i, std::true_type) const
    { return vectorFieldAccess_(mycomp, i, IsIndexable<decltype(field_[0][0])>()); }
 
    double accessChooser_(int mycomp, int i, std::false_type) const { return field_[i]; }
 
    double vectorFieldAccess_(int mycomp, int i, std::false_type) const { return field_[i][mycomp]; }
 
    double vectorFieldAccess_(int mycomp, int i, std::true_type) const { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
 
    const F& field_;
    const std::string name_;
    int nComps_;
    const Mapper& mapper_;
};
 
template <typename GridView, typename Mapper, typename F>
struct VectorP1NonConformingVTKFunction : Dune::VTKFunction<GridView>
{
    enum { dim = GridView::dimension };
    using ctype = typename GridView::ctype;
    using Element = typename GridView::template Codim<0>::Entity;
 
public:
 
    virtual int ncomps() const { return nComps_; }
 
    virtual std::string name() const { return name_; }
 
    virtual double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& xi) const
    {
        const unsigned int dim = Element::mydimension;
        const unsigned int nVertices = e.subEntities(dim);
 
        std::vector<Dune::FieldVector<ctype, 1>> cornerValues(nVertices);
        for (unsigned i = 0; i < nVertices; ++i)
            cornerValues[i] = accessChooser_(mycomp, mapper_.index(e), i, IsIndexable<decltype(field_[0])>());
 
        // (Ab)use the MultiLinearGeometry class to do multi-linear interpolation between scalars
        const Dune::MultiLinearGeometry<ctype, dim, 1> interpolation(e.type(), std::move(cornerValues));
        return interpolation.global(xi);
    }
 
    VectorP1NonConformingVTKFunction(const GridView& gridView, const Mapper& mapper, const F& field, const std::string& name, int nComps)
    : field_(field), name_(name), nComps_(nComps), mapper_(mapper)
    {
        if (field.size()!=(unsigned int)(mapper.size()))
            DUNE_THROW(Dune::IOError, "VectorP1NonConformingVTKFunction: size mismatch between field "
                                      << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")");
    }
private:
 
    double accessChooser_(int mycomp, int eIdx, int cornerIdx, std::true_type) const
    { return fieldAccess_(mycomp, eIdx, cornerIdx, IsIndexable<decltype(field_[0][0])>()); }
 
    double accessChooser_(int mycomp, int eIdx, int cornerIdx, std::false_type) const
    { DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); }
 
    double fieldAccess_(int mycomp, int eIdx, int cornerIdx, std::false_type) const
    { return field_[eIdx][cornerIdx]; }
 
    double fieldAccess_(int mycomp, int eIdx, int cornerIdx, std::true_type) const
    { return field_[eIdx][cornerIdx][mycomp]; }
 
    const F& field_;
    const std::string name_;
    int nComps_;
    const Mapper& mapper_;
};
 
template<class GridView>
class Field
{
    enum { dim = GridView::dimension };
    using ctype = typename GridView::ctype;
    using Element = typename GridView::template Codim<0>::Entity;
 
public:
    // template constructor selects the right VTKFunction implementation
    template <typename F, class Mapper>
    Field(const GridView& gridView, const Mapper& mapper, F const& f,
          const std::string& name, int numComp = 1, int codim = 0,
          Dune::VTK::DataMode dm = Dune::VTK::conforming)
    : codim_(codim)
    {
        if (codim == GridView::dimension)
        {
            if (dm == Dune::VTK::conforming)
                field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
            else
                field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
        }
        else if (codim == 0)
            field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp);
        else
            DUNE_THROW(Dune::NotImplemented, "Only element or vertex quantities allowed.");
    }
 
    virtual ~Field() {}
 
    virtual std::string name () const { return field_->name(); }
 
    virtual int ncomps() const { return field_->ncomps(); }
 
    int codim() const { return codim_; }
 
    virtual double evaluate(int mycomp,
                            const Element &element,
                            const Dune::FieldVector< ctype, dim > &xi) const
    { return field_->evaluate(mycomp, element, xi); }
 
    std::shared_ptr<const Dune::VTKFunction<GridView>> get() const
    { return field_; }
 
private:
    int codim_;
    // can point to anything fulfilling the VTKFunction interface
    std::shared_ptr<Dune::VTKFunction<GridView>> field_;
};
 
} // end namespace Vtk
 
} // end namespace Dumux
 
#endif