chart.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/geometry/mesh_part.hpp>
#include <kernel/util/tiny_algebra.hpp>
#include <kernel/util/xml_scanner.hpp>
 
namespace FEAT
{
  namespace Geometry
  {
    namespace Atlas
    {
      template<typename Mesh_>
      class ChartBase
      {
      public:
        typedef Mesh_ MeshType;
        typedef MeshPart<Mesh_> PartType;
        typedef typename MeshType::VertexSetType VertexSetType;
        typedef typename VertexSetType::VertexType WorldPoint;
        typedef typename VertexSetType::CoordType CoordType;
 
      public:
        ChartBase() = default;
        ChartBase(const ChartBase&) = default;
        ChartBase& operator=(const ChartBase&) = default;
        virtual ~ChartBase() = default;
 
        virtual std::size_t bytes() const
        {
          return std::size_t(0);
        }
 
        virtual bool can_explicit() const = 0;
 
        virtual bool can_implicit() const = 0;
 
        virtual void adapt(MeshType& mesh, const PartType& meshpart) const = 0;
 
        virtual void adapt(PartType& mesh, const PartType& meshpart) const = 0;
 
        virtual void transform(const WorldPoint& origin, const WorldPoint& angles, const WorldPoint& offset) = 0;
 
        virtual WorldPoint map(const WorldPoint& param) const = 0;
 
        virtual WorldPoint project(const WorldPoint& point) const = 0;
 
        virtual CoordType dist(const WorldPoint& point) const = 0;
 
        virtual CoordType dist(const WorldPoint& point, WorldPoint& grad_dist) const = 0;
 
        virtual CoordType signed_dist(const WorldPoint& point) const = 0;
 
        virtual CoordType signed_dist(const WorldPoint& point, WorldPoint& grad_dist) const = 0;
 
        virtual String get_type() const = 0;
 
        virtual void write(std::ostream& os, const String& sindent) const = 0;
      }; // class ChartBase<...>
 
      namespace Intern
      {
        template<bool enable_>
        struct ImplicitChartHelper
        {
          template<typename CT_, typename MT_, typename PT_>
          static bool adapt(const CT_&, MT_&, const PT_&)
          {
            return false;
          }
 
          template<typename CT_, typename WP_>
          static bool project(const CT_&, WP_&)
          {
            return false;
          }
        };
 
        template<>
        struct ImplicitChartHelper<true>
        {
          template<typename CT_, typename MT_, typename PT_>
          static bool adapt(const CT_& chart, MT_& mesh, const PT_& part)
          {
            // First of all, check whether the chart can really perform
            // implicit adaption
            if(!chart.can_implicit())
              return false;
 
            // Try to project the whole meshpart
            chart.project_meshpart(mesh, part);
 
            // okay
            return true;
          }
 
          template<typename CT_, typename WP_>
          static bool project(const CT_& chart, WP_& wp)
          {
            chart.project_point(wp);
            return true;
          }
        };
 
        template<bool enable_>
        struct ExplicitChartHelper
        {
          template<typename CT_, typename MT_, typename PT_>
          static bool adapt(const CT_&, MT_&, const PT_&)
          {
            return false;
          }
 
          template<typename CT_, typename WP_, typename PP_>
          static bool map(const CT_&, WP_&, const PP_&)
          {
            return false;
          }
        };
 
        template<>
        struct ExplicitChartHelper<true>
        {
          template<typename CT_, typename MT_, typename PT_>
          static bool adapt(const CT_& chart, MT_& mesh, const PT_& part)
          {
            // a world point in the mesh
            typedef typename CT_::WorldPoint WorldPoint;
            // a parameter point in the part
            typedef typename CT_::ParamPoint ParamPoint;
 
            // vertex set type of our mesh
            typedef typename MT_::VertexSetType VertexSetType;
 
            // attribute type of our mesh part
            typedef typename PT_::AttributeSetType AttributeSetType;
 
            // First of all, check whether the chart can really perform
            // explicit adaption
            if(!chart.can_explicit())
              return false;
 
            // Try to fetch the parametrization attribute.
            const AttributeSetType* attrib = part.find_attribute("param");
            if(attrib == nullptr)
              return false;
 
            // We have the attribute; check whether it matches our chart
            int attrib_dim = attrib->get_dimension();
            XASSERTM(attrib_dim == CT_::param_dim, "Invalid chart attribute dimension");
 
            // Get the vertex set of the mesh
            VertexSetType& vtx = mesh.get_vertex_set();
 
            // Get the vertex target set of the part
            const TargetSet& vidx = part.template get_target_set<0>();
 
            // loop over all vertices in the mesh part
            Index num_vtx = vidx.get_num_entities();
            for(Index i(0); i < num_vtx; ++i)
            {
              // apply the chart's map function
              chart.map_param(
                 reinterpret_cast<      WorldPoint&>(vtx[vidx[i]]),
                *reinterpret_cast<const ParamPoint*>(attrib->raw_at(i))
                );
            }
 
            // okay
            return true;
          }
 
          template<typename CT_, typename WP_, typename PP_>
          static bool map(const CT_& chart, WP_& wp, const PP_& pp)
          {
            chart.map_param(wp, pp);
            return true;
          }
        };
      } // namespace Intern
 
      template<typename Derived_, typename Mesh_, typename Traits_>
      class ChartCRTP :
        public ChartBase<Mesh_>
      {
      public:
        typedef ChartBase<Mesh_> BaseClass;
        typedef Traits_ TraitsType;
        typedef typename BaseClass::MeshType MeshType;
        typedef typename BaseClass::PartType PartType;
        typedef typename PartType::AttributeSetType AttributeSetType;
        typedef typename MeshType::VertexSetType VertexSetType;
        typedef typename VertexSetType::CoordType CoordType;
 
        static constexpr bool is_explicit = TraitsType::is_explicit;
        static constexpr bool is_implicit = TraitsType::is_implicit;
 
        static constexpr int world_dim = TraitsType::world_dim;
        static constexpr int param_dim = TraitsType::param_dim;
 
        //typedef Tiny::Vector<CoordType, world_dim> WorldPoint;
        typedef typename BaseClass::WorldPoint WorldPoint;
        typedef Tiny::Vector<CoordType, param_dim> ParamPoint;
 
      protected:
        Derived_& cast() {return static_cast<Derived_&>(*this);}
 
        const Derived_& cast() const {return static_cast<const Derived_&>(*this);}
 
      public:
        virtual bool can_explicit() const override
        {
          return is_explicit;
        }
 
        virtual bool can_implicit() const override
        {
          return is_implicit;
        }
 
        virtual void adapt(MeshType& mesh, const PartType& part) const override
        {
          // ensure that the mesh world dimension is compatible
          XASSERTM(MeshType::world_dim == world_dim, "Mesh/Chart world dimension mismatch");
 
          // Try to adapt explicitly
          if(Intern::ExplicitChartHelper<is_explicit>::adapt(cast(), mesh, part))
            return;
 
          // Try to adapt implicitly
          if(Intern::ImplicitChartHelper<is_implicit>::adapt(cast(), mesh, part))
            return;
 
          // If we come out here, we have no way of adaption...
          XABORTM("No mesh adaption possible!");
        }
 
        virtual void adapt(PartType& DOXY(parent_meshpart), const PartType& DOXY(meshpart)) const override
        {
          XABORTM("Adaption of MeshPart not possible yet");
        }
 
        virtual WorldPoint map(const WorldPoint& param) const override
        {
          ASSERTM(is_explicit, "cannot map point: Chart is not explicit");
          WorldPoint wp;
          if constexpr(!is_explicit)
            wp = CoordType(0);
          ParamPoint pp;
          pp.template copy_n<param_dim>(param);
          Intern::ExplicitChartHelper<is_explicit>::map(cast(), wp, pp);
          return wp;
        }
 
        virtual WorldPoint project(const WorldPoint& point) const override
        {
          ASSERTM(is_implicit, "cannot project point: Chart is not implicit");
          WorldPoint wp(point);
          if constexpr(!is_implicit)
            wp = CoordType(0);
          Intern::ImplicitChartHelper<is_implicit>::project(cast(), wp);
          return wp;
        }
 
        virtual CoordType dist(const WorldPoint& point) const override
        {
          return (this->cast()).compute_dist(point);
        }
 
        virtual CoordType dist(const WorldPoint& point, WorldPoint& grad_dist) const override
        {
          return (this->cast()).compute_dist(point, grad_dist);
        }
 
        virtual CoordType signed_dist(const WorldPoint& point) const override
        {
          return (this->cast()).compute_signed_dist(point);
        }
 
        virtual CoordType signed_dist(const WorldPoint& point, WorldPoint& grad_dist) const override
        {
          return (this->cast()).compute_signed_dist(point, grad_dist);
        }
 
      }; // class ChartCRTP<...>
 
    } // namespace Atlas
  } // namespace Geometry
} // namespace FEAT