curvatureoperator.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_INTERFACE_CURVATUREOPERATOR_HH
#define DUNE_MMESH_INTERFACE_CURVATUREOPERATOR_HH
 
#include <algorithm>
#include <dune/common/dynvector.hh>
#include <dune/common/dynmatrix.hh>
 
namespace Dune
{
 
enum CurvatureLayout
{
  Vertex,
  Element
};
 
template<class IGridView, class IGridMapper, CurvatureLayout CL = Vertex>
class CurvatureOperator
{
public:
  static constexpr int dim = IGridView::dimensionworld;
  using Scalar = typename IGridView::ctype;
  using Vector = Dune::FieldVector<Scalar, dim>;
 
  CurvatureOperator(const IGridView& iGridView, const IGridMapper& iGridMapper)
   : iGridView_(iGridView), iGridMapper_(iGridMapper) {};
 
  template<class Curvatures, class Centers, CurvatureLayout Layout = CL>
  std::enable_if_t<Layout == Vertex, void>
  operator() (Curvatures& curvatures, Centers& centers) const
  {
    //storage for incident interface vertices
    std::vector<Vector> vertexPoints;
    Vector center;
 
    for (const auto& vertex : vertices(iGridView_))
    {
      const int iVertexIdx = iGridMapper_.index(vertex);
      vertexPoints.clear();
      vertexPoints.push_back(vertex.geometry().center());
 
      for (const auto& incidentVertex : incidentInterfaceVertices(vertex))
      {
        vertexPoints.push_back(incidentVertex.geometry().center());
      }
 
      //determine curvature associated with iVertexMapper (approximated by a
      //sphere with center "center")
      curvatures[iVertexIdx] = 1.0 / getRadius(vertexPoints, center);
      centers[iVertexIdx] = center;
    }
  }
 
  template<class Curvatures, class Centers, CurvatureLayout Layout = CL>
  std::enable_if_t<Layout == Element, void>
  operator() (Curvatures& curvatures, Centers& centers) const
  {
   //storage for vertices of incident interface elements
   std::vector<Vector> vertexPoints;
   Vector center;
 
   for (const auto& iElem : elements(iGridView_))
   {
     int iElemIdx = iGridMapper_.index(iElem);
     vertexPoints.clear();
 
     for (const auto& intersct : intersections(iGridView_, iElem))
     {
       const auto& neighborGeo = intersct.outside().geometry();
 
       for (int i = 0; i < neighborGeo.corners(); i++)
       {
         if (std::find(vertexPoints.begin(), vertexPoints.end(),
           neighborGeo.corner(i)) == vertexPoints.end())
         {
           vertexPoints.push_back(neighborGeo.corner(i));
         }
       }
     }
 
     //determine curvature associated with iElem (approximated by a sphere
     //with center "center")
     curvatures[iElemIdx] = 1.0 / getRadius(vertexPoints, center);
     centers[iElemIdx] = center;
   }
 }
 
private:
  template< class Points >
  double getRadius (const Points& points, Vector& center) const
  {
    using LargeVector = Dune::FieldVector<Scalar, dim+1>;
    using Matrix = Dune::FieldMatrix<Scalar, dim+1, dim+1>;
    using LargeDynVector = Dune::DynamicVector<Scalar>;
    using DynMatrix = Dune::DynamicMatrix<Scalar>;
 
    LargeDynVector b(points.size(), 0.0);
    DynMatrix A(points.size(), dim+1, 0.0);
    LargeVector x(0.0);
    LargeVector ATb(0.0); //A.transpose()*b
    Matrix ATA(0.0); //A.transpose()*A
 
    for (unsigned int i = 0; i < points.size(); i++)
    {
      b[i] = -points[i].two_norm2();
      A[i][0] = 1.0;
 
      for (int j = 0; j < dim; j++)
      {
        A[i][j+1] = points[i][j];
      }
    }
 
    A.mtv(b, ATb);
 
    for (int i = 0; i <= dim; i++)
      for (int j = 0; j <= dim; j++)
        for (std::size_t k = 0; k < points.size(); k ++)
          ATA[i][j] += A[k][i] * A[k][j];
 
    if (std::abs( ATA.determinant() ) > epsilon)
    {
      ATA.solve(x, ATb);
 
      double r = 0.0;
      for (int j = 1; j <= dim; j++)
      {
        r += x[j]*x[j];
        center[j-1] = -0.5*x[j];
      }
 
      r = sqrt(std::abs(0.25*r - x[0]));
 
      return r;
    }
 
    return std::numeric_limits<double>::infinity();
  }
 
  //The grid view of the interface grid
  const IGridView& iGridView_;
  //Element mapper for the interface grid
  const IGridMapper& iGridMapper_;
 
  //floating point accuracy
  static constexpr double epsilon = std::numeric_limits<double>::epsilon();
};
 
} // end namespace Dune
 
#endif