From 15b4cfaac6b2aba3036d05ac155ec31db62ef6c6 Mon Sep 17 00:00:00 2001 From: tcclevenger Date: Fri, 8 Feb 2019 15:42:12 -0700 Subject: [PATCH] add no normal flux to MGConstrainedDoFs (box mesh) --- .../changes/minor/20190211ConradClevenger | 6 ++ .../deal.II/multigrid/mg_constrained_dofs.h | 75 +++++++++++++++++++ tests/multigrid/constrained_dofs_06.cc | 75 +++++++++++++++++++ tests/multigrid/constrained_dofs_06.output | 14 ++++ 4 files changed, 170 insertions(+) create mode 100644 doc/news/changes/minor/20190211ConradClevenger create mode 100644 tests/multigrid/constrained_dofs_06.cc create mode 100644 tests/multigrid/constrained_dofs_06.output diff --git a/doc/news/changes/minor/20190211ConradClevenger b/doc/news/changes/minor/20190211ConradClevenger new file mode 100644 index 0000000000..da07b583e6 --- /dev/null +++ b/doc/news/changes/minor/20190211ConradClevenger @@ -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. +
+(Conrad Clevenger, 2019/02/11) diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 40f2cd7c0d..452563eb3b 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -100,6 +100,26 @@ public: const std::set &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 GridGenerator::hyper_cube() + * function, setting colorize=true during mesh generation and calling + * make_no_normal_flux_constraints() for each no normal flux boundary is + * sufficient. + */ + template + void + make_no_normal_flux_constraints(const DoFHandler &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 +inline void +MGConstrainedDoFs::make_no_normal_flux_constraints( + const DoFHandler &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::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 index 0000000000..cc0101f599 --- /dev/null +++ b/tests/multigrid/constrained_dofs_06.cc @@ -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 + +#include + +#include +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + +template +void +test() +{ + Triangulation triangulation( + Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_cube(triangulation, 0., 1., true); + triangulation.refine_global(4 - dim); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(FESystem(FE_Q(1), dim)); + dof_handler.distribute_mg_dofs(); + + std::set 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 index 0000000000..8f3fffee1a --- /dev/null +++ b/tests/multigrid/constrained_dofs_06.output @@ -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} + -- 2.39.5