surface_mesh.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>
#include <kernel/geometry/conformal_mesh.hpp>
#include <kernel/geometry/index_calculator.hpp>
 
#include <kernel/util/math.hpp>
 
#include <deque>
#include <memory>
 
namespace FEAT
{
  namespace Geometry
  {
    namespace Atlas
    {
 
      struct SurfaceMeshTraits
      {
        static constexpr bool is_explicit = false;
        static constexpr bool is_implicit = true;
        static constexpr int world_dim = 3;
        static constexpr int param_dim = 2;
      }; // struct SurfaceMeshTraits
 
      template<typename Mesh_>
      class SurfaceMesh :
        public ChartCRTP<SurfaceMesh<Mesh_>, Mesh_, SurfaceMeshTraits>
      {
      public:
        typedef ChartCRTP<SurfaceMesh<Mesh_>, Mesh_, SurfaceMeshTraits> BaseClass;
        typedef typename BaseClass::CoordType CoordType;
        typedef typename BaseClass::WorldPoint WorldPoint;
        typedef Shape::Simplex<2> ShapeType;
        typedef ConformalMesh<ShapeType, 3, CoordType> SurfaceMeshType;
        std::unique_ptr<SurfaceMeshType> _surface_mesh;
 
      private:
        // For computing point projections, we need to know if an edge was already traversed by adding the neighbor
        // across said edge to the search stack. If the surface is concave, it might happen that triangle A says
        // "go to my neighbor, B", and then B says "go to my neighbor, A". It's then clear that the wanted point
        // lies on the edge between A and B.
        bool* _edge_marker;
        bool* _cell_marker;
 
      public:
        SurfaceMesh() = delete;
 
        explicit SurfaceMesh(std::unique_ptr<SurfaceMeshType> surf_mesh) :
          _surface_mesh(std::move(surf_mesh)),
          _edge_marker(new bool[_surface_mesh->get_num_entities(SurfaceMeshType::shape_dim-1)]),
          _cell_marker(new bool[_surface_mesh->get_num_entities(SurfaceMeshType::shape_dim)])
        {
        }
 
        SurfaceMesh(SurfaceMeshType&&) = delete;
 
        virtual ~SurfaceMesh()
        {
          delete[] _edge_marker;
          delete[] _cell_marker;
        }
 
        virtual std::size_t bytes() const override
        {
          if(_surface_mesh != nullptr)
            return _surface_mesh->bytes();
          else
            return std::size_t(0);
        }
 
        virtual String get_type() const override
        {
          return "SurfaceMesh";
        }
 
        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 all world points
          WorldPoint tmp;
          auto& vtx = _surface_mesh->get_vertex_set();
          for(Index i(0); i < vtx.get_num_vertices(); ++i)
          {
            tmp = vtx[i] - origin;
            vtx[i].set_mat_vec_mult(rot, tmp) += offset;
          }
        }
 
        virtual void write(std::ostream& os, const String& sindent) const override
        {
          String indent(sindent);
 
          os << indent << "<SurfaceMesh";
          os << " verts=\" " << _surface_mesh->get_num_entities(0) << "\"";
          os << " trias=\" " << _surface_mesh->get_num_entities(ShapeType::dimension) << "\"";
          os << ">\n";
 
          // increase indent
          indent.resize(indent.size()+2, ' ');
 
          const auto& vtx = _surface_mesh->get_vertex_set();
 
          os << indent << "<Vertices>\n";
          indent.resize(indent.size()+2, ' ');
 
          for(Index i(0); i < vtx.get_num_vertices(); ++i)
          {
            const auto& v = vtx[i];
            os << indent << v[0];
            for(int j(1); j < SurfaceMeshType::world_dim; ++j)
              os << ' ' << v[j];
            os << '\n';
          }
 
          indent.resize(indent.size()-2);
          os << indent << "</Vertices>\n";
 
          // Get vertex at cell index set
          const auto& idx = _surface_mesh->template get_index_set<ShapeType::dimension, 0>();
 
          os << indent << "<Triangles>\n";
          indent.resize(indent.size()+2, ' ');
 
          for(Index i(0); i < idx.get_num_entities(); ++i)
          {
            const auto& idx_loc = idx[i];
            os << indent << idx_loc[0];
            for(int j(1); j < idx.num_indices; ++j)
              os << ' ' << idx_loc[j];
            os << '\n';
          }
 
          indent.resize(indent.size()-2);
          os << indent << "</Triangles>\n";
 
          indent.resize(indent.size()-2);
          os << indent << "</SurfaceMesh>";
        }
 
        void project_point(WorldPoint& point) const
        {
          CoordType signed_distance(0);
          WorldPoint grad_distance(CoordType(0));
 
          project_point(point, signed_distance, grad_distance);
        }
 
        void project_point(WorldPoint& point, CoordType& signed_distance, WorldPoint& grad_distance) const
        {
          // There is nothing to do if the surface mesh does not have any cells
          if(_surface_mesh->get_num_entities(SurfaceMeshType::shape_dim) == Index(0))
            return;
 
          // Vertex at cell information
          const auto& idx(_surface_mesh->template get_index_set<SurfaceMeshType::shape_dim, 0>());
          // The mesh's vertex set so we can get at the coordinates
          const auto& vtx(_surface_mesh->get_vertex_set());
 
          // When a point is found, this will hold the coefficients of its image under the projection wrt. the
          // standard P1/Q1 transformation
          Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1> coeffs(CoordType(0));
          // This will hold the coefficients of the transformation mapping
          Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1> bary(CoordType(0));
 
          // Since we want to project the point an arbitrary distance, we set this to something huge
          CoordType acceptable_distance(Math::huge<CoordType>());
          // Find the facet in the SurfaceMesh we project to and compute the coefficients
          Index best_facet =
            find_cell(coeffs, point, Index(0), acceptable_distance, *_surface_mesh, true);
 
          // Clip the coefficients
          for(int j(0); j < SurfaceMeshType::shape_dim; ++j)
          {
            coeffs[j] = Math::max(coeffs[j], CoordType(0));
            coeffs[j] = Math::min(coeffs[j], CoordType(1));
          }
          // Now that we have the surface mesh cell with the lowest distance, we can compute the projected point
          WorldPoint projected_point(CoordType(0));
 
          // Evaluate the surface mesh trafo for the computed coefficients
          CoordType coeff0(CoordType(1));
          for(int j(0); j < SurfaceMeshType::shape_dim; ++j)
          {
            // Index of the local vertex i in the SurfaceMesh
            Index i(idx(best_facet, j+1));
            projected_point += coeffs[j]*vtx[i];
            coeff0 -= coeffs[j];
          }
 
          projected_point += coeff0*vtx[idx(best_facet, Index(0))];
 
          grad_distance = (projected_point - point);
          signed_distance = grad_distance.norm_euclid();
 
          // If the distance is too small, we set the gradient vector to zero
          if(signed_distance < Math::eps<CoordType>())
            grad_distance.format(CoordType(0));
          else
          {
            grad_distance.normalize();
            WorldPoint nu(get_normal_on_tria(best_facet, coeffs));
            signed_distance *= Math::signum(Tiny::dot(nu, grad_distance));
          }
 
          point = projected_point;
 
        }
 
        WorldPoint get_normal_on_tria(
          const Index facet, const Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1> DOXY(coeffs)) const
        {
          // Vertex at cell information
          const auto& idx(_surface_mesh->template get_index_set<SurfaceMeshType::shape_dim, 0>());
          // The mesh's vertex set so we can get at the coordinates
          const auto& vtx(_surface_mesh->get_vertex_set());
 
          Tiny::Matrix<CoordType, 2, SurfaceMeshType::world_dim> coords(CoordType(0));
          Tiny::Matrix<CoordType, SurfaceMeshType::world_dim, 2> coords_transpose(CoordType(0));
 
          Index i0(idx(facet, 0));
          Index i1(idx(facet, 1));
          Index i2(idx(facet, 2));
 
          coords[0] = vtx[i1] - vtx[i0];
          coords[1] = vtx[i2] - vtx[i0];
 
          for(int i(0); i < coords.m; ++i)
          {
            for(int j(0); j < coords.n; ++j)
              coords_transpose(j,i) = coords(i,j);
          }
 
          WorldPoint nu(Tiny::orthogonal(coords_transpose));
          nu.normalize();
 
          return nu;
        }
 
        void project_meshpart(Mesh_& mesh, const MeshPart<Mesh_>& meshpart) const
        {
          // We need the topology of the meshpart because searching for points is cell-wise and we need to know
          // about the vertices in each cell
          XASSERTM(meshpart.get_topology() != nullptr, "The meshpart needs a topology for SurfaceMesh::project()!");
 
          Index num_facets(meshpart.get_num_entities(SurfaceMeshType::shape_dim));
 
          // There is nothing to do if the meshpart does not have any cells (which refer to facets of the original
          // mesh)
          if(num_facets == Index(0))
            return;
 
          // The number of vertices that need to be projected
          const Index num_verts(meshpart.get_num_entities(0));
          // Mapping of vertices from the meshpart to the real mesh
          const auto& ts_verts(meshpart.template get_target_set<0>());
 
          // Vertex at surface mesh cell information
          const auto& idx_sm(_surface_mesh->template get_index_set<SurfaceMeshType::shape_dim, 0>());
          // Vertex set of the surface mesh
          const auto& vtx_sm(_surface_mesh->get_vertex_set());
 
          // Vertex at facet information for the meshpart
          const auto& idx_mp(meshpart.template get_index_set<SurfaceMeshType::shape_dim, 0>());
          // Subfacet at facet information for the meshpart. i.e. if Mesh is a Hypercube<3> mesh, the boundary is a
          // Hypercube<2> mesh, the facets are quads and THEIR subfacets are edges, so this is edge@quad
          // Needed for computing neighbors.
          const auto& facet_idx_mp(meshpart.template get_index_set
            <SurfaceMeshType::shape_dim, SurfaceMeshType::shape_dim-1>());
 
          // Neighbor information for all facets in the MeshPart. Needed for adding facets to the search stack
          // so we can exploit local information
          typename MeshPart<Mesh_>::template IndexSet<SurfaceMeshType::shape_dim, SurfaceMeshType::shape_dim-1>::Type
            neigh_mp(meshpart.get_num_entities(SurfaceMeshType::shape_dim));
 
          Geometry::Intern::FacetNeighbors::compute(neigh_mp, facet_idx_mp);
 
          // The mesh's vertex set, since we modify the coordinates by projection
          auto& vtx(mesh.get_vertex_set());
 
          // When a point is found, this will hold the coefficients of its image under the projection wrt. the
          // standard P1/Q1 transformation
          Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1> coeffs(CoordType(0));
 
          // For every vertex, we need to know if it was already projected
          bool* vert_todo(new bool[num_verts]);
          for(Index i(0); i < num_verts; ++i)
            vert_todo[i] = true;
 
          // For every facet, we need to know if it was already added to the search stack
          bool* facet_todo(new bool[num_facets]);
          for(Index i(0); i < num_facets; ++i)
            facet_todo[i] = true;
 
          // For every facet we keep a guess, a cell number in the surface mesh where we start looking for vertices
          // from our facet.
          Index* guesses(new Index[num_facets]);
          for(Index i(0); i < num_facets; ++i)
            guesses[i] = Index(0);
 
          // The search stack for all facets. We add neighbors to it to keep the search as local as possible
          std::deque<Index> facet_stack;
          facet_stack.push_front(0);
 
          // Treat all facets
          while(!facet_stack.empty())
          {
            // This is the facet from the mesh part that's being worked on
            Index current_facet(facet_stack.front());
            // This will be the facet from the surface mesh were the vertex was found, set it to the best guess
            // initially
            Index facet_sm(guesses[current_facet]);
 
            facet_stack.pop_front();
            facet_todo[current_facet] = false;
 
            // We want to orthogonally project a point, so we need to figure out an acceptable distance for this
            // projection. Otherwise, it might happen that we project to the far side of the domain
            CoordType acceptable_distance(CoordType(1));
 
            // Compute acceptable distance. Because we do not have the transformation available, we cannot compute
            // i.e. the volume of the facet to deduce an acceptable distance from that. So instead, we pick a
            // random vertex of the facet, compute the distances to all other points, and take the geometric mean.
            for(int j(0); j < idx_mp.num_indices-1; ++j)
            {
              Index i0(idx_mp(current_facet, 0));
              Index i_mp(idx_mp(current_facet, j+1));
              acceptable_distance *= ((vtx[i_mp] - vtx[i0]).norm_euclid());
            }
            acceptable_distance = Math::pow(
              acceptable_distance, CoordType(1)/CoordType(idx_mp.num_indices) - CoordType(1));
 
            // Find all vertices in the current facet
            for(int j(0); j < idx_mp.num_indices; ++j)
            {
              // The index of the local vertex j in the Meshpart
              Index i_mp(idx_mp(current_facet, j));
 
              if(vert_todo[i_mp])
              {
                vert_todo[i_mp] = false;
                const auto& x(vtx[ts_verts[i_mp]]);
 
                // Find the facet in the SurfaceMesh we project to and compute the coefficients. We break at the
                // first facet fulfilling the criteria
                facet_sm = find_cell(coeffs, x, guesses[current_facet], acceptable_distance, * _surface_mesh, false);
                // Update the guess for this facet because there might be more points to project
                guesses[current_facet] = facet_sm;
 
                // This is the index of the local vertex j in the real mesh
                Index i_parent(ts_verts[i_mp]);
                // Clear the old vertex coordinates so we can just add onto that
                vtx[i_parent].format(CoordType(0));
                // Evaluate the surface mesh trafo for the computed coefficients
                CoordType coeff0(CoordType(1));
                for(int i(0); i < SurfaceMeshType::shape_dim; ++i)
                {
                  // Index of the local vertex i in the SurfaceMesh
                  Index i_sm(idx_sm(facet_sm, i+1));
                  vtx[i_parent] += coeffs[i]*vtx_sm[i_sm];
                  coeff0 -= coeffs[i];
                }
                vtx[i_parent] += coeff0*vtx_sm[idx_sm(facet_sm, 0)];
 
              } // vert_todo[i_mp]
            } // vertices in current facet
 
            // Now add all neighbors that are still on the todo list to the search stack
            for(int l(0); l < neigh_mp.num_indices; ++l)
            {
              Index k(neigh_mp(current_facet, l));
              if(k != ~Index(0) && facet_todo[k])
              {
                facet_todo[k] = false;
                facet_stack.push_back(k);
                // The last vertex was found in facet_sm, and since k is a neighbor of current, it cannot be too
                // far off
                guesses[k] = facet_sm;
              }
            } // adding neighbors
 
          } // search stack
 
          // Clean up
          delete[] facet_todo;
          delete[] vert_todo;
          delete[] guesses;
        }
 
        CoordType compute_dist(const WorldPoint& point) const
        {
          WorldPoint grad_distance(CoordType(0));
          return compute_dist(point, grad_distance);
        }
 
        CoordType compute_dist(const WorldPoint& point, WorldPoint& grad_distance) const
        {
          CoordType signed_distance(0);
          WorldPoint projected_point(point);
          project_point(projected_point, signed_distance, grad_distance);
 
          return Math::abs(signed_distance);
        }
 
        CoordType compute_signed_dist(const WorldPoint& point) const
        {
          WorldPoint grad_distance(CoordType(0));
          return compute_signed_dist(point, grad_distance);
        }
 
        CoordType compute_signed_dist(const WorldPoint& point, WorldPoint& grad_distance) const
        {
          CoordType signed_distance(0);
          WorldPoint projected_point(point);
          project_point(projected_point, signed_distance, grad_distance);
 
          return signed_distance;
        }
 
      private:
        template<int world_dim, int sc_, int sx_>
        Index find_cell(
          Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1, sc_>& coeffs,
          const Tiny::Vector<CoordType, world_dim, sx_>& x, Index guess, CoordType acceptable_distance,
          const SurfaceMeshType& mesh, bool find_best_approx) const
        {
          // For every subfacet, we need to know if it was already traversed by adding the neighbor across said
          // facet to the search stack. If the surface is concave, it might happen that facet A says "go to my
          // neighbor, B", and then B says "go to my neighbor, A". It's then clear that the wanted point lies on
          // the facet between A and B.
 
          // Clear traversed information
          for(Index facet(0); facet < mesh.get_num_entities(SurfaceMeshType::shape_dim-1); ++facet)
            _edge_marker[facet] = false;
 
          Index num_cells(mesh.get_num_entities(SurfaceMeshType::shape_dim));
          // Clear cell todo information
          for(Index cell(0); cell < num_cells; ++cell)
            _cell_marker[cell] = false;
 
          // Type of the vertex at cell index set
          typedef typename SurfaceMeshType::template IndexSet<SurfaceMeshType::shape_dim, 0>::Type VertAtCellIdxType;
          // Number of vertices per cell
          static constexpr int num_vert_loc = VertAtCellIdxType::num_indices;
          // Vertex at cell information
          const VertAtCellIdxType& idx(mesh.template get_index_set<SurfaceMeshType::shape_dim, 0>());
          // Facet at cell information
          const auto& facet_idx(mesh.template
              get_index_set<SurfaceMeshType::shape_dim, SurfaceMeshType::shape_dim-1>());
          // The mesh's vertex set so we can get at the coordinates
          const auto& vtx(mesh.get_vertex_set());
 
          const auto& neigh(mesh.get_neighbors());
 
          // This will contain all vertices making up a cell in the mesh. Since the mesh is of shape
          // Simplex<shape_dim>, there are shape_dim+1 vertices.
          Tiny::Matrix<CoordType, num_vert_loc, world_dim> coords(CoordType(0));
 
          // This is the search stack
          std::deque<Index> cell_stack;
          // We start with our guess
          cell_stack.push_front(guess);
 
 
          Index best_approx_cell(~Index(0));
 
          CoordType best_approx_dist(Math::huge<CoordType>());
 
          Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1, sc_> current_coeffs(CoordType(0));
 
          // Search all cells
          while(!cell_stack.empty())
          {
            // Get new cell
            Index current(cell_stack.front());
            cell_stack.pop_front();
            _cell_marker[current] = true;
 
            // Collect all vertices that make up the cell in the coords matrix
            for(int j(0); j < SurfaceMeshType::shape_dim+1; ++j)
            {
              Index i(idx(current, j));
              coords[j] = vtx[i];
            }
 
            // Compute the coefficients of the inverse coordinate transformation. The last coefficient is the
            // signed distance of x to the facet made up by coords
            compute_inverse_coeffs(current_coeffs, x, coords);
 
            // From the coefficients, we can compute the barycentric coordinates (as the surface mesh is of Simplex
            // shape) so we get an idea about in which direction we need to search first
            Tiny::Vector<CoordType, SurfaceMeshType::shape_dim+1> bary(CoordType(1));
            CoordType ortho_dist = compute_bary_and_dist(bary, current_coeffs);
 
            // From all this, we can now determine if there is an orthogonal projection of x onto this facet. For
            // this, we need to know the traversed information as well and the neighbor information of the current
            // facet
            bool success = determine_success(ortho_dist, acceptable_distance, bary, facet_idx[current], _edge_marker);
 
            // Ok, maybe we are done
            if(success && !find_best_approx)
            {
              best_approx_cell = current;
              coeffs = current_coeffs;
              break;
            }
            // Check if we need to update the best approx data
            else if(success && find_best_approx)
            {
              if(ortho_dist < best_approx_dist)
              {
                best_approx_dist = ortho_dist;
                best_approx_cell = current;
                coeffs = current_coeffs;
 
              }
            }
 
            // If we are not done, add all neighbors that have not been searched to the search stack
            for(int l(0); l < SurfaceMeshType::shape_dim+1; ++l)
            {
              // The neighbor over the l-th facet
              Index n(neigh(current, l));
              // Check if the neighbor exists and is still on the to do list
              if(n != ~ Index(0) && !_cell_marker[n])
              {
                Index facet(facet_idx(current,l));
                // If the neighbor is still on the to do list and the facet it lies at is in the traversed list,
                // this means the barycentric coordinate for the corresponding direction was negative, meaning
                // that x lies in that direction. So we add that neighbor to the front.
                if(_edge_marker[facet])
                  cell_stack.push_front(n);
                else
                  cell_stack.push_back(n);
 
                // Since we added the neighbor to the stack, prevent it from being added again later.
                _cell_marker[n] = true;
              }
            } // adding neighbors
 
          } // search loop
 
          // If we get here without success being true, the point could not be found
          XASSERTM(best_approx_cell == ~Index(0), "Could not find point "+stringify(x)+" in SurfaceMesh");
 
          // Otherwise, it was found in current and coeffs is still set from that evaluation
          return best_approx_cell;
        } // find_cell()
 
        template<typename DT_, int n_, int sn_, int sb_>
        static DT_ compute_bary_and_dist(Tiny::Vector<DT_, n_, sb_>& bary, const Tiny::Vector<DT_, n_, sn_>& coeffs)
        {
          bary(0) = CoordType(1);
          for(int j(0); j < n_-1; ++j)
          {
            bary(0) -= coeffs(j);
            bary(j+1) = coeffs(j);
          }
 
          return Math::abs(coeffs(n_-1));
        }
 
        template<typename DT_, int n_, int sn_, typename IdxType_>
        static bool determine_success(
          DT_ ortho_dist,
          DT_ acceptable_distance,
          Tiny::Vector<DT_, n_, sn_>& bary,
          const IdxType_& facet_idx,
          bool* traversed)
        {
          static const DT_ tol = Math::sqrt(Math::eps<DT_>());
          bool success(true);
 
          // Determine success
          if(ortho_dist > acceptable_distance)
            success = false;
 
          int num_negative(0);
 
          for(int l(0); l < n_; ++l)
          {
            if(bary(l) < - tol)
            {
              Index facet(facet_idx[l]);
              ++num_negative;
              if(traversed[facet])
              {
                bary(l) = CoordType(0);
              }
              else
              {
                success = false;
                traversed[facet] = true;
              }
            }
          }
 
          if(num_negative > 1)
            success = false;
 
          return success;
        }
 
 
        template<typename DT_, int sb_, int sp_, int smx_, int snx_>
        void compute_inverse_coeffs(
          Tiny::Vector<DT_, 3, sb_>& coeffs,
          const Tiny::Vector<DT_, 3, sp_>& point,
          const Tiny::Matrix<DT_, 3, 3, smx_, snx_>& x) const
        {
          coeffs.format(DT_(0));
 
          // A will contain the transformation matrix
          Tiny::Matrix<DT_, 3, 3> A(DT_(0));
          // Fill the all rows in the first shape_dim = world_dim-1 columns
          for(int i(0); i < 3; ++i)
          {
            for(int j(0); j < 2; ++j)
              A(i,j) = x(j+1,i) - x(0,i);
          }
 
          // The last column is the additional direction for our augmented simplex and it is orthogonal to the rest
          Tiny::Vector<DT_, 3> ortho;
          Tiny::orthogonal_3x2(ortho, A);
          // Normalize this so the last coefficient is the signed distance
          ortho.normalize();
 
          // Set the last column in A
          for(int i(0); i < 3; ++i)
            A(i,2) = ortho(i);
 
          // This will be the inverse of the transformation matrix
          Tiny::Matrix<DT_, 3 , 3> Ainv(DT_(0));
          Ainv.set_inverse(A);
 
          // This is the solution of A u = point - x[0]
          coeffs = Ainv*(point - x[0]);
        }
 
        template<typename DT_, int world_dim, int sc_, int sp_, int smx_, int snx_>
        void compute_inverse_coeffs(
          Tiny::Vector<DT_, 2, sc_>& coeffs,
          const Tiny::Vector<DT_, world_dim, sp_>& point,
          const Tiny::Matrix<DT_, 2, world_dim, smx_, snx_>& x) const
        {
          static_assert( (world_dim == 2 || world_dim == 3),
          "world dim has to be 2 or 3 for complementary barycentric coordinates");
 
          auto tmp = x[1]-x[0];
          DT_ sp(Tiny::dot(point - x[0],tmp));
          DT_ nsqr(Math::sqr(tmp.norm_euclid()));
          coeffs[0] = sp/nsqr;
          tmp = point - (x[0] + coeffs[0]*(x[1]-x[0]));
          coeffs[1] = tmp.norm_euclid();
        }
      }; // class SurfaceMesh
 
      template<typename Mesh_>
      class SurfaceMeshVertsParser :
        public Xml::MarkupParser
      {
      public:
        typedef typename Mesh_::CoordType CoordType;
        typedef VertexSet<3, CoordType> VertexSetType;
 
      protected:
        VertexSetType& _vertex_set;
        Index _read;
 
      public:
        explicit SurfaceMeshVertsParser(VertexSetType& vertex_set) :
          _vertex_set(vertex_set),
          _read(0)
        {
        }
 
        virtual bool attribs(std::map<String,bool>&) const override
        {
          return true;
        }
 
        virtual void create(int iline, const String& sline, const String&, const std::map<String, String>&, bool closed) override
        {
          if(closed)
            throw Xml::GrammarError(iline, sline, "Invalid closed markup");
        }
 
        virtual void close(int iline, const String& sline) override
        {
          // ensure that we have read all vertices
          if(_read < _vertex_set.get_num_vertices())
            throw Xml::GrammarError(iline, sline, "Invalid terminator; expected point");
        }
 
        virtual std::shared_ptr<MarkupParser> markup(int, const String&, const String&) override
        {
          // no children allowed
          return nullptr;
        }
 
        virtual bool content(int iline, const String& sline) override
        {
          // make sure that we do not read more points than expected
          if(_read >= _vertex_set.get_num_vertices())
            throw Xml::ContentError(iline, sline, "Invalid content; expected terminator");
 
          // split line by whitespaces
          std::deque<String> scoords = sline.split_by_whitespaces();
 
          // check size
          if(scoords.size() != std::size_t(3))
            throw Xml::ContentError(iline, sline, "Invalid number of coordinates");
 
          // get our current vertex
          auto& vtx = _vertex_set[_read];
          for(int i(0); i < 3; ++i)
          {
            if(!scoords.at(std::size_t(i)).parse(vtx[i]))
              throw Xml::ContentError(iline, sline, "Failed to parse vertex coordinate");
          }
 
          // okay, another point done
          ++_read;
 
          return true;
        }
      };
 
      template<typename Mesh_>
      class SurfaceMeshTriasParser :
        public Xml::MarkupParser
      {
      public:
        typedef IndexSet<3> IndexSetType;
 
      protected:
        IndexSetType& _index_set;
        Index _read;
 
      public:
        explicit SurfaceMeshTriasParser(IndexSetType& index_set) :
          _index_set(index_set),
          _read(0)
        {
        }
 
        virtual bool attribs(std::map<String,bool>&) const override
        {
          return true;
        }
 
        virtual void create(int iline, const String& sline, const String&, const std::map<String, String>&, bool closed) override
        {
          if(closed)
            throw Xml::GrammarError(iline, sline, "Invalid closed markup");
        }
 
        virtual void close(int iline, const String& sline) override
        {
          // ensure that we have read all triangles
          if(_read < _index_set.get_num_entities())
            throw Xml::GrammarError(iline, sline, "Invalid terminator; expected triangle indices");
        }
 
        virtual std::shared_ptr<MarkupParser> markup(int, const String&, const String&) override
        {
          // no children allowed
          return nullptr;
        }
 
        virtual bool content(int iline, const String& sline) override
        {
          // make sure that we do not read more points than expected
          if(_read >= _index_set.get_num_entities())
            throw Xml::ContentError(iline, sline, "Invalid content; exprected terminator");
 
          // split line by whitespaces
          std::deque<String> sidx = sline.split_by_whitespaces();
 
          // check size
          if(sidx.size() != std::size_t(3))
            throw Xml::ContentError(iline, sline, "Invalid number of indices");
 
          // get our current vertex
          auto& idx = _index_set[_read];
          for(int i(0); i < 3; ++i)
          {
            if(!sidx.at(std::size_t(i)).parse(idx[i]))
              throw Xml::ContentError(iline, sline, "Failed to parse vertex index");
          }
 
          // okay, another point done
          ++_read;
 
          return true;
        }
      };
 
      template<typename Mesh_, bool enable_ = (Mesh_::shape_dim > 2)>
      class SurfaceMeshChartParser :
        public Xml::DummyParser
      {
      public:
        explicit SurfaceMeshChartParser(std::unique_ptr<ChartBase<Mesh_>>&)
        {
          XABORTM("Thou shall not arrive here");
        }
      };
 
      template<typename Mesh_>
      class SurfaceMeshChartParser<Mesh_, true> :
        public Xml::MarkupParser
      {
      protected:
        typedef typename SurfaceMesh<Mesh_>::SurfaceMeshType SurfaceMeshType;
        std::unique_ptr<ChartBase<Mesh_>>& _chart;
        std::unique_ptr<SurfaceMeshType> _surfmesh;
        Index _nverts, _ntrias;
        bool _have_verts, _have_trias;
 
      public:
        explicit SurfaceMeshChartParser(std::unique_ptr<ChartBase<Mesh_>>& chart) :
          _chart(chart),
          _surfmesh(),
          _nverts(0),
          _ntrias(0),
          _have_verts(false), _have_trias(false)
        {
        }
 
        virtual bool attribs(std::map<String,bool>& attrs) const override
        {
          attrs.emplace("verts", true);
          attrs.emplace("trias", true);
          return true;
        }
 
        virtual void create(
          int iline,
          const String& sline,
          const String&,
          const std::map<String, String>& attrs,
          bool closed) override
        {
          // make sure this one isn't closed
          if(closed)
            throw Xml::GrammarError(iline, sline, "Invalid closed SurfaceMesh markup");
 
          // try to parse the size
          if(!attrs.find("verts")->second.parse(_nverts))
            throw Xml::GrammarError(iline, sline, "Failed to parse number of vertices");
          if(!attrs.find("trias")->second.parse(_ntrias))
            throw Xml::GrammarError(iline, sline, "Failed to parse number of triangles");
 
          // create surface mesh
          Index num_entities[3] = {_nverts, 0, _ntrias};
          _surfmesh.reset(new SurfaceMeshType(num_entities));
        }
 
        virtual void close(int iline, const String& sline) override
        {
          if(!_have_verts)
            throw Xml::GrammarError(iline, sline, "Vertices of SurfaceMesh are missing");
          if(!_have_trias)
            throw Xml::GrammarError(iline, sline, "Triangles of SurfaceMesh are missing");
 
          // deduct topology
          _surfmesh->deduct_topology_from_top();
 
          // note: '_surfmesh' is a Geometry::ConformalMesh, not a Chart::SurfaceMesh
          _chart.reset(new SurfaceMesh<Mesh_>(std::move(_surfmesh)));
        }
 
        virtual bool content(int, const String&) override
        {
          return false;
        }
 
        virtual std::shared_ptr<Xml::MarkupParser> markup(int, const String&, const String& name) override
        {
          if(name == "Vertices")
          {
            _have_verts = true;
            return std::make_shared<SurfaceMeshVertsParser<Mesh_>>(_surfmesh->get_vertex_set());
          }
          if(name == "Triangles")
          {
            _have_trias = true;
            return std::make_shared<SurfaceMeshTriasParser<Mesh_>>(_surfmesh->template get_index_set<2,0>());
          }
          return nullptr;
        }
      };
 
    } // namespace Atlas
  } // namespace Geometry
} // namespace FEAT