--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_matching_quadrature_generator
+#define dealii_non_matching_quadrature_generator
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <type_traits>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace NonMatching
+{
+ namespace internal
+ {
+ namespace QuadratureGeneratorImplementation
+ {
+ /**
+ * Returns the center $(0.5, ..., 0.5)$ of the unit hypercube:
+ * $[0,1]^{dim}$.
+ */
+ template <int dim>
+ Point<dim>
+ unit_hypercube_center();
+
+
+ /**
+ * Describes a dim-dimensional hypercube aligned with the coordinate
+ * axes. That is, a region
+ *
+ * @f[
+ * [x_0-L/2, x_0+L/2] \times ... \times [x_{dim-1}-L/2, x_{dim-1}+L/2]
+ * @f]
+ *
+ * where L is the side length and $(x_0, ..., x_{dim-1})$ are the
+ * coordinates of the center of the hypercube.
+ *
+ * The purpose of this class is to help us keep track of where we are on
+ * the reference cell when generating immersed quadrature rules. There
+ * are two main things we want to be able to do.
+ *
+ * 1. Split a hypercube into its $2^{dim}$ children.
+ * 2. Take a cross section orthogonal to a given direction and get this
+ * as a Hypercube<dim-1>.
+ *
+ * These are needed because the algorithm is allowed to recurse both
+ * over children and over dimensions.
+ *
+ * In 3D, the 2 coordinates of the cross section of Hypercube<3> can be
+ * ordered in 2 ways. That is, if we take the cross section orthogonal
+ * to the y direction we could either order a 3D-coordinate into a
+ * 2D-coordinate as $(x,z)$ or as $(z,x)$. This class uses the second
+ * convention, corresponding to the coordinates being ordered cyclicly
+ * $x \rightarrow y \rightarrow z \rightarrow x \rightarrow ... $
+ * To be precise, if we take a cross section:
+ *
+ * | Orthogonal to | Cross section coordinates ordered as |
+ * |:-------------:|:------------------------------------:|
+ * | x | (y, z) |
+ * | y | (z, x) |
+ * | z | (x, y) |
+ *
+ * This is according to the convention set by the function
+ * <code>coordinate_to_one_dim_higher</code>.
+ */
+ template <int dim>
+ class Hypercube
+ {
+ public:
+ /**
+ * Constructor, takes the center of the hypercube and its side length.
+ * Passing no inarguments creates a hypercube between 0 and 1 in all
+ * directions: $[0, 1]^{dim}$.
+ */
+ Hypercube(
+ const Point<dim> ¢er_coordinate = unit_hypercube_center<dim>(),
+ const double length_of_side = 1);
+
+ /**
+ * Return the lower bound of the hypercube in @p direction.
+ */
+ double
+ lower_bound(const unsigned int direction) const;
+
+ /**
+ * Return the upper bound of the hypercube in @p direction.
+ */
+ double
+ upper_bound(const unsigned int direction) const;
+
+ /**
+ * Return the bounds of the hypercube in @p direction, as a one-dimensional
+ * hypercube.
+ */
+ Hypercube<1>
+ bounds(const unsigned int direction) const;
+
+ /**
+ * Returns the cross section of the hypercube orthogonal to @p direction.
+ * This is a hypercube in one dimension lower.
+ */
+ Hypercube<dim - 1>
+ cross_section(const unsigned int direction) const;
+
+ /**
+ * Returns the point in the center of the hypercube.
+ */
+ const Point<dim> &
+ center() const;
+
+ /**
+ * Returns the length of the side of the hypercube.
+ */
+ double
+ side_length() const;
+
+ /**
+ * Returns the volume of the hypercube.
+ */
+ double
+ volume() const;
+
+ /**
+ * Returns the indexth vertex of the hypercube.
+ */
+ Point<dim>
+ vertex(const unsigned int index) const;
+
+ /**
+ * Returns the indexth child of the hypercube. Child is meant in the
+ * same way as for a cell.
+ */
+ Hypercube<dim>
+ child(const unsigned int index) const;
+
+ private:
+ /**
+ * Coordinate in the center of the hypercube.
+ */
+ const Point<dim> center_coordinate;
+
+ /**
+ * Side length of the hypercube.
+ */
+ const double length_of_side;
+ };
+
+
+ /**
+ * This function defines a convention for how coordinates in dim
+ * dimensions should translate to the coordinates in dim + 1 dimensions,
+ * when one of the coordinates in dim + 1 dimensions is locked to a given
+ * value.
+ *
+ * The convention is the following: Starting from the locked coordinate we
+ * store the lower dimensional coordinates consecutively and wrapping
+ * around when going over the dimension. This relationship can be
+ * described by the following tables:
+ *
+ * 2D
+ * ------------------------------------------------
+ * | locked in 2D | 1D coordinate | 2D coordinate |
+ * |--------------------------------------------- |
+ * | x0 | (a) | (x0, a) |
+ * | x1 | (a) | (a , x1) |
+ * ------------------------------------------------
+ *
+ * 3D
+ * --------------------------------------------------
+ * | locked in 3D | 2D coordinates | 3D coordinates |
+ * |----------------------------------------------- |
+ * | x0 | (a, b) | (x0, a, b) |
+ * | x1 | (a, b) | ( b, x1, a) |
+ * | x2 | (a, b) | ( a, b, x2) |
+ * --------------------------------------------------
+ *
+ * Given a locked coordinate, this function maps a coordinate index in dim
+ * dimension to a coordinate index in dim + 1 dimensions.
+ *
+ * @param locked_coordinate should be in the range [0, dim+1).
+ * @param coordiante_in_dim should be in the range [0, dim).
+ * @return A coordinate index in the range [0, dim+1)
+ */
+ template <int dim>
+ inline unsigned int
+ coordinate_to_one_dim_higher(const unsigned int locked_coordinate,
+ const unsigned int coordiante_in_dim)
+ {
+ AssertIndexRange(locked_coordinate, dim + 1);
+ AssertIndexRange(coordiante_in_dim, dim);
+ return (locked_coordinate + coordiante_in_dim + 1) % (dim + 1);
+ }
+
+
+
+ /* -------------- declaration of explicit specializations ------------- */
+
+ template <>
+ Hypercube<-1>
+ Hypercube<0>::cross_section(const unsigned int) const;
+
+ template <>
+ Point<0>
+ Hypercube<0>::vertex(const unsigned int) const;
+
+ template <>
+ Hypercube<0>
+ Hypercube<0>::child(const unsigned int) const;
+
+ } // namespace QuadratureGeneratorImplementation
+ } // namespace internal
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
+
+#endif /* dealii_non_matching_quadrature_generator */
SET(_src
coupling.cc
immersed_surface_quadrature.cc
+ quadrature_generator.cc
)
SET(_inst
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/geometry_info.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace NonMatching
+{
+ namespace internal
+ {
+ namespace QuadratureGeneratorImplementation
+ {
+ template <int dim>
+ Point<dim>
+ unit_hypercube_center()
+ {
+ Point<dim> center;
+ for (unsigned int i = 0; i < dim; ++i)
+ center[i] = .5;
+
+ return center;
+ }
+
+
+
+ template <int dim>
+ Hypercube<dim>::Hypercube(const Point<dim> ¢er_coordinate,
+ const double length_of_side)
+ : center_coordinate(center_coordinate)
+ , length_of_side(length_of_side)
+ {
+ Assert(0 < dim, ExcImpossibleInDim(dim));
+ }
+
+
+
+ template <int dim>
+ double
+ Hypercube<dim>::lower_bound(const unsigned int direction) const
+ {
+ AssertIndexRange(direction, dim);
+
+ return center_coordinate[direction] - length_of_side / 2;
+ }
+
+
+
+ template <int dim>
+ double
+ Hypercube<dim>::upper_bound(const unsigned int direction) const
+ {
+ AssertIndexRange(direction, dim);
+
+ return center_coordinate[direction] + length_of_side / 2;
+ }
+
+
+
+ template <int dim>
+ Hypercube<1>
+ Hypercube<dim>::bounds(const unsigned int direction) const
+ {
+ AssertIndexRange(direction, dim);
+
+ const Point<1> center(center_coordinate[direction]);
+ return Hypercube<1>(center, side_length());
+ }
+
+
+
+ template <int dim>
+ const Point<dim> &
+ Hypercube<dim>::center() const
+ {
+ return center_coordinate;
+ }
+
+
+
+ template <int dim>
+ double
+ Hypercube<dim>::side_length() const
+ {
+ return length_of_side;
+ }
+
+
+
+ template <int dim>
+ double
+ Hypercube<dim>::volume() const
+ {
+ return std::pow(length_of_side, dim);
+ }
+
+
+
+ template <int dim>
+ Hypercube<dim - 1>
+ Hypercube<dim>::cross_section(const unsigned int direction) const
+ {
+ AssertIndexRange(direction, dim);
+
+ Point<dim - 1> center_of_cross_section;
+ for (unsigned int d = 0; d < dim - 1; ++d)
+ {
+ const int index_to_write_from =
+ coordinate_to_one_dim_higher<dim - 1>(direction, d);
+ center_of_cross_section[d] = center_coordinate[index_to_write_from];
+ }
+
+ return Hypercube<dim - 1>(center_of_cross_section, length_of_side);
+ }
+
+
+
+ template <int dim>
+ Point<dim>
+ Hypercube<dim>::vertex(const unsigned int index) const
+ {
+ AssertIndexRange(index, GeometryInfo<dim>::vertices_per_cell);
+
+ Point<dim> bottom_corner;
+ for (unsigned int i = 0; i < dim; ++i)
+ bottom_corner[i] = center_coordinate[i] - length_of_side / 2;
+
+ const Point<dim> point =
+ bottom_corner +
+ length_of_side * GeometryInfo<dim>::unit_cell_vertex(index);
+
+ return point;
+ }
+
+
+
+ template <int dim>
+ Hypercube<dim>
+ Hypercube<dim>::child(const unsigned int index) const
+ {
+ AssertIndexRange(index, GeometryInfo<dim>::max_children_per_cell);
+
+ const double child_side_length = length_of_side / 2;
+
+ const Point<dim> unit_cell_vertex =
+ GeometryInfo<dim>::unit_cell_vertex(index);
+
+ Point<dim> child_center;
+ for (unsigned int i = 0; i < dim; ++i)
+ child_center[i] = center_coordinate[i] +
+ child_side_length * (unit_cell_vertex[i] - .5);
+
+ return Hypercube<dim>(child_center, child_side_length);
+ }
+
+
+
+ template Point<0>
+ unit_hypercube_center<0>();
+ template Point<1>
+ unit_hypercube_center<1>();
+ template Point<2>
+ unit_hypercube_center<2>();
+ template Point<3>
+ unit_hypercube_center<3>();
+
+ template class Hypercube<0>;
+ template class Hypercube<1>;
+ template class Hypercube<2>;
+ template class Hypercube<3>;
+
+ } // namespace QuadratureGeneratorImplementation
+ } // namespace internal
+
+
+
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+/*
+ * Lock all possible coordinates in dim + 1. For each locked coordinate, print
+ * how the coordinates in dim dimensions map to the coordinates in dim + 1
+ * dimensions.
+ */
+template <int dim>
+void
+print_what_each_coordinate_maps_to()
+{
+ deallog << "dim = " << dim << std::endl;
+
+ for (unsigned int locked = 0; locked < dim + 1; ++locked)
+ {
+ deallog << "locked coordinate = " << locked << std::endl;
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ deallog << i << " -> "
+ << NonMatching::internal::QuadratureGeneratorImplementation::
+ coordinate_to_one_dim_higher<dim>(locked, i)
+ << std::endl;
+ }
+ }
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ print_what_each_coordinate_maps_to<1>();
+ deallog << std::endl;
+ print_what_each_coordinate_maps_to<2>();
+
+ return 0;
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::locked coordinate = 0
+DEAL::0 -> 1
+DEAL::locked coordinate = 1
+DEAL::0 -> 0
+DEAL::
+DEAL::dim = 2
+DEAL::locked coordinate = 0
+DEAL::0 -> 1
+DEAL::1 -> 2
+DEAL::locked coordinate = 1
+DEAL::0 -> 2
+DEAL::1 -> 0
+DEAL::locked coordinate = 2
+DEAL::0 -> 0
+DEAL::1 -> 1
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 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.
+//
+// ---------------------------------------------------------------------
+
+// Test all member functions of the HyperCube<dim> class. Each function is
+// tested in the function named test_{member_function_name}.
+
+#include <deal.II/base/geometry_info.h>
+
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using NonMatching::internal::QuadratureGeneratorImplementation::Hypercube;
+
+// Print the center and the side length of the incoming hypercube, which are
+// its degrees of freedom.
+template <int dim>
+void
+print_cube(const Hypercube<dim> &hyper_cube)
+{
+ deallog << "center = " << hyper_cube.center() << ", ";
+ deallog << "side = " << hyper_cube.side_length() << std::endl;
+}
+
+
+
+// Construct the cross section orthogonal to all axes and print them.
+template <int dim>
+void
+test_cross_section()
+{
+ deallog << "test_cross_section" << std::endl;
+
+ const double width = 1;
+ Point<dim> center;
+ for (unsigned int d = 0; d < dim; ++d)
+ center[d] = d;
+
+ const Hypercube<dim> hypercube(center, width);
+
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ deallog << "orthogonal to " << d << std::endl;
+ const Hypercube<dim - 1> down = hypercube.cross_section(d);
+ print_cube(down);
+ }
+}
+
+
+
+// 1D-specialization of the above test. We can't take a cross section in 1D, so
+// we do nothing.
+template <>
+void
+test_cross_section<1>()
+{}
+
+
+
+// Construct all possible children of a hypercube and print them.
+template <int dim>
+void
+test_child()
+{
+ deallog << "test_child" << std::endl;
+
+ const double width = 4;
+ Point<dim> center;
+
+ const Hypercube<dim> hypercube(center, width);
+
+ // Take the cross section in all directions.
+ for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell; ++i)
+ {
+ deallog << "child " << i << std::endl;
+ const Hypercube<dim> child = hypercube.child(i);
+ print_cube(child);
+ }
+}
+
+
+
+// Comptue and print the volume of a hypercube.
+template <int dim>
+void
+test_volume()
+{
+ deallog << "test_volume" << std::endl;
+
+ const double width = 2;
+ Point<dim> center;
+
+ const Hypercube<dim> hypercube(center, width);
+ deallog << "volume = " << hypercube.volume() << std::endl;
+}
+
+
+
+// Print the result of the lower_bound and upper_bound functions in all
+// coordinate directions.
+template <int dim>
+void
+test_lower_upper_bound()
+{
+ deallog << "test_lower_upper_bound" << std::endl;
+
+ const double width = 2;
+ Point<dim> center;
+ for (unsigned int d = 0; d < dim; ++d)
+ center[d] = d + 1;
+
+ const Hypercube<dim> hypercube(center, width);
+ for (int d = 0; d < dim; ++d)
+ {
+ deallog << "[" << hypercube.lower_bound(d) << ", "
+ << hypercube.upper_bound(d) << "]";
+ if (dim - 1 != d)
+ {
+ deallog << "x";
+ }
+ }
+}
+
+
+
+// Extract and print all the vertices of a hypercube.
+template <int dim>
+void
+test_vertex()
+{
+ deallog << "test_vertex" << std::endl;
+
+ const Hypercube<dim> hypercube;
+ for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ deallog << "vertex " << i << " = " << hypercube.vertex(i) << std::endl;
+ }
+}
+
+
+
+// Print the results of the bounds function for all coordinate directions.
+template <int dim>
+void
+test_bounds()
+{
+ deallog << "test_bounds" << std::endl;
+
+ const double width = 2;
+ Point<dim> center;
+ for (unsigned int d = 0; d < dim; ++d)
+ center[d] = d + 1;
+
+ const Hypercube<dim> hypercube(center, width);
+
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ const Hypercube<1> bounds = hypercube.bounds(i);
+ deallog << "Bounds in direction " << i << std::endl;
+ print_cube(bounds);
+ }
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ deallog << "dim = " << dim << std::endl;
+
+ test_cross_section<dim>();
+ deallog << std::endl;
+
+ test_child<dim>();
+ deallog << std::endl;
+
+ test_volume<dim>();
+ deallog << std::endl;
+
+ test_lower_upper_bound<dim>();
+ deallog << std::endl;
+
+ test_vertex<dim>();
+ deallog << std::endl;
+
+ test_bounds<dim>();
+ deallog << std::endl;
+
+ deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::
+DEAL::test_child
+DEAL::child 0
+DEAL::center = -1.00000, side = 2.00000
+DEAL::child 1
+DEAL::center = 1.00000, side = 2.00000
+DEAL::
+DEAL::test_volume
+DEAL::volume = 2.00000
+DEAL::
+DEAL::test_lower_upper_bound
+DEAL::[0.00000, 2.00000]
+DEAL::test_vertex
+DEAL::vertex 0 = 0.00000
+DEAL::vertex 1 = 1.00000
+DEAL::
+DEAL::test_bounds
+DEAL::Bounds in direction 0
+DEAL::center = 1.00000, side = 2.00000
+DEAL::
+DEAL::
+DEAL::dim = 2
+DEAL::test_cross_section
+DEAL::orthogonal to 0
+DEAL::center = 1.00000, side = 1.00000
+DEAL::orthogonal to 1
+DEAL::center = 0.00000, side = 1.00000
+DEAL::
+DEAL::test_child
+DEAL::child 0
+DEAL::center = -1.00000 -1.00000, side = 2.00000
+DEAL::child 1
+DEAL::center = 1.00000 -1.00000, side = 2.00000
+DEAL::child 2
+DEAL::center = -1.00000 1.00000, side = 2.00000
+DEAL::child 3
+DEAL::center = 1.00000 1.00000, side = 2.00000
+DEAL::
+DEAL::test_volume
+DEAL::volume = 4.00000
+DEAL::
+DEAL::test_lower_upper_bound
+DEAL::[0.00000, 2.00000]x[1.00000, 3.00000]
+DEAL::test_vertex
+DEAL::vertex 0 = 0.00000 0.00000
+DEAL::vertex 1 = 1.00000 0.00000
+DEAL::vertex 2 = 0.00000 1.00000
+DEAL::vertex 3 = 1.00000 1.00000
+DEAL::
+DEAL::test_bounds
+DEAL::Bounds in direction 0
+DEAL::center = 1.00000, side = 2.00000
+DEAL::Bounds in direction 1
+DEAL::center = 2.00000, side = 2.00000
+DEAL::
+DEAL::
+DEAL::dim = 3
+DEAL::test_cross_section
+DEAL::orthogonal to 0
+DEAL::center = 1.00000 2.00000, side = 1.00000
+DEAL::orthogonal to 1
+DEAL::center = 2.00000 0.00000, side = 1.00000
+DEAL::orthogonal to 2
+DEAL::center = 0.00000 1.00000, side = 1.00000
+DEAL::
+DEAL::test_child
+DEAL::child 0
+DEAL::center = -1.00000 -1.00000 -1.00000, side = 2.00000
+DEAL::child 1
+DEAL::center = 1.00000 -1.00000 -1.00000, side = 2.00000
+DEAL::child 2
+DEAL::center = -1.00000 1.00000 -1.00000, side = 2.00000
+DEAL::child 3
+DEAL::center = 1.00000 1.00000 -1.00000, side = 2.00000
+DEAL::child 4
+DEAL::center = -1.00000 -1.00000 1.00000, side = 2.00000
+DEAL::child 5
+DEAL::center = 1.00000 -1.00000 1.00000, side = 2.00000
+DEAL::child 6
+DEAL::center = -1.00000 1.00000 1.00000, side = 2.00000
+DEAL::child 7
+DEAL::center = 1.00000 1.00000 1.00000, side = 2.00000
+DEAL::
+DEAL::test_volume
+DEAL::volume = 8.00000
+DEAL::
+DEAL::test_lower_upper_bound
+DEAL::[0.00000, 2.00000]x[1.00000, 3.00000]x[2.00000, 4.00000]
+DEAL::test_vertex
+DEAL::vertex 0 = 0.00000 0.00000 0.00000
+DEAL::vertex 1 = 1.00000 0.00000 0.00000
+DEAL::vertex 2 = 0.00000 1.00000 0.00000
+DEAL::vertex 3 = 1.00000 1.00000 0.00000
+DEAL::vertex 4 = 0.00000 0.00000 1.00000
+DEAL::vertex 5 = 1.00000 0.00000 1.00000
+DEAL::vertex 6 = 0.00000 1.00000 1.00000
+DEAL::vertex 7 = 1.00000 1.00000 1.00000
+DEAL::
+DEAL::test_bounds
+DEAL::Bounds in direction 0
+DEAL::center = 1.00000, side = 2.00000
+DEAL::Bounds in direction 1
+DEAL::center = 2.00000, side = 2.00000
+DEAL::Bounds in direction 2
+DEAL::center = 3.00000, side = 2.00000
+DEAL::
+DEAL::
projects / dealii.git / commitdiff