+++ /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) 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