fe_interpolator.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
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/space/lagrange1/element.hpp>
#include <kernel/space/lagrange2/element.hpp>
#include <kernel/util/math.hpp>
 
namespace FEAT
{
  namespace Assembly
  {
 
    template<typename ToSpace_, typename FromSpace_>
    struct FEInterpolator
    {
    };
 
    template<typename Trafo_, int shape_dim_ = Trafo_::MeshType::shape_dim>
    struct Lagrange1To2DofAtEntity
    {
      typedef Space::Lagrange2::Element<Trafo_> ToSpace;
      typedef Space::Lagrange1::Element<Trafo_> FromSpace;
 
      template<typename Vector_>
      static void recurse(Vector_& to_coeffs, const Vector_& from_coeffs, const ToSpace& to_space)
      {
        typedef typename Vector_::DataType DataType;
 
        const auto& mesh = to_space.get_trafo().get_mesh();
 
        // We use this to identify the vertex DoF at our shape
        const auto& from_vertex_at_shape_idx = mesh.template get_index_set<shape_dim_, 0>();
 
        // This is for mapping to the correct Lagrange2 DoF
        typename ToSpace::template DofAssignment<shape_dim_, DataType>::Type to_dof_assignment(to_space);
 
        for(Index k(0); k < mesh.get_num_entities(shape_dim_); ++k)
        {
          to_dof_assignment.prepare(k);
 
          // There is only one Lagrange2 DoF at each entity, but this is easier to understand for future reference
          for(int j(0); j < to_dof_assignment.get_num_assigned_dofs(); ++j)
          {
 
            typename Vector_::ValueType tmp(0);
 
            for(int i(0); i < from_vertex_at_shape_idx.get_num_indices(); ++i)
            {
              tmp += from_coeffs(from_vertex_at_shape_idx(k,i));
            }
 
            // Linear interpolation
            tmp *= ( DataType(1)/DataType(from_vertex_at_shape_idx.get_num_indices()) );
 
            // Emplacement operator
            to_coeffs(to_dof_assignment.get_index(j), tmp);
          }
        }
 
        // recurse down
        Lagrange1To2DofAtEntity<Trafo_, shape_dim_-1>::recurse( to_coeffs, from_coeffs, to_space);
      }
    };
 
    template<typename Trafo_>
    struct Lagrange1To2DofAtEntity<Trafo_, 0>
    {
      typedef Space::Lagrange2::Element<Trafo_> ToSpace;
      typedef Space::Lagrange1::Element<Trafo_> FromSpace;
 
      template<typename Vector_>
      static void recurse(Vector_& to_coeffs, const Vector_& from_coeffs, const ToSpace& to_space)
      {
 
        const auto& mesh = to_space.get_trafo().get_mesh();
 
        for(Index i(0); i < mesh.get_num_entities(0); ++i)
        {
          to_coeffs(i, from_coeffs(i));
        }
 
      }
    };
 
    template<typename Trafo_>
    struct FEInterpolator<Space::Lagrange2::Element<Trafo_>, Space::Lagrange1::Element<Trafo_>>
    {
      typedef Trafo_ TrafoType;
      typedef Space::Lagrange2::Element<Trafo_> ToSpace;
      typedef Space::Lagrange1::Element<Trafo_> FromSpace;
 
      template<typename DT_, typename IT_>
      static void interpolate(LAFEM::DenseVector<DT_, IT_>& to_coeffs,
      const LAFEM::DenseVector<DT_, IT_>& from_coeffs,
      const ToSpace& to_space,
      const FromSpace& DOXY(from_space))
      {
        Lagrange1To2DofAtEntity<TrafoType>::recurse(to_coeffs, from_coeffs, to_space);
      }
 
      template<typename DT_, typename IT_, int blocksize_>
      static void interpolate(LAFEM::DenseVectorBlocked<DT_, IT_, blocksize_>& to_coeffs,
      const LAFEM::DenseVectorBlocked<DT_, IT_, blocksize_>& from_coeffs,
      const ToSpace& to_space,
      const FromSpace& DOXY(from_space))
      {
        Lagrange1To2DofAtEntity<TrafoType>::recurse(to_coeffs, from_coeffs, to_space);
      }
    }; //FEInterpolator Lagrange1 to Lagrange2
 
    template<typename Trafo_>
    struct FEInterpolator<Space::Lagrange1::Element<Trafo_>, Space::Lagrange2::Element<Trafo_>>
    {
      typedef Trafo_ TrafoType;
      typedef Space::Lagrange1::Element<Trafo_> ToSpace;
      typedef Space::Lagrange2::Element<Trafo_> FromSpace;
 
      template<typename DT_, typename IT_>
      static void interpolate(LAFEM::DenseVector<DT_, IT_>& to_coeffs,
      const LAFEM::DenseVector<DT_, IT_>& from_coeffs,
      const ToSpace& DOXY(to_space),
      const FromSpace& DOXY(from_space))
      {
        //sanity check for the vector sizes
        XASSERTM(to_coeffs.size() < from_coeffs.size(), "Coefficient vectors do not match!\n To_coeffs size is ");
        //due to the two level ordering of Lagrange1 and Lagrange2 elements it is always sufficient to just truncate our vector
        for(Index i(0); i < to_coeffs.size(); ++i)
        {
          to_coeffs(i, from_coeffs(i));
        }
        //And we are done
      }
 
      template<typename DT_, typename IT_, int blocksize_>
      static void interpolate(LAFEM::DenseVectorBlocked<DT_, IT_, blocksize_>& to_coeffs,
      const LAFEM::DenseVectorBlocked<DT_, IT_, blocksize_>& from_coeffs,
      const ToSpace& DOXY(to_space),
      const FromSpace& DOXY(from_space))
      {
        //sanity check for the vector sizes
        XASSERTM(to_coeffs.size() < from_coeffs.size(), "Coefficient vectors do not match!\n To_coeffs size is ");
        //due to the two level ordering of Lagrange1 and Lagrange2 elements it is always sufficient to only truncate our vector
        for(Index i(0); i < to_coeffs.size(); ++i)
        {
          to_coeffs(i, from_coeffs(i));
        }
        //And we are done
      }
 
    }; //FEInterpolator Lagrange2 to Lagrange1
 
    template<typename Space_>
    struct FEInterpolator<Space_, Space_>
    {
      typedef Space_ ToSpace;
 
      typedef Space_ FromSpace;
 
      template<typename Vector_>
      static void interpolate(Vector_& to_coeffs,
      const Vector_& from_coeffs,
      const ToSpace& DOXY(to_space),
      const FromSpace& DOXY(from_space))
      {
        to_coeffs.copy(from_coeffs, true);
      }
 
    }; //FEInterpolator same space
 
  } // namespace Assembly
} // namespace FEAT