]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add no normal flux to MGConstrainedDoFs (box mesh) 7711/head
authortcclevenger <tcleven@clemson.edu>
Fri, 8 Feb 2019 22:42:12 +0000 (15:42 -0700)
committertcclevenger <tcleven@clemson.edu>
Thu, 28 Feb 2019 17:31:23 +0000 (10:31 -0700)
doc/news/changes/minor/20190211ConradClevenger [new file with mode: 0644]
include/deal.II/multigrid/mg_constrained_dofs.h
tests/multigrid/constrained_dofs_06.cc [new file with mode: 0644]
tests/multigrid/constrained_dofs_06.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190211ConradClevenger b/doc/news/changes/minor/20190211ConradClevenger
new file mode 100644 (file)
index 0000000..da07b58
--- /dev/null
@@ -0,0 +1,6 @@
+New: Function MGConstrainedDoFs::make_no_normal_flux_constraints which adds 
+functionality for no normal flux constraints during geometric mutigrid computations. 
+Currently, this function is limited to meshes with no normal flux boundaries 
+normal to the x-, y-, or z-axis.
+<br>
+(Conrad Clevenger, 2019/02/11)
index 40f2cd7c0d1c5e8831ad095528ed7a990ced17a2..452563eb3b115ab42f7c6ec8971cd3ae76302bb6 100644 (file)
@@ -100,6 +100,26 @@ public:
     const std::set<types::boundary_id> &boundary_ids,
     const ComponentMask &               component_mask = ComponentMask());
 
+  /**
+   * Fill the internal data structures with information
+   * about no normal flux boundary dofs.
+   *
+   * This function is limited to meshes whose no normal flux boundaries
+   * have faces which are normal to the x-, y-, or z-axis. Also, for a
+   * specific boundary id, all faces must be facing in the same direction,
+   * i.e., a boundary normal to the x-axis must have a different boundary
+   * id than a boundary normal to the y- or z-axis and so on. If the mesh
+   * was produced, for example, using the <tt>GridGenerator::hyper_cube()</tt>
+   * function, setting <tt>colorize=true</tt> during mesh generation and calling
+   * make_no_normal_flux_constraints() for each no normal flux boundary is
+   * sufficient.
+   */
+  template <int dim, int spacedim>
+  void
+  make_no_normal_flux_constraints(const DoFHandler<dim, spacedim> &dof,
+                                  const types::boundary_id         bid,
+                                  const unsigned int first_vector_component);
+
   /**
    * Reset the data structures.
    */
@@ -321,6 +341,61 @@ MGConstrainedDoFs::make_zero_boundary_constraints(
 }
 
 
+template <int dim, int spacedim>
+inline void
+MGConstrainedDoFs::make_no_normal_flux_constraints(
+  const DoFHandler<dim, spacedim> &dof,
+  const types::boundary_id         bid,
+  const unsigned int               first_vector_component)
+{
+  // For a given boundary id, find which vector component is on the boundary
+  // and set a zero boundary constraint for those degrees of freedom.
+  const unsigned int n_components = DoFTools::n_components(dof);
+  Assert(first_vector_component + dim <= n_components,
+         ExcIndexRange(first_vector_component, 0, n_components - dim + 1));
+
+  ComponentMask comp_mask(n_components, false);
+
+
+  typename Triangulation<dim>::face_iterator
+    face = dof.get_triangulation().begin_face(),
+    endf = dof.get_triangulation().end_face();
+  for (; face != endf; ++face)
+    if (face->at_boundary() && face->boundary_id() == bid)
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          Tensor<1, dim, double> unit_vec;
+          unit_vec[d] = 1.0;
+
+          const Tensor<1, dim> normal_vec =
+            face->get_manifold().normal_vector(face, face->center());
+
+          if (std::abs(std::abs(unit_vec * normal_vec) - 1.0) < 1e-10)
+            comp_mask.set(d + first_vector_component, true);
+          else
+            Assert(
+              std::abs(unit_vec * normal_vec) < 1e-10,
+              ExcMessage(
+                "We can currently only support no normal flux conditions "
+                "for a specific boundary id if all faces are normal to the "
+                "x, y, or z axis."));
+        }
+
+  Assert(comp_mask.n_selected_components() == 1,
+         ExcMessage(
+           "We can currently only support no normal flux conditions "
+           "for a specific boundary id if all faces are facing in the "
+           "same direction, i.e., a boundary normal to the x-axis must "
+           "have a different boundary id than a boundary normal to the "
+           "y- or z-axis and so on. If the mesh here was produced using "
+           "GridGenerator::..., setting colorize=true during mesh generation "
+           "and calling make_no_normal_flux_constraints() for each no normal "
+           "flux boundary will fulfill the condition."));
+
+  this->make_zero_boundary_constraints(dof, {bid}, comp_mask);
+}
+
+
 inline void
 MGConstrainedDoFs::clear()
 {
diff --git a/tests/multigrid/constrained_dofs_06.cc b/tests/multigrid/constrained_dofs_06.cc
new file mode 100644 (file)
index 0000000..cc0101f
--- /dev/null
@@ -0,0 +1,75 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Check make_no_normal_flux_constraints function for a box mesh
+
+#include <deal.II/base/exceptions.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> triangulation(
+    Triangulation<dim>::limit_level_difference_at_vertices);
+  GridGenerator::hyper_cube(triangulation, 0., 1., true);
+  triangulation.refine_global(4 - dim);
+
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(FESystem<dim>(FE_Q<dim>(1), dim));
+  dof_handler.distribute_mg_dofs();
+
+  std::set<types::boundary_id> no_flux_boundary = {0, 1, 2};
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof_handler);
+  for (auto bid : no_flux_boundary)
+    mg_constrained_dofs.make_no_normal_flux_constraints(dof_handler, bid, 0);
+
+  for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
+    {
+      deallog << "Level " << level << ": " << std::flush;
+      mg_constrained_dofs.get_boundary_indices(level).print(
+        deallog.get_file_stream());
+      deallog << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog << "2D" << std::endl;
+  test<2>();
+  deallog << std::endl << "3D" << std::endl;
+  test<3>();
+  return 0;
+}
diff --git a/tests/multigrid/constrained_dofs_06.output b/tests/multigrid/constrained_dofs_06.output
new file mode 100644 (file)
index 0000000..8f3fffe
--- /dev/null
@@ -0,0 +1,14 @@
+
+DEAL::2D
+DEAL::Level 0: {[0,4], 6}
+
+DEAL::Level 1: {[0,1], [3,4], [8,10], 12, 16}
+
+DEAL::Level 2: {[0,1], [3,4], 9, 12, 19, [22,24], 28, 30, 36, 44, 48}
+
+DEAL::
+DEAL::3D
+DEAL::Level 0: {[0,1], [3,4], 6, 9, [12,13], [15,16], 18, 21}
+
+DEAL::Level 1: {[0,1], 4, 6, [12,13], 16, 18, [24,25], 27, [30,31], 33, 36, 42, 48, 51, [54,55], 58, 60, [66,67], 69, 72, 78}
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.