--- /dev/null
+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)
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.
*/
}
+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()
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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}
+
projects / dealii.git / commitdiff