sphere.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/geometry/atlas/chart.hpp>
 
namespace FEAT
{
  namespace Geometry
  {
    namespace Atlas
    {
      struct SphereTraits
      {
        static constexpr bool is_explicit = false;
        static constexpr bool is_implicit = true;
        static constexpr int world_dim = 3;
        static constexpr int param_dim = 1;
      };
 
      template<typename Mesh_>
      class Sphere :
        public ChartCRTP<Sphere<Mesh_>, Mesh_, SphereTraits>
      {
      public:
        typedef ChartCRTP<Sphere<Mesh_>, Mesh_, SphereTraits> BaseClass;
        typedef typename BaseClass::CoordType CoordType;
        typedef typename BaseClass::WorldPoint WorldPoint;
        typedef typename BaseClass::ParamPoint ParamPoint;
 
      protected:
        WorldPoint _midpoint;
        CoordType _radius;
 
      public:
        explicit Sphere(CoordType mid_x, CoordType mid_y, CoordType mid_z, CoordType radius) :
          _radius(radius)
        {
          XASSERTM(radius > CoordType(0), "invalid Sphere radius");
          _midpoint[0] = mid_x;
          _midpoint[1] = mid_y;
          _midpoint[2] = mid_z;
        }
 
        virtual String get_type() const override
        {
          return "sphere";
        }
 
        virtual void transform(const WorldPoint& origin, const WorldPoint& angles, const WorldPoint& offset) override
        {
          // create rotation matrix
          Tiny::Matrix<CoordType, 3, 3> rot;
          rot.set_rotation_3d(angles(0), angles(1), angles(2));
 
          // transform midpoint
          WorldPoint tmp;
          tmp = _midpoint - origin;
          _midpoint.set_mat_vec_mult(rot, tmp) += offset;
        }
 
        void project_point(WorldPoint& point) const
        {
          WorldPoint grad_dist(point - _midpoint);
          CoordType distance(grad_dist.norm_euclid());
 
          if(distance < Math::eps<CoordType>())
          {
            grad_dist(0) = _radius;
            point += grad_dist;
          }
          else
          {
            point = _midpoint + (_radius / distance)*grad_dist;
          }
        }
 
        void project_meshpart(Mesh_& mesh, const MeshPart<Mesh_>& meshpart) const
        {
          auto& vtx = mesh.get_vertex_set();
          const auto& target_vtx = meshpart.template get_target_set<0>();
 
          for(Index i(0); i < meshpart.get_num_entities(0); ++i)
          {
            project_point(reinterpret_cast<WorldPoint&>(vtx[target_vtx[i]]));
          }
        }
 
        CoordType compute_dist(const WorldPoint& point) const
        {
          return Math::abs(compute_signed_dist(point));
        }
 
        CoordType compute_dist(const WorldPoint& point, WorldPoint& grad_dist) const
        {
          WorldPoint projected(point);
          project_point(projected);
 
          CoordType my_dist(_radius - (point - _midpoint).norm_euclid());
 
          grad_dist = (point - projected);
          grad_dist.normalize();
 
          return Math::abs(my_dist);
        }
 
        CoordType compute_signed_dist(const WorldPoint& point) const
        {
          return (point - _midpoint).norm_euclid() - _radius;
        }
 
        CoordType compute_signed_dist(const WorldPoint& point, WorldPoint& grad_dist) const
        {
          WorldPoint projected(point);
          project_point(projected);
 
          CoordType my_dist(_radius - (point - _midpoint).norm_euclid());
 
          grad_dist = (point - projected);
          grad_dist.normalize();
          grad_dist *= Math::signum(my_dist);
 
          return my_dist;
        }
 
        virtual void write(std::ostream& os, const String& sindent) const override
        {
          os << sindent << "<Sphere";
          os << " radius=\"" << this->_radius << "\"";
          os << " midpoint=\"" << this->_midpoint[0] << " " << this->_midpoint[1] << " " << this->_midpoint[2] << "\"";
          os << " />\n";
        }
      };
 
      template<typename Mesh_, typename ChartReturn_ = ChartBase<Mesh_>>
      class SphereChartParser :
        public Xml::MarkupParser
      {
      private:
        typedef Sphere<Mesh_> ChartType;
        typedef typename ChartType::CoordType CoordType;
        std::unique_ptr<ChartReturn_>& _chart;
 
      public:
        explicit SphereChartParser(std::unique_ptr<ChartReturn_>& chart) :
          _chart(chart)
        {
        }
 
        virtual bool attribs(std::map<String,bool>& attrs) const override
        {
          attrs.emplace("radius", true);
          attrs.emplace("midpoint", true);
          return true;
        }
 
        virtual void create(
          int iline,
          const String& sline,
          const String&,
          const std::map<String, String>& attrs,
          bool) override
        {
          CoordType radius = CoordType(0);
          CoordType mid_x = CoordType(0);
          CoordType mid_y = CoordType(0);
          CoordType mid_z = CoordType(0);
 
          // try to parse the radius
          if(!attrs.find("radius")->second.parse(radius))
            throw Xml::GrammarError(iline, sline, "Failed to parse sphere radius");
          if(radius < CoordType(1E-5))
            throw Xml::GrammarError(iline, sline, "Invalid sphere radius");
 
          // try to parse midpoint
          std::deque<String> mids = attrs.find("midpoint")->second.split_by_whitespaces();
          if(mids.size() != std::size_t(3))
            throw Xml::GrammarError(iline, sline, "Invalid sphere midpoint string");
          if(!mids.front().parse(mid_x) || !mids.at(1).parse(mid_y) || !mids.back().parse(mid_z))
            throw Xml::GrammarError(iline, sline, "Failed to parse sphere midpoint");
 
          _chart.reset(new ChartType(mid_x, mid_y, mid_z, radius));
        }
 
        virtual void close(int, const String&) override
        {
        }
 
        virtual bool content(int, const String&) override
        {
          return false;
        }
 
        virtual std::shared_ptr<Xml::MarkupParser> markup(int, const String&, const String&) override
        {
          return nullptr;
        }
      }; // class SphereChartParser<...>
    } // namespace Atlas
  } // namespace Geometry
} // namespace FEAT