makesphere.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_MAKE_SPHERE_HH
#define DUNE_FUFEM_MAKE_SPHERE_HH
 
#include <cmath>
#include <memory>
#include <type_traits>
 
#include <dune/common/hybridutilities.hh>
 
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/common/boundarysegment.hh>
 
class SphereQuadSegment : public Dune::BoundarySegment<3>
{
public:
    SphereQuadSegment(const Dune::FieldVector<double,3>& a,
                      const Dune::FieldVector<double,3>& b,
                      const Dune::FieldVector<double,3>& c,
                      const Dune::FieldVector<double,3>& d,
                      const Dune::FieldVector<double,3>& center,
                      double radius)
        : a_(a), b_(b), c_(c), d_(d), center_(center), radius_(radius)
    {
        assert( std::abs((a_ - center_).two_norm() - radius) < 1e-6 );
        assert( std::abs((b_ - center_).two_norm() - radius) < 1e-6 );
        assert( std::abs((c_ - center_).two_norm() - radius) < 1e-6 );
        assert( std::abs((d_ - center_).two_norm() - radius) < 1e-6 );
    }
 
    virtual Dune::FieldVector<double, 3> operator()(const Dune::FieldVector<double,2>& local) const {
 
        Dune::FieldVector<double, 3> result = a_ - center_;
        result.axpy(local[0], b_-a_);
        result.axpy(local[1], d_-a_);
        result.axpy(local[0]*local[1], a_ + c_ -b_ -d_);
 
        result *= radius_ / result.two_norm();
        result += center_;
 
        return result;
    }
 
    Dune::FieldVector<double, 3> a_, b_, c_, d_;
 
    Dune::FieldVector<double, 3> center_;
 
    double radius_;
};
 
class SphereTriSegment
  : public Dune::BoundarySegment<3>
{
public:
    SphereTriSegment(const Dune::FieldVector<double,3>& a,
                     const Dune::FieldVector<double,3>& b,
                     const Dune::FieldVector<double,3>& c,
                     const Dune::FieldVector<double,3>& center,
                     double radius)
        : a_(a), b_(b), c_(c), center_(center), radius_(radius)
    {
        assert( std::abs((a_ - center_).two_norm() - radius) < 1e-6 );
        assert( std::abs((b_ - center_).two_norm() - radius) < 1e-6 );
        assert( std::abs((c_ - center_).two_norm() - radius) < 1e-6 );
    }
 
    virtual Dune::FieldVector<double, 3> operator()(const Dune::FieldVector<double,2>& local) const {
 
        Dune::FieldVector<double, 3> result = a_ - center_;
        result.axpy(local[0],b_-a_);
        result.axpy(local[1],c_-a_);
 
        result *= radius_ / result.two_norm();
        result += center_;
 
        return result;
    }
 
    Dune::FieldVector<double, 3> a_, b_, c_;
 
    Dune::FieldVector<double, 3> center_;
 
    double radius_;
};
 
template <class GridType, class field_type = double>
std::unique_ptr<GridType> makeSphere(
    const Dune::FieldVector<field_type, 3>& center,
    field_type radius)
{
    using namespace Dune;
 
    Dune::GridFactory<GridType> factory;
 
    // ////////////////////////
    //   Insert vertices
    // ////////////////////////
    Dune::FieldVector<field_type,3> pos[8];
 
    double root = radius / std::sqrt(3);
 
    for (int i=0; i<8; i++) {
        pos[i] = center;
        pos[i][0] += (i&1) ? root : -root;
        pos[i][1] += (i&2) ? root : -root;
        pos[i][2] += (i&4) ? root : -root;
        factory.insertVertex(pos[i]);
    }
 
    factory.insertVertex(center);
 
    unsigned int segments[12][3] = {{2, 0, 6}, {6, 0, 4},
                                    {1, 3, 5}, {5, 3, 7},
                                    {0, 1, 4}, {4, 1, 5},
                                    {3, 2, 7}, {7, 2, 6},
                                    {1, 0, 3}, {3, 0, 2},
                                    {4, 5, 6}, {6, 5, 7}};
 
    std::vector<unsigned int> v(3);
    std::vector<unsigned int> cornerIDs(4);
 
    for (int i=0; i<12; i++) {
        for (int j=0; j<3; j++)
            v[j] = segments[i][j];
 
        auto boundarySegment = std::make_shared<SphereTriSegment>(pos[v[0]], pos[v[1]], pos[v[2]], center, radius);
        factory.insertBoundarySegment(v,boundarySegment);
 
        // /////////////////
        // Insert elements
        // /////////////////
 
        cornerIDs = {segments[i][0], segments[i][1], segments[i][2], 8};
        factory.insertElement(Dune::GeometryTypes::simplex(3), cornerIDs);
    }
 
    return std::unique_ptr<GridType>(factory.createGrid());
}
 
template <class GridType, class field_type = double>
std::unique_ptr<GridType> makeSphereOnOctahedron(const Dune::FieldVector<field_type, 3>& center,
                                 field_type radius) {
 
  static_assert(GridType::dimensionworld==3, "Only spheres in R^3 can be created!");
 
  Dune::GridFactory<GridType> factory;
 
  //   Insert vertices
 
  std::vector<Dune::FieldVector<double, 3> > pos = {{{1, 0, 0}}, {{0, 1, 0}}, {{-1, 0, 0}},
                                                    {{0, -1, 0}}, {{0, 0, 1}}, {{0, 0, -1}}};
 
  for (int i=0; i<6; i++)
  {
    pos[i] *= radius;
    pos[i] += center;
    factory.insertVertex(pos[i]);
  }
 
  // 3d grids get a vertex at the sphere center
  if (GridType::dimension==3)
    factory.insertVertex(center);
 
  // Insert elements
  std::vector<std::vector<unsigned int> > segments = {{0, 1, 4}, {1, 2, 4},
                                                      {2, 3, 4}, {3, 0, 4},
                                                      {1, 0, 5}, {2, 1, 5},
                                                      {3, 2, 5}, {0, 3, 5}};
 
  std::vector<unsigned int> cornerIDs;
 
  for (int i=0; i<8; i++)
  {
    // If the grid is sd the element boundaries carry a parametrization
    Dune::Hybrid::ifElse(std::bool_constant<GridType::dimension==3>{},
    [&](auto id)
    {
        auto boundarySegment = std::make_shared<SphereTriSegment>(pos[segments[i][0]],
                                                                  pos[segments[i][1]],
                                                                  pos[segments[i][2]], center, radius);
        id(factory).insertBoundarySegment(segments[i], id(boundarySegment));
 
    });
 
    // Insert elements
    if (GridType::dimension==2)
      cornerIDs = {segments[i][0], segments[i][1], segments[i][2]};
    else
      cornerIDs = {segments[i][0], segments[i][1], segments[i][2], 6};
 
    Dune::Hybrid::ifElse(std::bool_constant<GridType::dimension==2>{},
    [&](auto id)
    {
      // Insert parametrized element
      auto elementParametrization = std::make_shared<SphereTriSegment>(pos[segments[i][0]],
                                                                       pos[segments[i][1]],
                                                                       pos[segments[i][2]],
                                                                       center,
                                                                       radius);
      id(factory).insertElement(Dune::GeometryTypes::simplex(GridType::dimension), cornerIDs, elementParametrization);
    },
    [&](auto id)
    {
      id(factory).insertElement(Dune::GeometryTypes::simplex(GridType::dimension), cornerIDs);
    });
  }
  return std::unique_ptr<GridType>(factory.createGrid());
}
 
#endif