discrete_evaluator.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/assembly/asm_traits.hpp>
#include <kernel/trafo/inverse_mapping.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/util/dist.hpp>
 
#include <vector>
 
namespace FEAT
{
  namespace Assembly
  {
    template<typename DT_, int dim_>
    class ScalarDiscreteEvalData
    {
    public:
      static constexpr int dim = dim_;
      typedef DT_ DataType;
      typedef DataType ValueType;
 
      std::vector<ValueType> values;
 
      bool empty() const
      {
        return values.empty();
      }
 
      ValueType mean_value() const
      {
        if(values.empty())
          return DT_(0);
 
        // compute and return mean value
        DT_ r = DT_(0);
        for(const auto& v : values)
          r += v;
        return r / DT_(values.size());
      }
 
      ValueType mean_value_dist(const Dist::Comm& comm) const
      {
        // In a distributed domain, we need to compute the mean value
        // over all patches, as the evaluation point may intersect
        // with several patches. Therefore, we also need to sum up
        // the number of values in addition to the values themselves,
        // so that we can compute the mean value over all patches.
        DT_ val[2] =
        {
          DT_(values.size()),
          DT_(0)
        };
 
        // add all values from this patch
        for(const auto& v : values)
          val[1] += v;
 
        // sum up over all patches
        comm.allreduce(val, val, std::size_t(2), Dist::op_sum);
 
        // compute global mean value
        return (val[0] > DT_(0)) ? (val[1] / val[0]) : DT_(0);
      }
    }; // class ScalarDiscreteEvalData
 
    template<typename DT_, int dim_>
    class VectorDiscreteEvalData
    {
    public:
      static constexpr int dim = dim_;
      typedef DT_ DataType;
      typedef Tiny::Vector<DT_, dim_> ValueType;
 
      std::vector<ValueType> values;
 
      bool empty() const
      {
        return values.empty();
      }
 
      ValueType mean_value() const
      {
        ValueType r;
        r.format();
        if(values.empty())
          return r;
 
        for(const auto& v : values)
          r += v;
        r  *= (DT_(1) / DT_(values.size()));
 
        return r;
      }
 
      ValueType mean_value_dist(const Dist::Comm& comm) const
      {
        ValueType r;
        r.format();
 
        DT_ val[1+dim_] =
        {
          DT_(values.size()),
        };
 
        for(int i(0); i < dim_; ++i)
          val[i+1] = DT_(0);
 
        for(const auto& v : values)
        {
          for(int i(0); i < dim_; ++i)
            val[i+1] += v[i];
        }
 
        comm.allreduce(val, val, std::size_t(dim_+1), Dist::op_sum);
 
        for(int i(0); i < dim_; ++i)
          r[i] = val[i+1];
 
        if(val[0] > DT_(0))
          r *= (DT_(1) / val[0]);
 
        return r;
      }
    }; // class VectorDiscreteEvalData
 
    template<typename DT_, int dim_i_, int dim_j_>
    class MatrixDiscreteEvalData
    {
    public:
      static constexpr int dim_i = dim_i_;
      static constexpr int dim_j = dim_j_;
      typedef DT_ DataType;
      typedef Tiny::Matrix<DT_, dim_i_, dim_j_> ValueType;
 
      std::vector<ValueType> values;
 
      bool empty() const
      {
        return values.empty();
      }
 
      ValueType mean_value() const
      {
        ValueType r;
        r.format();
        if(values.empty())
          return r;
 
        for(const auto& v : values)
          r += v;
        r  *= (DT_(1) / DT_(values.size()));
 
        return r;
      }
 
      ValueType mean_value_dist(const Dist::Comm& comm) const
      {
        ValueType r;
        r.format();
 
        DT_ val[1+dim_i_*dim_j_] =
        {
          DT_(values.size()),
        };
 
        for(int i(0); i < dim_i_; ++i)
        {
          for(int j(0); j < dim_j_; ++j)
          {
            val[i*dim_j_ + j + 1] = DT_(0);
          }
        }
 
        for(const auto& v : values)
        {
          for(int i(0); i < dim_i_; ++i)
          {
            for(int j(0); j < dim_j_; ++j)
            {
              val[i*dim_j_ + j + 1] += v[i][j];
            }
          }
        }
 
        comm.allreduce(val, val, std::size_t(dim_i_*dim_j_+1), Dist::op_sum);
 
        for(int i(0); i < dim_i_; ++i)
        {
          for(int j(0); j < dim_j_; ++j)
          {
            r[i][j] = val[i*dim_j_ + j + 1];
          }
        }
 
        if(val[0] > DT_(0))
          r *= (DT_(1) / val[0]);
 
        return r;
      }
    }; // class MatrixDiscreteEvalData
 
    class DiscreteEvaluator
    {
    public:
      template<typename DT_, typename DTP_, typename IT_, int dim_, int s_, typename Space_>
      static ScalarDiscreteEvalData<DT_, dim_> eval_fe_function(
        const Tiny::Vector<DTP_, dim_, s_>& point,
        const LAFEM::DenseVector<DT_, IT_>& vector,
        const Space_& space)
      {
        // create inverse trafo mapping
        Trafo::InverseMapping<typename Space_::TrafoType, DT_> inv_mapping(space.get_trafo());
 
        // unmap point
        auto inv_map_data = inv_mapping.unmap_point(point);
 
        // evaluate FE function
        return eval_fe_function(inv_map_data, vector, space);
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, typename Space_>
      static ScalarDiscreteEvalData<DT_, dim_> eval_fe_function(
        const Trafo::InverseMappingData<DTP_, dim_, dim_>& inv_map_data,
        const LAFEM::DenseVector<DT_, IT_>& vector,
        const Space_& space)
      {
        // vector type
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
        // space type
        typedef Space_ SpaceType;
        // assembly traits
        typedef AsmTraits1<DT_, SpaceType, TrafoTags::none, SpaceTags::value> AsmTraits;
        // data type
        typedef typename AsmTraits::DataType DataType;
 
        // fetch the trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // create local vector data
        typename AsmTraits::template TLocalVector<DT_> local_vector;
 
        // create matrix scatter-axpy
        typename VectorType::GatherAxpy gather_axpy(vector);
 
        // create evaluation data
        ScalarDiscreteEvalData<DT_, dim_> eval_data;
 
        // loop over all cells of the mesh
        for(std::size_t i(0); i < inv_map_data.size(); ++i)
        {
          const Index& cell = inv_map_data.cells.at(i);
          const auto& dom_point = inv_map_data.dom_points.at(i);
 
          // format local vector
          local_vector.format();
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // finish dof-mapping
          dof_mapping.finish();
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // compute trafo data
          trafo_eval(trafo_data, dom_point);
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // compute function value
          DataType value = DataType(0);
          for(int j(0); j < num_loc_dofs; ++j)
            value += local_vector[j] * space_data.phi[j].value;
 
          // finally, add this result to the eval data
          eval_data.values.push_back(value);
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
 
          // continue with next cell
        }
 
        // return evaluation data
        return eval_data;
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, int s_, typename Space_>
      static VectorDiscreteEvalData<DT_, dim_> eval_fe_function(
        const Tiny::Vector<DTP_, dim_, s_>& point,
        const LAFEM::DenseVectorBlocked<DT_, IT_, dim_>& vector,
        const Space_& space)
      {
        // create inverse trafo mapping
        Trafo::InverseMapping<typename Space_::TrafoType, DT_> inv_mapping(space.get_trafo());
 
        // unmap point
        auto inv_map_data = inv_mapping.unmap_point(point);
 
        // evaluate FE function
        return eval_fe_function(inv_map_data, vector, space);
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, typename Space_>
      static VectorDiscreteEvalData<DT_, dim_> eval_fe_function(
        const Trafo::InverseMappingData<DTP_, dim_, dim_>& inv_map_data,
        const LAFEM::DenseVectorBlocked<DT_, IT_, dim_>& vector,
        const Space_& space)
      {
        // vector type
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim_> VectorType;
        // space type
        typedef Space_ SpaceType;
        // assembly traits
        typedef AsmTraits1<DT_, SpaceType, TrafoTags::none, SpaceTags::value> AsmTraits;
        // data type
        typedef typename AsmTraits::DataType DataType;
 
        // fetch the trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // get maximum number of local dofs
        static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
        // create local vector data
        typedef Tiny::Vector<DataType, dim_> ValueType;
        typedef Tiny::Vector<ValueType, max_local_dofs> LocalVectorType;
        LocalVectorType local_vector;
 
        // create matrix scatter-axpy
        typename VectorType::GatherAxpy gather_axpy(vector);
 
        // create evaluation data
        VectorDiscreteEvalData<DT_, dim_> eval_data;
 
        ValueType value;
 
        // loop over all cells of the mesh
        for(std::size_t i(0); i < inv_map_data.size(); ++i)
        {
          const Index& cell = inv_map_data.cells.at(i);
          const auto& dom_point = inv_map_data.dom_points.at(i);
 
          // format local vector
          local_vector.format();
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // finish dof-mapping
          dof_mapping.finish();
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // compute trafo data
          trafo_eval(trafo_data, dom_point);
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // compute function value
          value.format();
          for(int j(0); j < num_loc_dofs; ++j)
            value.axpy(space_data.phi[j].value, local_vector[j]);
 
          // finally, add this result to the eval data
          eval_data.values.push_back(value);
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
 
          // continue with next cell
        }
 
        // return evaluation data
        return eval_data;
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, int s_, typename Space_>
      static VectorDiscreteEvalData<DT_, dim_> eval_fe_gradient(
        const Tiny::Vector<DTP_, dim_, s_>& point,
        const LAFEM::DenseVector<DT_, IT_>& vector,
        const Space_& space)
      {
        // create inverse trafo mapping
        Trafo::InverseMapping<typename Space_::TrafoType, DT_> inv_mapping(space.get_trafo());
 
        // unmap point
        auto inv_map_data = inv_mapping.unmap_point(point);
 
        // evaluate FE function
        return eval_fe_gradient(inv_map_data, vector, space);
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, typename Space_>
      static VectorDiscreteEvalData<DT_, dim_> eval_fe_gradient(
        const Trafo::InverseMappingData<DTP_, dim_, dim_>& inv_map_data,
        const LAFEM::DenseVector<DT_, IT_>& vector,
        const Space_& space)
      {
        // vector type
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
        // space type
        typedef Space_ SpaceType;
        // assembly traits
        typedef AsmTraits1<DT_, SpaceType, TrafoTags::none, SpaceTags::value|SpaceTags::grad> AsmTraits;
        // data type
        typedef typename AsmTraits::DataType DataType;
 
        // fetch the trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // get maximum number of local dofs
        static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
        // create local vector data
        typedef Tiny::Vector<DataType, max_local_dofs> LocalVectorType;
        LocalVectorType local_vector;
 
        // create matrix scatter-axpy
        typename VectorType::GatherAxpy gather_axpy(vector);
 
        // our local gradient
        Tiny::Vector<DataType, dim_> loc_grad;
        loc_grad.format();
 
        // create evaluation data
        VectorDiscreteEvalData<DT_, dim_> eval_data;
 
        // loop over all cells of the mesh
        for(std::size_t i(0); i < inv_map_data.size(); ++i)
        {
          const Index& cell = inv_map_data.cells.at(i);
          const auto& dom_point = inv_map_data.dom_points.at(i);
 
          // format local vector
          local_vector.format();
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // finish dof-mapping
          dof_mapping.finish();
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // compute trafo data
          trafo_eval(trafo_data, dom_point);
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // compute function value
          loc_grad.format();
          for(int j(0); j < num_loc_dofs; ++j)
          {
            // update gradient
            Tiny::axpy(loc_grad, space_data.phi[j].grad, local_vector[j]);
          }
 
          // finally, add this result to the eval data
          eval_data.values.push_back(loc_grad);
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
 
          // continue with next cell
        }
 
        // return evaluation data
        return eval_data;
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, int s_, typename Space_>
      static MatrixDiscreteEvalData<DT_, dim_, dim_> eval_fe_gradient(
        const Tiny::Vector<DTP_, dim_, s_>& point,
        const LAFEM::DenseVectorBlocked<DT_, IT_, dim_>& vector,
        const Space_& space)
      {
        // create inverse trafo mapping
        Trafo::InverseMapping<typename Space_::TrafoType, DT_> inv_mapping(space.get_trafo());
 
        // unmap point
        auto inv_map_data = inv_mapping.unmap_point(point);
 
        // evaluate FE function
        return eval_fe_gradient(inv_map_data, vector, space);
      }
 
      template<typename DT_, typename DTP_, typename IT_, int dim_, typename Space_>
      static MatrixDiscreteEvalData<DT_, dim_, dim_> eval_fe_gradient(
        const Trafo::InverseMappingData<DTP_, dim_, dim_>& inv_map_data,
        const LAFEM::DenseVectorBlocked<DT_, IT_, dim_>& vector,
        const Space_& space)
      {
        // vector type
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim_> VectorType;
        // space type
        typedef Space_ SpaceType;
        // assembly traits
        typedef AsmTraits1<DT_, SpaceType, TrafoTags::none, SpaceTags::value|SpaceTags::grad> AsmTraits;
        // data type
        typedef typename AsmTraits::DataType DataType;
 
        // fetch the trafo
        const typename AsmTraits::TrafoType& trafo = space.get_trafo();
 
        // create a trafo evaluator
        typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
 
        // create a space evaluator and evaluation data
        typename AsmTraits::SpaceEvaluator space_eval(space);
 
        // create a dof-mapping
        typename AsmTraits::DofMapping dof_mapping(space);
 
        // create trafo evaluation data
        typename AsmTraits::TrafoEvalData trafo_data;
 
        // create space evaluation data
        typename AsmTraits::SpaceEvalData space_data;
 
        // get maximum number of local dofs
        static constexpr int max_local_dofs = AsmTraits::max_local_test_dofs;
 
        // create local vector data
        typedef Tiny::Vector<DataType, dim_> ValueType;
        typedef Tiny::Vector<ValueType, max_local_dofs> LocalVectorType;
        LocalVectorType local_vector;
 
        // create matrix scatter-axpy
        typename VectorType::GatherAxpy gather_axpy(vector);
 
        // our local gradient
        Tiny::Matrix<DataType, dim_, dim_> loc_grad;
        loc_grad.format();
 
        // create evaluation data
        MatrixDiscreteEvalData<DT_, dim_, dim_> eval_data;
 
        // loop over all cells of the mesh
        for(std::size_t i(0); i < inv_map_data.size(); ++i)
        {
          const Index& cell = inv_map_data.cells.at(i);
          const auto& dom_point = inv_map_data.dom_points.at(i);
 
          // format local vector
          local_vector.format();
 
          // initialize dof-mapping
          dof_mapping.prepare(cell);
 
          // gather local vector data
          gather_axpy(local_vector, dof_mapping);
 
          // finish dof-mapping
          dof_mapping.finish();
 
          // prepare trafo evaluator
          trafo_eval.prepare(cell);
 
          // prepare space evaluator
          space_eval.prepare(trafo_eval);
 
          // compute trafo data
          trafo_eval(trafo_data, dom_point);
 
          // compute basis function data
          space_eval(space_data, trafo_data);
 
          // fetch number of local dofs
          const int num_loc_dofs = space_eval.get_num_local_dofs();
 
          // compute function value
          loc_grad.format();
          for(int j(0); j < num_loc_dofs; ++j)
          {
            // update gradient
            loc_grad.add_outer_product(local_vector[j], space_data.phi[j].grad);
          }
 
          // finally, add this result to the eval data
          eval_data.values.push_back(loc_grad);
 
          // finish evaluators
          space_eval.finish();
          trafo_eval.finish();
 
          // continue with next cell
        }
 
        // return evaluation data
        return eval_data;
      }
    }; // class DiscreteEvaluator<...>
  } // namespace Assembly
} // namespace FEAT