diff --git a/dumux/io/vtk/function.hh b/dumux/io/vtk/function.hh index d62f4e79fa9d17f8a5a35c020a0857087224df0b..47f3744861c3995b4654066eef5f05f33fc82ef6 100644 --- a/dumux/io/vtk/function.hh +++ b/dumux/io/vtk/function.hh @@ -28,12 +28,15 @@ #include <memory> #include <dune/common/typetraits.hh> +#include <dune/common/fvector.hh> +#include <dune/common/reservedvector.hh> #include <dune/grid/common/mcmgmapper.hh> #include <dune/grid/io/file/vtk/common.hh> #include <dune/grid/io/file/vtk/function.hh> #include <dumux/io/vtk/precision.hh> +#include <dumux/discretization/nonconformingfecache.hh> namespace Dumux::Vtk { @@ -253,6 +256,96 @@ private: }; +/*! + * \ingroup InputOutput + * \brief A VTK function that supports both scalar and vector values for each face. + * This expects the data to be organized by a two-dimensional field storing for + * each element the element-local nodal values. Used for non-conforming spaces + * such as Rannacher-Turek and Crouzeix-Raviart. + * + * \tparam GridView The Dune grid view type + * \tparam Mapper The type used for mapping elements to indices in the field + * \tparam F The field type (either vector of scalars or vectors) + */ +template <typename GridView, typename Mapper, typename F> +struct VectorP1FaceNonConformingVTKFunction : Dune::VTKFunction<GridView> +{ + static constexpr int dim = GridView::dimension; + using ctype = typename GridView::ctype; + using Element = typename GridView::template Codim<0>::Entity; + +public: + + //! return number of components + int ncomps() const final { return nComps_; } + + //! get name + std::string name() const final { return name_; } + + //! evaluate + double evaluate(int mycomp, const Element& e, const Dune::FieldVector<ctype, dim>& localPos) const final + { + const unsigned int dim = Element::mydimension; + const unsigned int nFaces = e.subEntities(1); + + // assemble coefficient vector + using Coefficients = Dune::ReservedVector<ctype, 2*dim>; + Coefficients faceValues; + faceValues.resize(nFaces); + for (int i = 0; i < nFaces; ++i) + faceValues[i] = accessEntry_(mycomp, mapper_.index(e), i); + + using ShapeValues = std::vector<Dune::FieldVector<ctype, 1>>; + ShapeValues N; + feCache_.get(e.type()).localBasis().evaluateFunction(localPos, N); + double result = 0.0; + for (int i = 0; i < N.size(); ++i) + result += faceValues[i]*N[i]; + return result; + } + + //! get output precision for the field + Dumux::Vtk::Precision precision() const final + { return precision_; } + + //! Constructor + VectorP1FaceNonConformingVTKFunction( + const GridView& gridView, + const Mapper& mapper, + const F& field, + const std::string& name, + int nComps, + Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32 + ) + : field_(field), name_(name), nComps_(nComps), mapper_(mapper), precision_(precision) + { + if (field.size()!=(unsigned int)(mapper.size())) + DUNE_THROW(Dune::IOError, "VectorP1FaceNonConformingVTKFunction: size mismatch between field " + << name << " (" << field.size() << ") and mapper (" << mapper.size() << ")"); + } +private: + //! access to the field + double accessEntry_([[maybe_unused]] int mycomp, [[maybe_unused]] int eIdx, [[maybe_unused]] int faceIdx) const + { + if constexpr (Dune::IsIndexable<decltype(std::declval<F>()[0])>{}) + { + if constexpr (Dune::IsIndexable<decltype(std::declval<F>()[0][0])>{}) + return field_[eIdx][faceIdx][mycomp]; + else + return field_[eIdx][faceIdx]; + } + else + DUNE_THROW(Dune::InvalidStateException, "Invalid field type"); + } + + const F& field_; + const std::string name_; + int nComps_; + const Mapper& mapper_; + Dumux::Vtk::Precision precision_; + NonconformingFECache<ctype, ctype, dim> feCache_ = {}; +}; + /*! * \ingroup InputOutput * \brief struct that can hold any field that fulfills the VTKFunction interface @@ -262,7 +355,7 @@ private: template<class GridView> class Field { - enum { dim = GridView::dimension }; + static constexpr int dim = GridView::dimension; using ctype = typename GridView::ctype; using Element = typename GridView::template Codim<0>::Entity; @@ -275,17 +368,41 @@ public: Dumux::Vtk::Precision precision = Dumux::Vtk::Precision::float32) : codim_(codim) { - if (codim == GridView::dimension) + if constexpr (dim > 1) { - if (dm == Dune::VTK::conforming) - field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + if (codim == dim) + { + if (dm == Dune::VTK::conforming) + field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + else + field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + } + else if (codim == 1) + { + if (dm == Dune::VTK::conforming) + DUNE_THROW(Dune::NotImplemented, "Only non-conforming output is implemented"); + else + field_ = std::make_shared< VectorP1FaceNonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + } + else if (codim == 0) + field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); else - field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + DUNE_THROW(Dune::NotImplemented, "Only element, face or vertex quantities allowed."); } - else if (codim == 0) - field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); else - DUNE_THROW(Dune::NotImplemented, "Only element or vertex quantities allowed."); + { + if (codim == dim) + { + if (dm == Dune::VTK::conforming) + field_ = std::make_shared< VectorP1VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + else + field_ = std::make_shared< VectorP1NonConformingVTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + } + else if (codim == 0) + field_ = std::make_shared< VectorP0VTKFunction<GridView, Mapper, F> >(gridView, mapper, f, name, numComp, precision); + else + DUNE_THROW(Dune::NotImplemented, "Only element or vertex quantities allowed."); + } } //! return the name of this field