From 0446a91f3e261648245bb39d519828340077a835 Mon Sep 17 00:00:00 2001
From: Melanie Gerault <mgerault@mit.edu>
Date: Thu, 8 Aug 2019 15:48:19 -0600
Subject: [PATCH] Add GridGenerator::eccentric_hyper_shell

Add an eccentric hyper shell geometry that is formed by two spheres
with different center points.
---
 source/grid/grid_generator.cc | 52 +++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 83ef28df08..c42eb2acbf 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -3003,6 +3003,58 @@ namespace GridGenerator
     tria.set_manifold(0, SphericalManifold<2>(center));
   }
 
+  template <int dim>
+  void
+  eccentric_hyper_shell(Triangulation<dim> &tria,
+                        const Point<dim> &  inner_center,
+                        const Point<dim> &  outer_center,
+                        const double        inner_radius,
+                        const double        outer_radius,
+                        const unsigned int  n_cells)
+  {
+    GridGenerator::hyper_shell(
+      tria, outer_center, inner_radius, outer_radius, n_cells, true);
+
+    // check the consistency of the dimensions provided
+    Assert(
+      outer_radius - inner_radius > outer_center.distance(inner_center),
+      ExcInternalError(
+        "Inner radius superior or equal to outer radius plus eccentricity."));
+
+    // shift nodes along the inner boundary according to the position of
+    // inner_circle
+    std::set<Point<dim> *> vertices_to_move;
+
+    for (const auto &cell : tria.active_cell_iterators())
+      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+        if (cell->face(f)->boundary_id() == 0)
+          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
+               ++v)
+            vertices_to_move.insert(&cell->face(f)->vertex(v));
+
+    const auto shift = inner_center - outer_center;
+    for (const auto &p : vertices_to_move)
+      (*p) += shift;
+
+    // the original hyper_shell function assigns the same manifold id
+    // to all cells and faces. Set all manifolds ids to a different
+    // value (2), then use boundary ids to assign different manifolds to
+    // the inner (0) and outer manifolds (1). Use a transfinite manifold
+    // for all faces and cells aside from the boundaries.
+    tria.set_all_manifold_ids(2);
+    GridTools::copy_boundary_to_manifold_id(tria);
+
+    SphericalManifold<dim> inner_manifold(inner_center);
+    SphericalManifold<dim> outer_manifold(outer_center);
+
+    TransfiniteInterpolationManifold<dim> transfinite;
+    transfinite.initialize(tria);
+
+    tria.set_manifold(0, inner_manifold);
+    tria.set_manifold(1, outer_manifold);
+    tria.set_manifold(2, transfinite);
+  }
+
 
 
   template <int dim>
-- 
2.39.5