p1_vector_vtk_function.hpp

// Copyright (C) 2011-2012 by the BEM++ Authors
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

#ifndef bempp_p1_vector_vtk_function_hpp
#define bempp_p1_vector_vtk_function_hpp

#include "../common/common.hpp"

#include <dune/grid/io/file/vtk/function.hh>

namespace Bempp
{

template<typename GV, typename V>
class P1VectorVTKFunction
        : public Dune::VTKFunction<GV>
{
    typedef Dune::VTKFunction<GV> Base;
    typedef Dune::MultipleCodimMultipleGeomTypeMapper<
    GV, Dune::MCMGVertexLayout> Mapper;

    const V& v;
    std::string s;
    int ncomps_;
    Mapper mapper;

public:
    typedef typename Base::Entity Entity;
    typedef typename Base::ctype ctype;
    using Base::dim;

    virtual int ncomps() const
    {
        return ncomps_;
    }

    virtual double evaluate(int comp, const Entity& e,
                            const Dune::FieldVector<ctype, dim>& xi) const
    {
        double min = 1E100;
        int imin = -1;
        Dune::GeometryType gt = e.type();
        for (int i = 0; i < e.template count<dim>(); ++i)
        {
            Dune::FieldVector<ctype,dim> local =
                    Dune::GenericReferenceElements<ctype,dim>::general(gt)
                    .position(i,dim);
            local -= xi;
            if (local.infinity_norm() < min)
            {
                min = local.infinity_norm();
                imin = i;
            }
        }
        return v[mapper.map(e, imin, dim) * ncomps_ + comp];
    }

    virtual std::string name() const
    {
        return s;
    }

    P1VectorVTKFunction(const GV& gv, const V &v_, const std::string &s_,
                        int ncomps = 1)
        : v(v_), s(s_), ncomps_(ncomps), mapper(gv)
    {
        if (v.size() != (unsigned int)(mapper.size() * ncomps_))
            DUNE_THROW(Dune::IOError, "P1VectorVTKFunction: size mismatch");
    }
};

} // namespace Bempp

#endif