Merge pull request #31 from fdrmrc/fix_2d3

Fix some 2d3d tests

Author: fdrmrc

README.md

@@ -1,18 +1,16 @@
 # Non Matching Test Suite
 
-Collection of C++ finite element programs testing different coupling strategies with Non Matching techniques using composite quadrature rules on mesh intersections. This is done by integrating the C++ library `CGAL` (https://www.cgal.org/) into `deal.II` (www.dealii.org). Currently, all examples require a `deal.II` version greater or equal to `9.4.2`, a C++17-compliant compiler and a `CGAL` version greater than 5.5.2 installed.
+C++ application consisting of finite element programs testing different coupling strategies with via non matching techniques such as using quadrature rules on mesh intersections. Part of this work builds on top of the integration of the C++ library `CGAL` (https://www.cgal.org/) into `deal.II` (www.dealii.org). Currently, a `deal.II` version greater or equal than `9.4.2`, and `CGAL` versions greater than 5.5.2 are required, along with a C++17 compliant compiler.
 
-![Screenshot](grids/Flower_interface.png)
+![Screenshot](grids/mesh_3d.png)
 
-## How to compile and run
-We provide parameter files (.prm) for test cases, allowing the user to change rhs, boundary conditions, number of refinement cycles and mesh related parameters. To compile and run, move into the desired folder (e.g. `/LagrangeMultiplier/1d2d`).
-
-- Set the `DEAL_II_DIR` to the directory where you built the branch above. Alternatively, you can pass it as one of the CMake flags
+![Screenshot](grids/iso_contour_3D_ns.png)
 
-- `cmake .` , or `cmake -DDEAL_II_DIR=/your/local/path/to/deal` if you want to pass it as a flag
+## How to compile and run
 
-- Compile with `make ` (or `make -jN`)
+- `mkdir build && cd build`
+- `cmake .` , or `cmake -DDEAL_II_DIR=/your/local/path/to/deal` 
+- `make`, or `make -j N`, begin `N` the number of make jobs you may want to ask.
 
-- Run `./lagrange_multiplier disk_parameters.prm`
 
 ## Using Docker image [TODO]

examples/interface_penalisation/2d3d/interface_penalisation_non_smooth_disk.cc

@@ -91,11 +91,9 @@ template <>
 double Solution<3>::value(const Point<3> &p,
                           const unsigned int component) const {
   (void)component;
-  // const Point<3> xc{Cx, Cy, Cz}; // center of the sphere
-  // const double   r = (p - xc).norm();
-  // return r <= R ? 1. / R : 1. / r;
-  return std::sin(2. * numbers::PI * p[0]) * std::sin(2. * numbers::PI * p[1]) *
-         std::sin(2. * numbers::PI * p[2]);
+  const Point<3> xc{Cx, Cy, Cz}; // center of the sphere
+  const double r = (p - xc).norm();
+  return r <= R ? 1. / R : 1. / r;
 }
 
 template <>
@@ -120,6 +118,8 @@ template <int dim, int spacedim> class PoissonNitscheInterface {
 private:
   void generate_grids();
 
+  void adjust_grids();
+
   void setup_system();
 
   void assemble_system();
@@ -211,6 +211,112 @@ void PoissonNitscheInterface<dim, spacedim>::generate_grids() {
       std::make_unique<GridTools::Cache<dim, spacedim>>(embedded_triangulation);
 }
 
+template <int dim, int spacedim>
+void PoissonNitscheInterface<dim, spacedim>::adjust_grids() {
+  // Adjust grid diameters to satisfy ratio suggested by the theory.
+
+  std::cout << "Adjusting the grids..." << std::endl;
+  namespace bgi = boost::geometry::index;
+
+  auto refine = [&]() {
+    bool done = false;
+
+    double min_embedded = 1e10;
+    double max_embedded = 0;
+    double min_space = 1e10;
+    double max_space = 0;
+
+    while (done == false) {
+      // Bounding boxes of the space grid
+      const auto &tree =
+          space_cache->get_locally_owned_cell_bounding_boxes_rtree();
+
+      // Bounding boxes of the embedded grid
+      const auto &embedded_tree =
+          embedded_cache->get_cell_bounding_boxes_rtree();
+
+      // Let's check all cells whose bounding box contains an embedded
+      // bounding box
+      done = true;
+
+      const bool use_space = false;
+
+      const bool use_embedded = false;
+
+      AssertThrow(!(use_embedded && use_space),
+                  ExcMessage("You can't refine both the embedded and "
+                             "the space grid at the same time."));
+
+      for (const auto &[embedded_box, embedded_cell] : embedded_tree) {
+        const auto &[p1, p2] = embedded_box.get_boundary_points();
+        const auto diameter = p1.distance(p2);
+        min_embedded = std::min(min_embedded, diameter);
+        max_embedded = std::max(max_embedded, diameter);
+
+        for (const auto &[space_box, space_cell] :
+             tree | bgi::adaptors::queried(bgi::intersects(embedded_box))) {
+          const auto &[sp1, sp2] = space_box.get_boundary_points();
+          const auto space_diameter = sp1.distance(sp2);
+          min_space = std::min(min_space, space_diameter);
+          max_space = std::max(max_space, space_diameter);
+
+          if (use_embedded && space_diameter < diameter) {
+            embedded_cell->set_refine_flag();
+            done = false;
+          }
+          if (use_space && diameter < space_diameter) {
+            space_cell->set_refine_flag();
+            done = false;
+          }
+        }
+      }
+      if (done == false) {
+        if (use_embedded) {
+          // Compute again the embedded displacement grid
+          embedded_triangulation.execute_coarsening_and_refinement();
+        }
+        if (use_space) {
+          // Compute again the embedded displacement grid
+          space_triangulation.execute_coarsening_and_refinement();
+        }
+      }
+    }
+    return std::make_tuple(min_space, max_space, min_embedded, max_embedded);
+  };
+
+  // Do the refinement loop once, to make sure we satisfy our criterions
+  refine();
+
+  // Pre refine the space grid according to the delta refinement
+  const bool apply_delta_refinements = true;
+  unsigned int space_pre_refinement_cycles = 1;
+  if (apply_delta_refinements)
+    for (unsigned int i = 0; i < space_pre_refinement_cycles; ++i) {
+      const auto &tree =
+          space_cache->get_locally_owned_cell_bounding_boxes_rtree();
+
+      const auto &embedded_tree =
+          embedded_cache->get_cell_bounding_boxes_rtree();
+
+      for (const auto &[embedded_box, embedded_cell] : embedded_tree)
+        for (const auto &[space_box, space_cell] :
+             tree | bgi::adaptors::queried(bgi::intersects(embedded_box)))
+          space_cell->set_refine_flag();
+      space_triangulation.execute_coarsening_and_refinement();
+
+      // Make sure again we satisfy our criterion after the space refinement
+      refine();
+    }
+
+  // Check once again we satisfy our criterion, and record min/max
+  const auto [sm, sM, em, eM] = refine();
+
+  std::cout << "Space local min/max diameters   : " << sm << "/" << sM
+            << std::endl
+            << "Embedded space min/max diameters: " << em << "/" << eM
+            << std::endl;
+}
+
 template <int dim, int spacedim>
 void PoissonNitscheInterface<dim, spacedim>::setup_system() {
   TimerOutput::Scope timer_section(timer, "Setup system");
@@ -357,15 +463,18 @@ template <int dim, int spacedim>
 void PoissonNitscheInterface<dim, spacedim>::output_results(
     const unsigned cycle) const {
   TimerOutput::Scope timer_section(timer, "Output results");
-  // std::cout << "Output results" << std::endl;
-  data_out.clear();
-  data_out.attach_dof_handler(space_dh);
-  data_out.add_data_vector(solution, "solution");
-  data_out.build_patches();
-  std::ofstream output("solution_nitsche" + std::to_string(dim) +
-                       std::to_string(spacedim) + std::to_string(cycle) +
-                       ".vtu");
-  data_out.write_vtu(output);
+
+  if (cycle < 2) {
+    data_out.clear();
+    data_out.attach_dof_handler(space_dh);
+    data_out.add_data_vector(solution, "solution");
+    data_out.build_patches();
+    std::ofstream output("solution_nitsche" + std::to_string(dim) +
+                         std::to_string(spacedim) + std::to_string(cycle) +
+                         ".vtu");
+    data_out.write_vtu(output);
+  }
+
   {
     Vector<double> difference_per_cell(space_triangulation.n_active_cells());
     VectorTools::integrate_difference(
@@ -385,19 +494,11 @@ void PoissonNitscheInterface<dim, spacedim>::output_results(
 
     convergence_table.add_value("cycle", cycle);
     convergence_table.add_value("cells", space_triangulation.n_active_cells());
-    convergence_table.add_value("dofs", space_dh.n_dofs());
+    convergence_table.add_value("dofs", space_dh.n_dofs() -
+                                            space_constraints.n_constraints());
     convergence_table.add_value("L2", L2_error);
     convergence_table.add_value("H1", H1_error);
   }
-
-  {
-    if (cycle == 3) {
-      std::ofstream output_test_space("space_grid_cycle3.vtk");
-      GridOut().write_vtk(space_triangulation, output_test_space);
-      std::ofstream output_test_embedded("embedded_grid_cycle3.vtk");
-      GridOut().write_vtk(embedded_triangulation, output_test_embedded);
-    }
-  }
 }
 
 // The run() method here differs only in the call to
@@ -407,7 +508,8 @@ void PoissonNitscheInterface<dim, spacedim>::run() {
   generate_grids();
   for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle) {
     std::cout << "Cycle: " << cycle << std::endl;
-
+    if (cycle == 0u)
+      adjust_grids();
     // Compute all the things we need to assemble the Nitsche's
     // contributions, namely the two cached triangulations and a degree to
     // integrate over the intersections.

examples/interface_penalisation/2d3d/interface_penalisation_smooth_disk.cc

@@ -124,6 +124,8 @@ template <int dim, int spacedim> class PoissonNitscheInterface {
 private:
   void generate_grids();
 
+  void adjust_grids();
+
   void setup_system();
 
   void assemble_system();
@@ -216,6 +218,112 @@ void PoissonNitscheInterface<dim, spacedim>::generate_grids() {
       std::make_unique<GridTools::Cache<dim, spacedim>>(embedded_triangulation);
 }
 
+template <int dim, int spacedim>
+void PoissonNitscheInterface<dim, spacedim>::adjust_grids() {
+  // Adjust grid diameters to satisfy ratio suggested by the theory.
+
+  std::cout << "Adjusting the grids..." << std::endl;
+  namespace bgi = boost::geometry::index;
+
+  auto refine = [&]() {
+    bool done = false;
+
+    double min_embedded = 1e10;
+    double max_embedded = 0;
+    double min_space = 1e10;
+    double max_space = 0;
+
+    while (done == false) {
+      // Bounding boxes of the space grid
+      const auto &tree =
+          space_cache->get_locally_owned_cell_bounding_boxes_rtree();
+
+      // Bounding boxes of the embedded grid
+      const auto &embedded_tree =
+          embedded_cache->get_cell_bounding_boxes_rtree();
+
+      // Let's check all cells whose bounding box contains an embedded
+      // bounding box
+      done = true;
+
+      const bool use_space = false;
+
+      const bool use_embedded = false;
+
+      AssertThrow(!(use_embedded && use_space),
+                  ExcMessage("You can't refine both the embedded and "
+                             "the space grid at the same time."));
+
+      for (const auto &[embedded_box, embedded_cell] : embedded_tree) {
+        const auto &[p1, p2] = embedded_box.get_boundary_points();
+        const auto diameter = p1.distance(p2);
+        min_embedded = std::min(min_embedded, diameter);
+        max_embedded = std::max(max_embedded, diameter);
+
+        for (const auto &[space_box, space_cell] :
+             tree | bgi::adaptors::queried(bgi::intersects(embedded_box))) {
+          const auto &[sp1, sp2] = space_box.get_boundary_points();
+          const auto space_diameter = sp1.distance(sp2);
+          min_space = std::min(min_space, space_diameter);
+          max_space = std::max(max_space, space_diameter);
+
+          if (use_embedded && space_diameter < diameter) {
+            embedded_cell->set_refine_flag();
+            done = false;
+          }
+          if (use_space && diameter < space_diameter) {
+            space_cell->set_refine_flag();
+            done = false;
+          }
+        }
+      }
+      if (done == false) {
+        if (use_embedded) {
+          // Compute again the embedded displacement grid
+          embedded_triangulation.execute_coarsening_and_refinement();
+        }
+        if (use_space) {
+          // Compute again the embedded displacement grid
+          space_triangulation.execute_coarsening_and_refinement();
+        }
+      }
+    }
+    return std::make_tuple(min_space, max_space, min_embedded, max_embedded);
+  };
+
+  // Do the refinement loop once, to make sure we satisfy our criterions
+  refine();
+
+  // Pre refine the space grid according to the delta refinement
+  const bool apply_delta_refinements = true;
+  unsigned int space_pre_refinement_cycles = 1;
+  if (apply_delta_refinements)
+    for (unsigned int i = 0; i < space_pre_refinement_cycles; ++i) {
+      const auto &tree =
+          space_cache->get_locally_owned_cell_bounding_boxes_rtree();
+
+      const auto &embedded_tree =
+          embedded_cache->get_cell_bounding_boxes_rtree();
+
+      for (const auto &[embedded_box, embedded_cell] : embedded_tree)
+        for (const auto &[space_box, space_cell] :
+             tree | bgi::adaptors::queried(bgi::intersects(embedded_box)))
+          space_cell->set_refine_flag();
+      space_triangulation.execute_coarsening_and_refinement();
+
+      // Make sure again we satisfy our criterion after the space refinement
+      refine();
+    }
+
+  // Check once again we satisfy our criterion, and record min/max
+  const auto [sm, sM, em, eM] = refine();
+
+  std::cout << "Space local min/max diameters   : " << sm << "/" << sM
+            << std::endl
+            << "Embedded space min/max diameters: " << em << "/" << eM
+            << std::endl;
+}
+
 template <int dim, int spacedim>
 void PoissonNitscheInterface<dim, spacedim>::setup_system() {
   TimerOutput::Scope timer_section(timer, "Setup system");
@@ -388,7 +496,8 @@ void PoissonNitscheInterface<dim, spacedim>::output_results(
 
     convergence_table.add_value("cycle", cycle);
     convergence_table.add_value("cells", space_triangulation.n_active_cells());
-    convergence_table.add_value("dofs", space_dh.n_dofs());
+    convergence_table.add_value("dofs", space_dh.n_dofs() -
+                                            space_constraints.n_constraints());
     convergence_table.add_value("L2", L2_error);
     convergence_table.add_value("H1", H1_error);
   }
@@ -411,6 +520,9 @@ void PoissonNitscheInterface<dim, spacedim>::run() {
   for (unsigned int cycle = 0; cycle < n_refinement_cycles; ++cycle) {
     std::cout << "Cycle: " << cycle << std::endl;
 
+    if (cycle == 0u) {
+      adjust_grids();
+    }
     // Compute all the things we need to assemble the Nitsche's
     // contributions, namely the two cached triangulations and a degree to
     // integrate over the intersections.
@@ -437,7 +549,6 @@ void PoissonNitscheInterface<dim, spacedim>::run() {
       solve();
     }
 
-    // error_table.error_from_exact(space_dh, solution, exact_solution);
     output_results(cycle);
 
     if (cycle < n_refinement_cycles - 1) {