mapping);
}
+ template <int dim, int spacedim>
+ void
+ compute_no_normal_flux_constraints_on_level(
+ const DoFHandler<dim, spacedim> & dof_handler,
+ const MGConstrainedDoFs & mg_constrained_dofs,
+ const unsigned int level,
+ const unsigned int first_vector_component,
+ const std::set<types::boundary_id> &boundary_ids,
+ AffineConstraints<double> & constraints,
+ const Mapping<dim> & mapping)
+ {
+ // Copied from compute_nonzero_normal_flux_constraints()
+ // Only changed active_cell_iterator to cell_iterator in DoFToNormalsMap,
+ // CellToNormalsMap, CellContributions, and set inhomogeneity to 0.
+ const IndexSet &refinement_edge_indices =
+ mg_constrained_dofs.get_refinement_edge_indices(level);
+
+ const auto & fe = dof_handler.get_fe();
+ const std::vector<Point<dim - 1>> &unit_support_points =
+ fe.get_unit_face_support_points();
+ const Quadrature<dim - 1> quadrature(unit_support_points);
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+ std::vector<types::global_dof_index> face_dofs(dofs_per_face);
+
+
+ FEFaceValues<dim, spacedim> fe_face_values(mapping,
+ fe,
+ quadrature,
+ update_quadrature_points |
+ update_normal_vectors);
+
+ std::set<types::boundary_id>::iterator b_id;
+ using DoFToNormalsMap = std::multimap<
+ internal::VectorDoFTuple<dim>,
+ std::pair<Tensor<1, dim>,
+ typename DoFHandler<dim, spacedim>::cell_iterator>>;
+ DoFToNormalsMap dof_to_normals_map;
+ for (const auto &cell : dof_handler.cell_iterators_on_level(level))
+ if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
+ cell->level_subdomain_id() != numbers::invalid_subdomain_id)
+ for (const unsigned int face_no : cell->face_indices())
+ if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) !=
+ boundary_ids.end())
+ {
+ typename DoFHandler<dim, spacedim>::level_face_iterator face =
+ cell->face(face_no);
+ face->get_mg_dof_indices(level, face_dofs);
+ fe_face_values.reinit(cell, face_no);
+
+ for (unsigned int i = 0; i < face_dofs.size(); ++i)
+ if (fe.face_system_to_component_index(i).first ==
+ first_vector_component)
+ // Refinement edge indices are going to be constrained to 0
+ // during a multigrid cycle and do not need no-normal-flux
+ // constraints, so skip them:
+ if (!refinement_edge_indices.is_element(face_dofs[i]))
+ {
+ const Point<dim> position =
+ fe_face_values.quadrature_point(i);
+ internal::VectorDoFTuple<dim> vector_dofs;
+ vector_dofs.dof_indices[0] = face_dofs[i];
+ for (unsigned int k = 0; k < dofs_per_face; ++k)
+ if ((k != i) &&
+ (quadrature.point(k) == quadrature.point(i)) &&
+ (fe.face_system_to_component_index(k).first >=
+ first_vector_component) &&
+ (fe.face_system_to_component_index(k).first <
+ first_vector_component + dim))
+ vector_dofs.dof_indices
+ [fe.face_system_to_component_index(k).first -
+ first_vector_component] = face_dofs[k];
+
+ Tensor<1, dim> normal_vector =
+ cell->face(face_no)->get_manifold().normal_vector(
+ cell->face(face_no), position);
+
+ // make sure the normal vector is pointing to the right
+ // direction, more dietails can be found in
+ // compute_nonzero_normal_flux_constraints()
+ if (normal_vector * fe_face_values.normal_vector(i) < 0)
+ normal_vector *= -1;
+ Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
+ ExcInternalError());
+
+ // remove small entries:
+ for (unsigned int d = 0; d < dim; ++d)
+ if (std::fabs(normal_vector[d]) < 1e-13)
+ normal_vector[d] = 0;
+ normal_vector /= normal_vector.norm();
+
+ dof_to_normals_map.insert(
+ std::make_pair(vector_dofs,
+ std::make_pair(normal_vector, cell)));
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << "Adding normal vector:" << std::endl
+ << " dofs=" << vector_dofs << std::endl
+ << " cell=" << cell << " at "
+ << cell->center() << std::endl
+ << " normal=" << normal_vector << std::endl;
+#endif
+ }
+ }
+ // Now do something with the collected information. To this end, loop
+ // through all sets of pairs (dofs,normal_vector) and identify which
+ // entries belong to the same set of dofs and then do as described in the
+ // documentation, i.e. either average the normal vector or don't for this
+ // particular set of dofs
+ typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
+
+ while (p != dof_to_normals_map.end())
+ {
+ // first find the range of entries in the multimap that corresponds to
+ // the same vector-dof tuple. as usual, we define the range
+ // half-open. the first entry of course is 'p'
+ typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
+ for (++p; p != dof_to_normals_map.end(); ++p)
+ if (p->first != same_dof_range[0]->first)
+ {
+ same_dof_range[1] = p;
+ break;
+ }
+ if (p == dof_to_normals_map.end())
+ same_dof_range[1] = dof_to_normals_map.end();
+
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << "For dof indices <" << p->first
+ << ">, found the following normals" << std::endl;
+ for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+ q != same_dof_range[1];
+ ++q)
+ std::cout << " " << q->second.first << " from cell "
+ << q->second.second << std::endl;
+#endif
+
+
+ // now compute the reverse mapping: for each of the cells that
+ // contributed to the current set of vector dofs, add up the normal
+ // vectors. the values of the map are pairs of normal vectors and
+ // number of cells that have contributed
+ using CellToNormalsMap =
+ std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+ std::pair<Tensor<1, dim>, unsigned int>>;
+
+ CellToNormalsMap cell_to_normals_map;
+ for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+ q != same_dof_range[1];
+ ++q)
+ if (cell_to_normals_map.find(q->second.second) ==
+ cell_to_normals_map.end())
+ cell_to_normals_map[q->second.second] =
+ std::make_pair(q->second.first, 1U);
+ else
+ {
+ const Tensor<1, dim> old_normal =
+ cell_to_normals_map[q->second.second].first;
+ const unsigned int old_count =
+ cell_to_normals_map[q->second.second].second;
+
+ Assert(old_count > 0, ExcInternalError());
+
+ // in the same entry, store again the now averaged normal vector
+ // and the new count
+ cell_to_normals_map[q->second.second] =
+ std::make_pair((old_normal * old_count + q->second.first) /
+ (old_count + 1),
+ old_count + 1);
+ }
+ Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
+
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout << " cell_to_normals_map:" << std::endl;
+ for (typename CellToNormalsMap::const_iterator x =
+ cell_to_normals_map.begin();
+ x != cell_to_normals_map.end();
+ ++x)
+ std::cout << " " << x->first << " -> (" << x->second.first << ','
+ << x->second.second << ')' << std::endl;
+#endif
+
+ // count the maximum number of contributions from each cell
+ unsigned int max_n_contributions_per_cell = 1;
+ for (typename CellToNormalsMap::const_iterator x =
+ cell_to_normals_map.begin();
+ x != cell_to_normals_map.end();
+ ++x)
+ max_n_contributions_per_cell =
+ std::max(max_n_contributions_per_cell, x->second.second);
+
+ // verify that each cell can have only contributed at most dim times,
+ // since that is the maximum number of faces that come together at a
+ // single place
+ Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
+
+ switch (max_n_contributions_per_cell)
+ {
+ case 1:
+ {
+ // compute the average normal vector from all the ones that
+ // have the same set of dofs. we could add them up and divide
+ // them by the number of additions, or simply normalize them
+ // right away since we want them to have unit length anyway
+ Tensor<1, dim> normal;
+ for (typename CellToNormalsMap::const_iterator x =
+ cell_to_normals_map.begin();
+ x != cell_to_normals_map.end();
+ ++x)
+ normal += x->second.first;
+ normal /= normal.norm();
+
+ // normalize again
+ for (unsigned int d = 0; d < dim; ++d)
+ if (std::fabs(normal[d]) < 1e-13)
+ normal[d] = 0;
+ normal /= normal.norm();
+
+ const auto &dof_indices = same_dof_range[0]->first;
+ internal::add_constraint<dim>(dof_indices,
+ normal,
+ constraints,
+ 0.0);
+
+ break;
+ }
+
+ case dim:
+ {
+ // assert that indeed only a single cell has contributed
+ Assert(cell_to_normals_map.size() == 1, ExcInternalError());
+
+ // check linear independence by computing the determinant of
+ // the matrix created from all the normal vectors. if they are
+ // linearly independent, then the determinant is nonzero. if
+ // they are orthogonal, then the matrix is in fact equal to 1
+ // (since they are all unit vectors); make sure the
+ // determinant is larger than 1e-3 to avoid cases where cells
+ // are degenerate
+ {
+ Tensor<2, dim> t;
+
+ typename DoFToNormalsMap::const_iterator x =
+ same_dof_range[0];
+ for (unsigned int i = 0; i < dim; ++i, ++x)
+ for (unsigned int j = 0; j < dim; ++j)
+ t[i][j] = x->second.first[j];
+
+ Assert(
+ std::fabs(determinant(t)) > 1e-3,
+ ExcMessage(
+ "Found a set of normal vectors that are nearly collinear."));
+ }
+
+ // so all components of this vector dof are constrained. enter
+ // this into the AffineConstraints object
+ //
+ // ignore dofs already constrained
+ const auto &dof_indices = same_dof_range[0]->first;
+
+ for (unsigned int i = 0; i < dim; ++i)
+ if (!constraints.is_constrained(
+ same_dof_range[0]->first.dof_indices[i]) &&
+ constraints.can_store_line(
+ same_dof_range[0]->first.dof_indices[i]))
+ {
+ const types::global_dof_index line =
+ dof_indices.dof_indices[i];
+ constraints.add_line(line);
+ }
+
+ break;
+ }
+
+ // this is the case of an edge contribution in 3d, i.e. the vector
+ // is constrained in two directions but not the third.
+ default:
+ {
+ Assert(dim >= 3, ExcNotImplemented());
+ Assert(max_n_contributions_per_cell == 2, ExcInternalError());
+
+ // as described in the documentation, let us first collect
+ // what each of the cells contributed at the current point. we
+ // use a std::list instead of a std::set (which would be more
+ // natural) because std::set requires that the stored elements
+ // are comparable with operator<
+ using CellContributions =
+ std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+ std::list<Tensor<1, dim>>>;
+ CellContributions cell_contributions;
+
+ for (typename DoFToNormalsMap::const_iterator q =
+ same_dof_range[0];
+ q != same_dof_range[1];
+ ++q)
+ cell_contributions[q->second.second].push_back(
+ q->second.first);
+ Assert(cell_contributions.size() >= 1, ExcInternalError());
+
+ // now for each cell that has contributed determine the number
+ // of normal vectors it has contributed. we currently only
+ // implement if this is dim-1 for all cells (if a single cell
+ // has contributed dim, or if all adjacent cells have
+ // contributed 1 normal vector, this is already handled
+ // above).
+ //
+ // we only implement the case that all cells contribute
+ // dim-1 because we assume that we are following an edge
+ // of the domain (think: we are looking at a vertex
+ // located on one of the edges of a refined cube where the
+ // boundary indicators of the two adjacent faces of the
+ // cube are both listed in the set of boundary indicators
+ // passed to this function). in that case, all cells along
+ // that edge of the domain are assumed to have contributed
+ // dim-1 normal vectors. however, there are cases where
+ // this assumption is not justified (see the lengthy
+ // explanation in test no_flux_12.cc) and in those cases
+ // we simply ignore the cell that contributes only
+ // once. this is also discussed at length in the
+ // documentation of this function.
+ //
+ // for each contributing cell compute the tangential vector
+ // that remains unconstrained
+ std::list<Tensor<1, dim>> tangential_vectors;
+ for (typename CellContributions::const_iterator contribution =
+ cell_contributions.begin();
+ contribution != cell_contributions.end();
+ ++contribution)
+ {
+#ifdef DEBUG_NO_NORMAL_FLUX
+ std::cout
+ << " Treating edge case with dim-1 contributions."
+ << std::endl
+ << " Looking at cell " << contribution->first
+ << " which has contributed these normal vectors:"
+ << std::endl;
+ for (typename std::list<Tensor<1, dim>>::const_iterator t =
+ contribution->second.begin();
+ t != contribution->second.end();
+ ++t)
+ std::cout << " " << *t << std::endl;
+#endif
+
+ // as mentioned above, simply ignore cells that only
+ // contribute once
+ if (contribution->second.size() < dim - 1)
+ continue;
+
+ Tensor<1, dim> normals[dim - 1];
+ {
+ unsigned int index = 0;
+ for (typename std::list<Tensor<1, dim>>::const_iterator
+ t = contribution->second.begin();
+ t != contribution->second.end();
+ ++t, ++index)
+ normals[index] = *t;
+ Assert(index == dim - 1, ExcInternalError());
+ }
+
+ // calculate the tangent as the outer product of the
+ // normal vectors. since these vectors do not need to be
+ // orthogonal (think, for example, the case of the
+ // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
+ // elements, where we have constraints from the two normal
+ // vectors of two faces of the sheared cube that are not
+ // perpendicular to each other), we have to normalize the
+ // outer product
+ Tensor<1, dim> tangent;
+ switch (dim)
+ {
+ case 3:
+ // take cross product between normals[0] and
+ // normals[1]. write it in the current form (with
+ // [dim-2]) to make sure that compilers don't warn
+ // about out-of-bounds accesses -- the warnings are
+ // bogus since we get here only for dim==3, but at
+ // least one isn't quite smart enough to notice this
+ // and warns when compiling the function in 2d
+ tangent =
+ cross_product_3d(normals[0], normals[dim - 2]);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+
+ Assert(
+ std::fabs(tangent.norm()) > 1e-12,
+ ExcMessage(
+ "Two normal vectors from adjacent faces are almost "
+ "parallel."));
+ tangent /= tangent.norm();
+
+ tangential_vectors.push_back(tangent);
+ }
+
+ // go through the list of tangents and make sure that they all
+ // roughly point in the same direction as the first one (i.e.
+ // have an angle less than 90 degrees); if they don't then
+ // flip their sign
+ {
+ const Tensor<1, dim> first_tangent =
+ tangential_vectors.front();
+ typename std::list<Tensor<1, dim>>::iterator t =
+ tangential_vectors.begin();
+ ++t;
+ for (; t != tangential_vectors.end(); ++t)
+ if (*t * first_tangent < 0)
+ *t *= -1;
+ }
+
+ // now compute the average tangent and normalize it
+ Tensor<1, dim> average_tangent;
+ for (typename std::list<Tensor<1, dim>>::const_iterator t =
+ tangential_vectors.begin();
+ t != tangential_vectors.end();
+ ++t)
+ average_tangent += *t;
+ average_tangent /= average_tangent.norm();
+
+ const auto &dof_indices = same_dof_range[0]->first;
+ internal::add_tangentiality_constraints(dof_indices,
+ average_tangent,
+ constraints);
+ }
+ }
+ }
+ }
+
+
+
template <int dim, int spacedim>
void
compute_normal_flux_constraints(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// add tests for compute_no_normal_flux_constraints_on_level()
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.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/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_constraints.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+run(const Triangulation<dim> & triangulation,
+ const std::set<types::boundary_id> &no_flux_boundary)
+{
+ FESystem<dim> fe(FE_Q<dim>(1), dim);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+ dof_handler.distribute_mg_dofs();
+
+ MGConstrainedDoFs mg_constrained_dofs;
+ mg_constrained_dofs.initialize(dof_handler);
+
+ const unsigned int n_levels = triangulation.n_global_levels();
+
+ MappingQ<dim> mapping(4);
+
+ for (unsigned int level = 0; level < n_levels; ++level)
+ {
+ AffineConstraints<double> user_level_constraints;
+
+ VectorTools::compute_no_normal_flux_constraints_on_level(
+ dof_handler,
+ mg_constrained_dofs,
+ level,
+ 0,
+ no_flux_boundary,
+ user_level_constraints,
+ mapping);
+
+ user_level_constraints.print(deallog.get_file_stream());
+
+ deallog.get_file_stream() << std::flush;
+ user_level_constraints.close();
+
+ deallog << "Level " << level << " OK" << std::endl;
+ }
+}
+
+
+int
+main()
+{
+ initlog();
+ deallog.get_file_stream().precision(7);
+ deallog.get_file_stream().setf(std::ios::fixed);
+
+ {
+ const unsigned int dim = 2;
+ Triangulation<dim> triangulation(
+ Triangulation<dim>::limit_level_difference_at_vertices);
+ GridGenerator::quarter_hyper_ball(triangulation);
+ triangulation.refine_global(1);
+ std::set<types::boundary_id> no_flux_boundary{0, 1};
+ run<dim>(triangulation, no_flux_boundary);
+ }
+ {
+ const unsigned int dim = 3;
+ Triangulation<dim> triangulation(
+ Triangulation<dim>::limit_level_difference_at_vertices);
+ GridGenerator::quarter_hyper_ball(triangulation);
+ triangulation.refine_global(1);
+ std::set<types::boundary_id> no_flux_boundary{0, 1};
+ run<dim>(triangulation, no_flux_boundary);
+ }
+}
projects / dealii.git / commitdiff