mesh_concentration_function.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/base_header.hpp>
#include <kernel/shape.hpp>
 
#include <kernel/analytic/common.hpp>
#include <kernel/geometry/mesh_node.hpp>
#include <kernel/geometry/intern/face_index_mapping.hpp>
#include <kernel/meshopt/rumpf_trafo.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/property_map.hpp>
#include <kernel/util/dist.hpp>
 
#include <deque>
 
namespace FEAT
{
  namespace Meshopt
  {
 
 
    // Forward declarations
    template<typename DT_>
    class ConcentrationFunctionDefault;
 
    template<typename DT_>
    class ConcentrationFunctionPowOfDist;
 
    template<typename DT_, typename ShapeType_>
#ifndef DOXYGEN
    struct AlignmentPenalty;
#else
    // Note: The following block is only visible for Doxygen. The implementation has to be provided by specializations
    // in ShapeType_
    struct AlignmentPenalty
    {
      typedef DT_ DataType;
      typedef Shape::Simplex<shape_dim> ShapeType;
 
      template<typename Mesh_, typename Dist_, typename EdgeFreqs_>
      static DataType compute_constraint(const Mesh_& DOXY(mesh), const Dist_& DOXY(dist),
      const EdgeFreqs_& DOXY(edge_freqs))
      {
      }
 
      template<typename Mesh_, typename Dist_>
      static DataType compute_constraint(DataType* constraint_vec, const Mesh_& mesh, const Dist_& dist)
      {
      }
 
      template<typename Vector_, typename Mesh_, typename Dist_, typename GradDist_, typename EdgeFreqs_>
      static void add_constraint_grad(
        Vector_& DOXY(grad), const DataType DOXY(alignment_fval), const DataType DOXY(fac),
        const Mesh_& DOXY(mesh), const Dist_& DOXY(dist), const GradDist_& DOXY(grad_dist),
        const EdgeFreqs_& DOXY(edge_freqs))
      {
      }
 
    };
#endif
 
    template<typename DT_, int shape_dim>
    struct AlignmentPenalty<DT_, Shape::Simplex<shape_dim>>
    {
      typedef DT_ DataType;
      typedef Shape::Simplex<shape_dim> ShapeType;
 
      template<typename Mesh_, typename Dist_, typename EdgeFreqs_>
      static DataType compute_constraint(const Mesh_& mesh, const Dist_& dist, const EdgeFreqs_& edge_freqs)
      {
        DataType constraint(0);
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          constraint += edge_freqs(edge)*FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(edge_idx(edge,0)) * dist(edge_idx(edge,1)));
        }
 
        return constraint;
      }
 
      template<typename Mesh_, typename Dist_>
      static DataType compute_constraint(DataType* constraint_vec, const Mesh_& mesh, const Dist_& dist)
      {
        XASSERT(constraint_vec != nullptr);
        typedef Geometry::Intern::FaceIndexMapping<ShapeType, 1, 0> FimType;
 
        DataType constraint(0);
        const auto& idx = mesh.template get_index_set<ShapeType::dimension,0>();
        for(Index cell(0); cell < mesh.get_num_entities(ShapeType::dimension); ++cell)
        {
          constraint_vec[cell] = DataType(0);
          for(int edge(0); edge < Shape::FaceTraits<ShapeType,1>::count; ++edge)
          {
            //These are the indices of the vertices on edge
            int i(FimType::map(edge,0));
            int j(FimType::map(edge,1));
            DataType my_constraint(FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
                eval(- DataType(1)*dist(idx(cell, i)) * dist(idx(cell, j))) );
 
            constraint_vec[cell] += my_constraint;
 
            constraint += my_constraint;
          }
 
        }
 
        return constraint;
      }
 
      template<typename Vector_, typename Mesh_, typename Dist_, typename GradDist_, typename EdgeFreqs_>
      static void add_constraint_grad(
        Vector_& grad, const DataType alignment_fval, const DataType fac, const Mesh_& mesh, const Dist_& dist, const GradDist_& grad_dist, const EdgeFreqs_& edge_freqs)
      {
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        typedef Tiny::Vector<DataType, Mesh_::world_dim> WorldPoint;
 
        WorldPoint grad_loc(DataType(0));
 
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          Index i(edge_idx(edge,0));
          Index j(edge_idx(edge,1));
 
          auto dist_prod =  dist(i) * dist(j);
          // Derivative of the heaviside function
          auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i) + edge_freqs(edge)*grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j) + edge_freqs(edge)*grad_loc);
        }
      }
 
    };
 
    template<typename DT_>
    struct AlignmentPenalty<DT_, Shape::Hypercube<2>>
    {
 
      typedef DT_ DataType;
      typedef Shape::Hypercube<2> ShapeType;
 
      template<typename Mesh_, typename Dist_, typename EdgeFreqs_>
      static DataType compute_constraint(const Mesh_& mesh, const Dist_& dist, const EdgeFreqs_& edge_freqs)
      {
        DataType constraint(0);
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          constraint += edge_freqs(edge)*FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(edge_idx(edge,0)) * dist(edge_idx(edge,1)));
        }
 
        const auto& cell_idx = mesh.template get_index_set<ShapeType::dimension,0>();
        for(Index cell(0); cell < mesh.get_num_entities(ShapeType::dimension); ++cell)
        {
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,0)) * dist(cell_idx(cell,3)));
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,1)) * dist(cell_idx(cell,2)));
        }
 
        return constraint;
      }
 
      template<typename Mesh_, typename Dist_>
      static DataType compute_constraint(DataType* constraint_vec, const Mesh_& mesh, const Dist_& dist)
      {
        XASSERT(constraint_vec != nullptr);
        typedef Geometry::Intern::FaceIndexMapping<ShapeType, 1, 0> FimType;
 
        DataType constraint(0);
        const auto& idx = mesh.template get_index_set<ShapeType::dimension,0>();
        for(Index cell(0); cell < mesh.get_num_entities(ShapeType::dimension); ++cell)
        {
          constraint_vec[cell] = DataType(0);
          for(int edge(0); edge < Shape::FaceTraits<ShapeType,1>::count; ++edge)
          {
            //These are the indices of the vertices on edge
            int i(FimType::map(edge,0));
            int j(FimType::map(edge,1));
            DataType my_constraint(FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
                eval(- DataType(1)*dist(idx(cell, i)) * dist(idx(cell, j))) );
 
            constraint_vec[cell] += my_constraint;
 
            constraint += my_constraint;
          }
 
          DataType my_constraint(FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
              eval( - dist(idx(cell,0)) * dist(idx(cell,3))));
          my_constraint = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(idx(cell,1)) * dist(idx(cell,2)));
 
          constraint_vec[cell] += my_constraint;
 
          constraint += my_constraint;
 
        }
 
        return constraint;
      }
 
      template<typename Vector_, typename Mesh_, typename Dist_, typename GradDist_, typename EdgeFreqs_>
      static void add_constraint_grad(
        Vector_& grad, const DataType alignment_fval, const DataType fac, const Mesh_& mesh, const Dist_& dist,
        const GradDist_& grad_dist, const EdgeFreqs_& edge_freqs)
      {
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        typedef Tiny::Vector<DataType, Mesh_::world_dim> WorldPoint;
 
        WorldPoint grad_loc(DataType(0));
 
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          Index i(edge_idx(edge,0));
          Index j(edge_idx(edge,1));
 
          auto dist_prod =  dist(i) * dist(j);
          // Derivative of the heaviside function
          auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i) + edge_freqs(edge)*grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j) + edge_freqs(edge)*grad_loc);
        }
 
        const auto& cell_idx = mesh.template get_index_set<ShapeType::dimension,0>();
        for(Index cell(0); cell < mesh.get_num_entities(ShapeType::dimension); ++cell)
        {
          Index i(cell_idx(cell,Index(0)));
          Index j(cell_idx(cell,Index(3)));
 
          auto dist_prod = dist(i) * dist(j);
          // Derivative of the heaviside function
          auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i)+grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j)+grad_loc);
 
          i = cell_idx(cell,Index(1));
          j = cell_idx(cell,Index(2));
 
          dist_prod = dist(i) * dist(j);
          // Derivative of the heaviside function
          heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i)+grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j)+grad_loc);
        }
      }
    };
 
    template<typename DT_>
    struct AlignmentPenalty<DT_, Shape::Hypercube<3>>
    {
 
      typedef DT_ DataType;
      typedef Shape::Hypercube<3> ShapeType;
 
      template<typename Mesh_, typename Dist_, typename EdgeFreqs_>
      static DataType compute_constraint(const Mesh_& mesh, const Dist_& dist, const EdgeFreqs_& edge_freqs)
      {
        DataType constraint(0);
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          constraint += edge_freqs(edge)*FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(edge_idx(edge,0)) * dist(edge_idx(edge,1)));
        }
 
        const auto& face_idx = mesh.template get_index_set<2,0>();
        for(Index face(0); face < mesh.get_num_entities(2); ++face)
        {
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(face_idx(face,0)) * dist(face_idx(face,3)));
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(face_idx(face,1)) * dist(face_idx(face,2)));
        }
 
        const auto& cell_idx = mesh.template get_index_set<3,0>();
        for(Index cell(0); cell < mesh.get_num_entities(3); ++cell)
        {
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,0)) * dist(cell_idx(cell,7)));
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,1)) * dist(cell_idx(cell,6)));
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,2)) * dist(cell_idx(cell,5)));
          constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(cell_idx(cell,3)) * dist(cell_idx(cell,4)));
        }
 
        return constraint;
      }
 
      template<typename Mesh_, typename Dist_>
      static DataType compute_constraint(DataType* constraint_vec, const Mesh_& mesh, const Dist_& dist)
      {
        XASSERT(constraint_vec != nullptr);
        typedef Geometry::Intern::FaceIndexMapping<ShapeType, 1, 0> VertAtEdge;
        typedef Geometry::Intern::FaceIndexMapping<ShapeType, 2, 0> VertAtFace;
 
        DataType constraint(0);
        const auto& idx = mesh.template get_index_set<ShapeType::dimension,0>();
        for(Index cell(0); cell < mesh.get_num_entities(ShapeType::dimension); ++cell)
        {
          constraint_vec[cell] = DataType(0);
 
          for(int edge(0); edge < Shape::FaceTraits<ShapeType,1>::count; ++edge)
          {
            //These are the indices of the vertices on edge
            int i(VertAtEdge::map(edge,0));
            int j(VertAtEdge::map(edge,1));
            DataType my_constraint(FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
                eval(- DataType(1)*dist(idx(cell, Index(i))) * dist(idx(cell, Index(j)))) );
 
            constraint_vec[cell] += my_constraint;
 
            constraint += my_constraint;
          }
 
          for(int face(0); face < Shape::FaceTraits<ShapeType,2>::count; ++face)
          {
            // These are the indices of the vertices on edge
            int i0(VertAtFace::map(face,0));
            int i1(VertAtFace::map(face,1));
            int i2(VertAtFace::map(face,2));
            int i3(VertAtFace::map(face,3));
 
            DataType my_constraint(0);
            my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
                eval(- DataType(1)*dist(idx(cell, Index(i0))) * dist(idx(cell, Index(i3))));
            my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
                eval(- DataType(1)*dist(idx(cell, Index(i1))) * dist(idx(cell, Index(i2))));
 
            constraint_vec[cell] += my_constraint;
 
            constraint += my_constraint;
          }
 
          DataType my_constraint(0);
          my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(idx(cell,0)) * dist(idx(cell,7)));
          my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(idx(cell,1)) * dist(idx(cell,6)));
          my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(idx(cell,2)) * dist(idx(cell,5)));
          my_constraint += FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::
            eval( - dist(idx(cell,3)) * dist(idx(cell,4)));
 
          constraint_vec[cell] += my_constraint;
 
          constraint += my_constraint;
 
        }
 
        return constraint;
      }
 
      template<typename Vector_, typename Mesh_, typename Dist_, typename GradDist_, typename EdgeFreqs_>
      static void add_constraint_grad(
        Vector_& grad, const DataType alignment_fval, const DataType fac, const Mesh_& mesh, const Dist_& dist,
        const GradDist_& grad_dist, const EdgeFreqs_& edge_freqs)
      {
        const auto& edge_idx = mesh.template get_index_set<1,0>();
        typedef Tiny::Vector<DataType, Mesh_::world_dim> WorldPoint;
 
        WorldPoint grad_loc(DataType(0));
 
        for(Index edge(0); edge < mesh.get_num_entities(1); ++edge)
        {
          Index i(edge_idx(edge,0));
          Index j(edge_idx(edge,1));
 
          auto dist_prod =  dist(i) * dist(j);
          // Derivative of the heaviside function
          auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i) + edge_freqs(edge)*grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j) + edge_freqs(edge)*grad_loc);
        }
 
        const auto& face_idx = mesh.template get_index_set<2,0>();
        for(Index face(0); face < mesh.get_num_entities(2); ++face)
        {
          Index i(face_idx(face,Index(0)));
          Index j(face_idx(face,Index(3)));
 
          auto dist_prod = dist(i) * dist(j);
          // Derivative of the heaviside function
          auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i)+grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j)+grad_loc);
 
          i = face_idx(face,Index(1));
          j = face_idx(face,Index(2));
 
          dist_prod = dist(i) * dist(j);
          // Derivative of the heaviside function
          heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j)) * grad_dist(i);
          grad(i, grad(i)+grad_loc);
 
          grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i)) * grad_dist(j);
          grad(j, grad(j)+grad_loc);
        }
 
        const auto& cell_idx = mesh.template get_index_set<3,0>();
        for(Index cell(0); cell < mesh.get_num_entities(2); ++cell)
        {
          Index i[4]{cell_idx(cell,Index(0)), cell_idx(cell,Index(1)),cell_idx(cell,Index(2)),cell_idx(cell,Index(3))};
          Index j[4]{cell_idx(cell,Index(7)), cell_idx(cell,Index(6)),cell_idx(cell,Index(5)),cell_idx(cell,Index(2))};
 
          for(int k(0); k < 4; ++k)
          {
 
            auto dist_prod = dist(i[k]) * dist(j[k]);
            // Derivative of the heaviside function
            auto heaviside_der = FEAT::Analytic::Common::template HeavisideRegStatic<DataType>::der_x(-dist_prod);
 
            grad_loc =  (-fac * alignment_fval * heaviside_der * dist(j[k])) * grad_dist(i[k]);
            grad(i[k], grad(i[k]) +grad_loc);
 
            grad_loc =  (-fac * alignment_fval * heaviside_der * dist(i[k])) * grad_dist(j[k]);
            grad(j[k], grad(j[k]) + grad_loc);
          }
 
        }
      }
    };
 
    template
    <
      typename Trafo_,
      typename RefCellTrafo_ = RumpfTrafo<Trafo_, typename Trafo_::CoordType>
    >
    class MeshConcentrationFunctionBase
    {
      public:
        typedef Trafo_ TrafoType;
 
        typedef typename TrafoType::MeshType MeshType;
        typedef typename MeshType::ShapeType ShapeType;
        typedef typename MeshType::CoordType CoordType;
 
        typedef Index IndexType;
 
        typedef LAFEM::DenseVector<CoordType, IndexType> ScalarVectorType;
        typedef LAFEM::DenseVectorBlocked<CoordType, IndexType, MeshType::world_dim> VectorType;
        typedef LAFEM::DenseVectorBlocked<CoordType, IndexType,
        MeshType::world_dim*Shape::FaceTraits<ShapeType,0>::count> GradHType;
        typedef Tiny::Vector<CoordType, MeshType::world_dim> WorldPoint;
      public:
        virtual ~MeshConcentrationFunctionBase()
        {
        }
 
        virtual String info() const = 0;
 
        virtual std::shared_ptr<MeshConcentrationFunctionBase> create_empty_clone() const = 0;
 
        virtual void compute_dist() = 0;
 
        virtual void compute_conc() = 0;
 
        virtual void compute_grad_conc() = 0;
 
        virtual void compute_grad_h(const CoordType& DOXY(sum_det)) = 0;
 
        virtual const ScalarVectorType& get_conc() const = 0;
 
        virtual const ScalarVectorType& get_dist() const = 0;
 
        virtual const VectorType& get_grad_dist() const = 0;
 
        virtual GradHType& get_grad_h() = 0;
 
        virtual const GradHType& get_grad_h() const = 0;
 
        virtual void set_mesh_node(const Geometry::RootMeshNode<MeshType>* DOXY(mesh_node_)) = 0;
 
        virtual bool use_derivative() const = 0;
 
        virtual CoordType compute_constraint() const = 0;
 
        virtual CoordType compute_constraint(CoordType* DOXY(constraint_at_vtx)) const = 0;
 
        virtual void add_constraint_grad(VectorType& DOXY(grad), const CoordType DOXY(constraint), const CoordType DOXY(fac)) const = 0;
 
        virtual void compute_grad_sum_det(const VectorType& DOXY(coords)) = 0;
 
        virtual void add_sync_vecs(std::set<VectorType*>& DOXY(sync_vecs)) = 0;
 
        virtual CoordType& get_sum_conc() = 0;
 
    };
 
    template
    <
      typename ElementalFunction_,
      typename Trafo_,
      typename RefCellTrafo_ = RumpfTrafo<Trafo_, typename Trafo_::CoordType>
    >
    class MeshConcentrationFunction : public MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_>
    {
      public:
        typedef MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_> BaseClass;
 
        typedef ElementalFunction_ ElementalFunction;
        typedef Trafo_ TrafoType;
 
        typedef typename TrafoType::MeshType MeshType;
        typedef typename MeshType::CoordType CoordType;
 
        typedef Index IndexType;
 
        typedef typename MeshType::ShapeType ShapeType;
        typedef LAFEM::DenseVector<CoordType, IndexType> ScalarVectorType;
        typedef LAFEM::DenseVectorBlocked<CoordType, IndexType, MeshType::world_dim> VectorType;
        typedef LAFEM::DenseVectorBlocked
        <CoordType, Index, MeshType::world_dim*Shape::FaceTraits<ShapeType,0>::count> GradHType;
        typedef Tiny::Vector<CoordType, MeshType::world_dim> WorldPoint;
        typedef AlignmentPenalty<CoordType, ShapeType> PenaltyFunction;
 
      protected:
        const Geometry::RootMeshNode<MeshType>* _mesh_node;
        ElementalFunction _func;
 
        ScalarVectorType _edge_freqs;
        ScalarVectorType _dist;
        VectorType _grad_dist;
        CoordType _sum_conc;
        ScalarVectorType _conc;
        VectorType _grad_conc;
        VectorType _grad_sum_det;
        // Each entry in the DenseVectorBlocked represents one cell. Each cell's block contains
        // world_dim*(number of local vertices) entries. Has to be serialized like this because there is no
        // DenseVector that saves a Tiny::Matrix
        GradHType _grad_h;
 
        MeshConcentrationFunction(const ElementalFunction& func_) :
          _mesh_node(nullptr),
          _func(func_),
          _edge_freqs(),
          _dist(),
          _grad_dist(),
          _sum_conc(CoordType(0)),
          _conc(),
          _grad_conc(),
          _grad_sum_det(),
          _grad_h()
          {
          }
 
        // *
        // */
        //explicit MeshConcentrationFunction(const MeshConcentrationFunction& other) :
        //  _mesh_node(other._mesh_node),
        //  _func(other._func),
        //  _sum_conc(other._sum_conc)
        //  {
        //    _dist.clone(other._dist, LAFEM::CloneMode::Deep);
        //    _grad_dist.clone(other._grad_dist, LAFEM::CloneMode::Deep);
        //    _grad_conc.clone(other._grad_conc, LAFEM::CloneMode::Deep);
        //    _grad_sum_det.clone(other._grad_sum_det, LAFEM::CloneMode::Deep);
        //    _grad_h.clone(other._grad_h, LAFEM::CloneMode::Deep);
        //  }
 
      public:
        virtual ~MeshConcentrationFunction()
        {
        }
 
        virtual void add_sync_vecs(std::set<VectorType*>& sync_vecs) override
        {
          if(_func.use_derivative)
          {
            sync_vecs.insert(&_grad_sum_det);
            sync_vecs.insert(&_grad_conc);
          }
        }
 
        virtual const ScalarVectorType& get_conc() const override
        {
          return _conc;
        }
 
        virtual const ScalarVectorType& get_dist() const override
        {
          return _dist;
        }
 
        virtual const VectorType& get_grad_dist() const override
        {
          return _grad_dist;
        }
 
        virtual GradHType& get_grad_h() override
        {
          return _grad_h;
        }
 
        virtual const GradHType& get_grad_h() const override
        {
          return _grad_h;
        }
 
        virtual CoordType& get_sum_conc() override
        {
          return _sum_conc;
        }
 
        virtual String info() const override
        {
          return _func.info();
        }
 
        virtual void set_mesh_node(const Geometry::RootMeshNode<MeshType>* mesh_node_) override
        {
          XASSERT(_mesh_node == nullptr);
 
          _mesh_node = mesh_node_;
 
          _sum_conc = CoordType(0);
 
          Index ndofs(_mesh_node->get_mesh()->get_num_entities(0));
          Index nedges(_mesh_node->get_mesh()->get_num_entities(1));
          Index ncells(_mesh_node->get_mesh()->get_num_entities(ShapeType::dimension));
 
          _edge_freqs = ScalarVectorType(nedges, CoordType(1));
          _dist = ScalarVectorType(ndofs, CoordType(0));
          _grad_dist = VectorType(ndofs,CoordType(0));
          _conc = ScalarVectorType(ncells);
          _grad_conc = VectorType(ndofs, CoordType(0));
          _grad_h = GradHType(ncells,CoordType(0));
          _grad_sum_det = VectorType(ndofs, CoordType(0));
 
          for(const auto& it: _mesh_node->get_mesh_part_names())
          {
            if(it.starts_with("_halo"))
            {
              const auto& edge_ts = _mesh_node->find_mesh_part(it)->template get_target_set<1>();
              for(Index edge(0); edge < edge_ts.get_num_entities(); ++edge)
                _edge_freqs(edge_ts[edge], _edge_freqs(edge_ts[edge])+CoordType(1));
            }
          }
          _edge_freqs.component_invert(_edge_freqs, CoordType(1));
        }
 
        virtual bool use_derivative() const override
        {
          return _func.use_derivative;
        }
 
        virtual void compute_conc() override
        {
          // Total number of cells in the mesh
          const Index ncells(_mesh_node->get_mesh()->get_num_entities(ShapeType::dimension));
          // Index set for local/global numbering
          const auto& idx = _mesh_node->get_mesh()->template get_index_set<ShapeType::dimension,0>();
 
          this->_conc.format(CoordType(0));
 
          _sum_conc = CoordType(0);
 
          for(Index cell(0); cell < ncells; ++cell)
          {
            CoordType avg_dist(0);
 
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              avg_dist += this->_dist(idx(cell,j));
            }
 
            avg_dist = avg_dist/CoordType(Shape::FaceTraits<ShapeType,0>::count);
 
            this->_conc(cell, _func.conc_val(avg_dist));
            _sum_conc += this->_conc(cell);
          }
        }
 
        virtual void compute_grad_sum_det(const VectorType& coords) override
        {
          if(_func.use_derivative)
          {
            RefCellTrafo_::compute_grad_sum_det(_grad_sum_det, coords, *(_mesh_node->get_mesh()));
          }
        }
 
        template<typename Tgrad_, typename Tl_, typename Tgradl_>
        void compute_grad_conc_local(Tgrad_& grad_loc_, const Tl_& dist_loc_, const Tgradl_& grad_dist_loc_)
        {
          grad_loc_.format(CoordType(0));
 
          // This will be the average of the levelset values at the vertices
          CoordType val(0);
          for(int i(0); i < Shape::FaceTraits<ShapeType,0>::count; ++i)
          {
            val += dist_loc_(i);
          }
 
          val = val/CoordType(Shape::FaceTraits<ShapeType,0>::count);
 
          for(int i(0); i < Shape::FaceTraits<ShapeType,0>::count; ++i)
          {
            grad_loc_[i] = _func.conc_der(val)/CoordType(Shape::FaceTraits<ShapeType,0>::count) * grad_dist_loc_[i];
          }
        }
 
        virtual void compute_grad_h(const CoordType& sum_det) override
        {
          XASSERT(_mesh_node != nullptr);
 
          _grad_h.format();
 
          if(_func.use_derivative)
          {
 
            // Index set for local/global numbering
            auto& idx = _mesh_node->get_mesh()->template get_index_set<ShapeType::dimension,0>();
            // This will hold the levelset values at the mesh vertices for one element
            FEAT::Tiny::Vector<CoordType, Shape::FaceTraits<ShapeType,0>::count> dist_loc;
            // This will hold the levelset gradient values for one element for passing to other routines
            FEAT::Tiny::Matrix<CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> grad_dist_loc;
            // This will hold the local gradient values for one element for passing to other routines
            FEAT::Tiny::Matrix<CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> grad_loc(0);
            // This will hold the computed local gradient values for one element for copy assigning to the blocked
            // datatype
            FEAT::Tiny::Vector<CoordType, MeshType::world_dim*Shape::FaceTraits<ShapeType,0>::count> tmp(0);
 
            CoordType exponent(CoordType(1)/CoordType(MeshType::world_dim) - CoordType(1));
 
            for(Index cell(0); cell < _mesh_node->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
            {
              grad_loc.format();
              for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
              {
                Index i(idx(cell, j));
                // Get levelset
                dist_loc(j) = _dist(i);
                // Get levelset gradient
                grad_dist_loc[j] = _grad_dist(i);
              }
 
              compute_grad_conc_local(grad_loc, dist_loc, grad_dist_loc);
 
              for(int d(0); d < MeshType::world_dim; ++d)
              {
                for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
                {
                  Index i(idx(cell, j));
 
                  tmp(j*MeshType::world_dim +d) =
                    CoordType(1)/CoordType(MeshType::world_dim)*Math::pow(_conc(cell)/_sum_conc*sum_det,exponent)
                    *( _conc(cell)*(_grad_sum_det(i)(d)*_sum_conc + sum_det*_grad_conc(i)(d) )
                        + grad_loc(j,d) * sum_det *_sum_conc)
                    / Math::sqr(_sum_conc);
                }
              }
              _grad_h(cell, tmp + _grad_h(cell));
 
            }
          } // _func.use_derivative
 
        } // compute_grad_h
 
        virtual void compute_grad_conc() override
        {
          // Clear the old gradient
          _grad_conc.format();
 
          // Index set for local/global numbering
          auto& idx = _mesh_node->get_mesh()->template get_index_set<ShapeType::dimension,0>();
 
          // This will hold the coordinates for one element for passing to other routines
          FEAT::Tiny::Matrix <CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> x;
          // Local cell dimensions for passing to other routines
          FEAT::Tiny::Vector <CoordType,MeshType::world_dim> h;
          // This will hold the local gradient for one element for passing to other routines
          FEAT::Tiny::Matrix <CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> grad_loc;
          // This will hold the levelset values at the mesh vertices for one element
          FEAT::Tiny::Vector <CoordType, Shape::FaceTraits<ShapeType,0>::count> grad_dist_loc;
          // This will hold the levelset gradient values for one element for passing to other routines
          FEAT::Tiny::Matrix <CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> dist_loc;
 
          // Compute the functional value for each cell
          for(Index cell(0); cell < _mesh_node->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            // Collect levelset and levelset grad values
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              // Global vertex/dof index
              Index i(idx(cell, j));
              // Get levelset
              grad_dist_loc(j) = _dist(i);
              // Get levelset gradient
              dist_loc[j] = _grad_dist(i);
            }
 
            // Compute gradient of the concentration on this cell
            compute_grad_conc_local(grad_loc, grad_dist_loc, dist_loc);
 
            // Add local contributions to global gradient vector
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              Index i(idx(cell, j));
              this->_grad_conc(i, _grad_conc(i) + grad_loc[j]);
            }
          }
 
        } // compute_grad_conc()
 
        virtual CoordType compute_constraint() const override
        {
            return PenaltyFunction::compute_constraint(*(_mesh_node->get_mesh()), this->_dist, this->_edge_freqs);
        }
 
        virtual CoordType compute_constraint(CoordType* constraint_at_vtx) const override
        {
          XASSERT(constraint_at_vtx != nullptr);
          return PenaltyFunction::compute_constraint(constraint_at_vtx, *(_mesh_node->get_mesh()), this->_dist);
        }
 
        virtual void add_constraint_grad(
          VectorType& grad, const CoordType constraint, const CoordType fac) const override
        {
          PenaltyFunction::add_constraint_grad(
            grad, constraint, fac, *(_mesh_node->get_mesh()), _dist, _grad_dist, this->_edge_freqs);
        }
 
    }; // class MeshConcentrationFunction
 
    template
    <
      typename ElementalFunction_,
      typename Trafo_,
      typename RefCellTrafo_ = RumpfTrafo<Trafo_, typename Trafo_::CoordType>
    >
    class ChartDistanceFunction : public MeshConcentrationFunction<ElementalFunction_, Trafo_, RefCellTrafo_>
    {
      public:
        typedef ElementalFunction_ ElementalFunction;
        typedef Trafo_ TrafoType;
        typedef typename TrafoType::MeshType MeshType;
        typedef typename MeshType::CoordType CoordType;
        typedef Tiny::Vector<CoordType, MeshType::world_dim> WorldPoint;
 
        typedef MeshConcentrationFunction<ElementalFunction_, Trafo_, RefCellTrafo_> DirectBaseClass;
        typedef MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_> BaseClass;
        typedef Index IndexType;
 
        typedef typename MeshType::ShapeType ShapeType;
        typedef LAFEM::DenseVector<CoordType, IndexType> ScalarVectorType;
        typedef LAFEM::DenseVectorBlocked<CoordType, IndexType, MeshType::world_dim> VectorType;
 
      protected:
        std::deque<String> _chart_list;
        const String _operation;
 
      public:
        ChartDistanceFunction() = delete;
 
        ChartDistanceFunction(const ElementalFunction& func_, const std::deque<String>& chart_list_, const String& operation_):
          DirectBaseClass(func_),
          _chart_list(),
          _operation(operation_)
          {
            XASSERTM(chart_list_.size() > size_t(0), "Empty chart list.");
 
            for(const auto& it:chart_list_)
              _chart_list.push_back(it);
          }
 
        //explicit ChartDistanceFunction(const ChartDistanceFunction& other) :
        //  BaseClass(other.BaseClass),
        //  _chart_list(other._chart_list)
        //  {
        //  }
 
        virtual ~ChartDistanceFunction()
        {
        }
 
        virtual std::shared_ptr<BaseClass> create_empty_clone() const override
        {
          std::shared_ptr<BaseClass> result(nullptr);
          result = std::make_shared<ChartDistanceFunction>(this->_func, _chart_list, _operation);
          return result;
        }
 
        virtual void compute_dist() override
        {
          if(_operation == "add")
            compute_dist_add();
          else if(_operation == "max")
            compute_dist_max();
          else if(_operation == "min")
            compute_dist_min();
          else
            XABORTM("Unknown operation "+_operation);
        }
 
        void compute_dist_add()
        {
          XASSERT(this->_mesh_node != nullptr);
 
          const auto& vtx = this->_mesh_node->get_mesh()->get_vertex_set();
 
          WorldPoint my_dist_vec(CoordType(0));
          WorldPoint tmp(CoordType(0));
 
          for(Index i(0); i < this->_mesh_node->get_mesh()->get_num_entities(0); ++i)
          {
            CoordType my_dist(0);
            my_dist_vec.format(CoordType(0));
 
            for(const auto& it:_chart_list)
            {
              auto* chart = this->_mesh_node->get_atlas()->find_mesh_chart(it);
              if(chart == nullptr)
                XABORTM("Could not find chart "+it);
 
              my_dist += Math::abs(chart->signed_dist(vtx[i], tmp));
              my_dist_vec += tmp;
 
            }
            this->_dist(i, my_dist);
 
            // Because we added distance function gradient vectors, we have to normalize again if possible
            CoordType my_norm(my_dist_vec.norm_euclid());
            if(my_norm > Math::eps<CoordType>())
              my_dist_vec *= (CoordType(1)/my_norm);
 
            this->_grad_dist(i, my_dist_vec);
          }
        }
 
        void compute_dist_max()
        {
          XASSERT(this->_mesh_node != nullptr);
 
          const auto& vtx = this->_mesh_node->get_mesh()->get_vertex_set();
 
          WorldPoint my_dist_vec(CoordType(0));
          WorldPoint tmp(CoordType(0));
 
          for(Index i(0); i < this->_mesh_node->get_mesh()->get_num_entities(0); ++i)
          {
            CoordType my_dist(-Math::huge<CoordType>());
            my_dist_vec.format(CoordType(0));
 
            for(const auto& it:_chart_list)
            {
              auto* chart = this->_mesh_node->get_atlas()->find_mesh_chart(it);
              if(chart == nullptr)
                XABORTM("Could not find chart "+it);
 
              CoordType this_dist = chart->signed_dist(vtx[i], tmp);
              if(this_dist > my_dist)
              {
                my_dist = this_dist;
                my_dist_vec = tmp;
              }
 
            }
 
            this->_dist(i, my_dist);
            this->_grad_dist(i, my_dist_vec);
          }
        }
        void compute_dist_min()
        {
          XASSERT(this->_mesh_node != nullptr);
 
          const auto& vtx = this->_mesh_node->get_mesh()->get_vertex_set();
 
          WorldPoint my_dist_vec(CoordType(0));
          WorldPoint tmp(CoordType(0));
 
          for(Index i(0); i < this->_mesh_node->get_mesh()->get_num_entities(0); ++i)
          {
            CoordType my_dist(Math::huge<CoordType>());
            my_dist_vec.format(CoordType(0));
 
            for(const auto& it:_chart_list)
            {
              auto* chart = this->_mesh_node->get_atlas()->find_mesh_chart(it);
              if(chart == nullptr)
                XABORTM("Could not find chart "+it);
 
              CoordType this_dist = chart->signed_dist(vtx[i], tmp);
              if(Math::abs(this_dist) < Math::abs(my_dist))
              {
                my_dist = this_dist;
                my_dist_vec = tmp;
              }
 
            }
 
            this->_dist(i, my_dist);
            this->_grad_dist(i, my_dist_vec);
          }
        }
 
        static String name()
        {
          return "ChartDistanceFunction<"+ElementalFunction::name()+">";
        }
 
        virtual String info() const override
        {
          const Index width(30);
 
          String msg = name() + ":\n";
 
          msg += DirectBaseClass::info() + "\n";
          msg += String("Operation").pad_back(width, '.') + String(": ") + _operation;
 
          for(const auto& it:_chart_list)
          {
            msg += String("\nDistance chart").pad_back(width, '.') + String(": ") + it;
          }
 
          return msg;
        }
 
    }; // class ChartDistanceFunction
 
    template<typename Trafo_, typename RefCellTrafo_>
    struct MeshConcentrationFunctionFactory
    {
      typedef typename Trafo_::MeshType MeshType;
      typedef typename MeshType::CoordType CoordType;
 
      static std::shared_ptr<MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_>>
        create(const String& section_key, PropertyMap* config)
        {
          XASSERT(config != nullptr);
 
          std::shared_ptr<MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_>> result(nullptr);
 
          // Get the configuration section
          auto my_section = config->query_section(section_key);
          XASSERTM(my_section != nullptr, "Config is missing the referenced "+section_key+" section.");
 
          auto type_p = my_section->query("type");
          XASSERTM(type_p.second, "Section "+section_key+" is missing the mandatory type key.");
 
          if(type_p.first== "ChartDistance")
          {
            auto chart_list_p = my_section->query("chart_list");
            XASSERTM(chart_list_p.second, "Section "+section_key+" is missing the mandatory chart_list key.");
 
            String operation("min");
            auto operation_p = my_section->query("operation");
            if(operation_p.second)
              operation = operation_p.first;
 
            std::deque<String> chart_list = chart_list_p.first.split_by_whitespaces();
 
            auto function_type_p = my_section->query("function_type");
            if(function_type_p.second)
            {
              // At the moment, there are only default and PowOfDist
              if(function_type_p.first == "PowOfDist")
              {
                CoordType minval(0);
                CoordType exponent(1);
                bool use_derivative(true);
 
                auto minval_p = my_section->query("minval");
                if(minval_p.second)
                  minval = std::stod(minval_p.first);
 
                auto exponent_p = my_section->query("exponent");
                if(exponent_p.second)
                  exponent = std::stod(exponent_p.first);
 
                auto use_derivative_p = my_section->query("use_derivative");
                if(use_derivative_p.second)
                  use_derivative = (std::stoi(use_derivative_p.first) == 1);
 
                typedef ConcentrationFunctionPowOfDist<CoordType> ElementalFunction;
                ElementalFunction my_func(minval, exponent, use_derivative);
 
                auto real_result = std::make_shared<ChartDistanceFunction<ElementalFunction, Trafo_, RefCellTrafo_>>
                  (my_func, chart_list, operation);
 
                result = real_result;
              }
              // Default case
              else
              {
                typedef ConcentrationFunctionDefault<CoordType> ElementalFunction;
                ElementalFunction my_func;
 
                auto real_result = std::make_shared<ChartDistanceFunction<ElementalFunction, Trafo_, RefCellTrafo_>>
                  (my_func, chart_list, operation);
 
                result = real_result;
              }
            }
 
          }
 
          XASSERTM(result != nullptr, "Unknown conc_function type!");
 
          return result;
        }
    };
 
    template<typename DT_>
    class ConcentrationFunctionDefault
    {
      public:
        typedef DT_ DataType;
 
        const bool use_derivative;
 
        explicit ConcentrationFunctionDefault():
          use_derivative(false)
        {
        }
 
        ConcentrationFunctionDefault(const ConcentrationFunctionDefault& other) :
          use_derivative(other.use_derivative)
        {
        }
 
        virtual ~ConcentrationFunctionDefault()
        {
        }
 
        static String name()
        {
          return "ConcentrationFunctionDefault";
        }
 
        virtual String info() const
        {
          const Index pad_width(30);
 
          String msg = name() + String(":\n");
 
          msg += String("Function").pad_back(pad_width, '.') + String(": c(d) = |d|\n");
          msg += String("use_derivative").pad_back(pad_width, '.') + String(": ") += stringify(use_derivative);
 
          return msg;
        }
 
        DataType conc_val(DT_ dist) const
        {
          return Math::abs(dist);
        }
 
        DataType conc_der(DT_ DOXY(dist)) const
        {
          return DataType(0);
        }
    };
 
    template<typename DT_>
    class ConcentrationFunctionPowOfDist
    {
      public:
        typedef DT_ DataType;
        const bool use_derivative;
 
      private:
        const DataType _minval;
        const DataType _exponent;
 
      public:
        explicit ConcentrationFunctionPowOfDist(DataType minval_, DataType exponent_, bool use_derivative_ = true):
          use_derivative(use_derivative_),
          _minval(minval_),
          _exponent(exponent_)
          {
          }
 
        //ConcentrationFunctionPowOfDist(const ConcentrationFunctionPowOfDist& other) :
        //  use_derivative(other.use_derivative),
        //  _minval(other._minval),
        //  _exponent(other._exponent)
        //  {
        //  }
 
        ConcentrationFunctionPowOfDist(const ConcentrationFunctionPowOfDist&) = default;
        ConcentrationFunctionPowOfDist& operator=(const ConcentrationFunctionPowOfDist&) = default;
 
        virtual ~ConcentrationFunctionPowOfDist()
        {
        }
 
        static String name()
        {
          return "ConcentrationFunctionPowOfDist";
        }
 
        virtual String info() const
        {
          const Index pad_width(30);
 
          String msg = name() + String(":\n");
 
          msg += String("Function").pad_back(pad_width, '.') + String(": c(d) = (alpha + |d|^beta)\n");
          msg += String("alpha").pad_back(pad_width, '.') + String(": ") + stringify_fp_sci(_minval) + "\n";
          msg += String("beta").pad_back(pad_width, '.') + String(": ") + stringify_fp_sci(_exponent) + "\n";
          msg += String("use_derivative").pad_back(pad_width, '.') + String(": ") + stringify(use_derivative);
 
          return msg;
        }
 
        DataType conc_val(DataType dist) const
        {
          return Math::pow(_minval + Math::abs(dist), _exponent);
        }
 
        DataType conc_der(DataType dist) const
        {
          return _exponent*Math::pow(_minval + Math::abs(dist), _exponent - DataType(1))*Math::signum(dist);
        }
    };
  } // namespace Meshopt
} // namespace FEAST