From 7a6a1a7c1e59b554e6cb324f47ee84aa357ab973 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Mon, 24 Apr 2017 19:49:04 +0200 Subject: [PATCH] Add class ImmersedSurfaceQuadrature Add a class representing an quadrature rule over a surface cutting arbitrary through the unit cell. In order to map a quadrature rule over a surface from the unit cell to physical space it is not sufficient to know only points and scalar weights, since information about the normal is also required. In addition to points and weights this class also stores the normal vectors. --- .../images/immersed_surface_quadrature.svg | 440 ++++++++++++++++++ .../immersed_surface_quadrature.h | 134 ++++++ source/CMakeLists.txt | 1 + source/non_matching/CMakeLists.txt | 30 ++ .../immersed_surface_quadrature.cc | 81 ++++ tests/non_matching/CMakeLists.txt | 4 + .../immersed_surface_quadrature.cc | 101 ++++ .../immersed_surface_quadrature.output | 13 + 8 files changed, 804 insertions(+) create mode 100644 doc/doxygen/images/immersed_surface_quadrature.svg create mode 100644 include/deal.II/non_matching/immersed_surface_quadrature.h create mode 100644 source/non_matching/CMakeLists.txt create mode 100644 source/non_matching/immersed_surface_quadrature.cc create mode 100644 tests/non_matching/CMakeLists.txt create mode 100644 tests/non_matching/immersed_surface_quadrature.cc create mode 100644 tests/non_matching/immersed_surface_quadrature.output diff --git a/doc/doxygen/images/immersed_surface_quadrature.svg b/doc/doxygen/images/immersed_surface_quadrature.svg new file mode 100644 index 0000000000..cbe33ab73f --- /dev/null +++ b/doc/doxygen/images/immersed_surface_quadrature.svg @@ -0,0 +1,440 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/deal.II/non_matching/immersed_surface_quadrature.h b/include/deal.II/non_matching/immersed_surface_quadrature.h new file mode 100644 index 0000000000..9d168bcd08 --- /dev/null +++ b/include/deal.II/non_matching/immersed_surface_quadrature.h @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__non_matching_immersed_surface_quadrature +#define dealii__non_matching_immersed_surface_quadrature + +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + + /** + * Defines a quadrature formula for integration over an oriented surface, + * $\hat{S}$, immersed in the unit cell. By immersed it is meant that the + * surface may intersect the unit cell in an arbitrary way. The quadrature + * formula is described by a set of quadrature points, $\hat{x}_q$, weights, + * $w_q$, and normalized surface normals, $\hat{n}_q$. + * + * We typically want to compute surface integrals in real space. + * A surface $S$ intersecting a cell $K$ in real space, can be mapped onto a + * surface $\hat{S}$ intersecting the unit cell $\hat{K}$. + * Thus a surface integral over $S\cap K$ in real space can be transformed to + * a surface integral over $\hat{S} \cap \hat{K}$ according to + * @f[ + * \int_{S\cap K} f(x) dS = + * \int_{S\cap K} f(x) |d\bar{S}| = + * \int_{\hat{S}\cap\hat{K}} f(F_{K}(\hat{x})) \det(J) |\left( J^{-1} \right )^T d\hat{S}|, + * @f] + * where $F_K$ is the mapping from reference to real space and $J$ is its + * Jacobian. This transformation is possible since the continuous surface + * elements are vectors: $d\bar{S}, d\hat{S} \in \mathbb{R}^{dim}$ which are + * parallel to the normals of $S$ and $\hat{S}$. So in order to compute the + * integral in real space one needs information about the normal to do the + * transformation. + * + * Thus, in addition to storing points and weights, this quadrature stores + * also the normalized normal for each quadrature point. This can be viewed + * as storing a discrete surface element, + * @f[ + * \Delta \hat{S}_q := w_q \hat{n}_q \approx d\hat{S}(\hat{x}_q), + * @f] + * for each quadrature point. The surface integral in real space would then be + * approximated as + * @f[ + * \int_{S\cap K} f(x) dS \approx + * \sum_{q} f \left(F_{K}(\hat{x}_{q}) \right) \det(J_q) + * |\left( J_q^{-1} \right)^T \hat{n}_q| w_q. + * @f] + * + * @image html immersed_surface_quadrature.svg + * + * @author Simon Sticko, 2017 + */ + template + class ImmersedSurfaceQuadrature : public Quadrature + { + public: + + /** + * Default constructor to initialize the quadrature with no quadrature + * points. + */ + ImmersedSurfaceQuadrature()=default; + + /** + * Construct a quadrature formula from vectors of points, weights and + * surface normals. The points, weights and normals should be with respect + * to reference space, and the normals should be normalized. + */ + ImmersedSurfaceQuadrature( + const std::vector> &points, + const std::vector &weights, + const std::vector> &normals); + + /** + * Extend the given formula by an additional quadrature point. + * The point, weight and normal should be with respect to reference space, + * and the normal should be normalized. + * + * This function exists since immersed quadrature rules can be rather + * complicated to construct. Often the construction is done by + * partitioning the cell into regions and constructing points on each + * region separately. This can make it cumbersome to create the quadrature + * from the constructor since all quadrature points have to be known at + * time of creation of the object. + * + * @note This function should only be used during construction of the + * quadrature formula. + */ + void push_back(const Point &point, + const double weight, + const Tensor<1, dim> &normal); + + /** + * Return a reference to the ith surface normal. + */ + const Tensor<1,dim> &normal_vector(const unsigned int i) const; + + /** + * Return a reference to the whole %vector of normals. + */ + const std::vector> &get_normal_vectors() const; + + protected: + + /** + * %Vector of surface normals at each quadrature point. + */ + std::vector> normals; + }; + +} +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index b78b0a0f1b..36e209f704 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -49,6 +49,7 @@ ADD_SUBDIRECTORY(matrix_free) ADD_SUBDIRECTORY(meshworker) ADD_SUBDIRECTORY(opencascade) ADD_SUBDIRECTORY(physics) +ADD_SUBDIRECTORY(non_matching) FOREACH(build ${DEAL_II_BUILD_TYPES}) STRING(TOLOWER ${build} build_lowercase) diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt new file mode 100644 index 0000000000..085f7f8ccb --- /dev/null +++ b/source/non_matching/CMakeLists.txt @@ -0,0 +1,30 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2012 - 2014 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 at +## the top level of the deal.II distribution. +## +## --------------------------------------------------------------------- + +INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_src + immersed_surface_quadrature.cc + ) + +SET(_inst + ) + +FILE(GLOB _header + ${CMAKE_SOURCE_DIR}/include/deal.II/non_matching/*.h + ) + +DEAL_II_ADD_LIBRARY(obj_non_matching OBJECT ${_src} ${_header} ${_inst}) +EXPAND_INSTANTIATIONS(obj_non_matching "${_inst}") diff --git a/source/non_matching/immersed_surface_quadrature.cc b/source/non_matching/immersed_surface_quadrature.cc new file mode 100644 index 0000000000..288620d14a --- /dev/null +++ b/source/non_matching/immersed_surface_quadrature.cc @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + + template + ImmersedSurfaceQuadrature::ImmersedSurfaceQuadrature( + const std::vector > &points, + const std::vector &weights, + const std::vector> &normals): + Quadrature(points,weights), + normals(normals) + { + AssertDimension (weights.size(), points.size()); + AssertDimension (normals.size(), points.size()); + for (auto normal : normals) + { + Assert(std::abs(normal.norm() - 1.0) < 1e-9, + ExcMessage("Normal is not normalized.")); + } + } + + + + template + void ImmersedSurfaceQuadrature::push_back( + const Point &point, + const double weight, + const Tensor<1, dim> &normal) + { + this->quadrature_points.push_back(point); + this->weights.push_back(weight); + this->normals.push_back(normal); + Assert(std::abs(normal.norm() - 1.0) < 1e-9, + ExcMessage("Normal is not normalized.")); + } + + + + template + const Tensor<1,dim> & + ImmersedSurfaceQuadrature::normal_vector( + const unsigned int i) const + { + AssertIndexRange(i, this->size()); + return normals[i]; + } + + + + template + const std::vector > & + ImmersedSurfaceQuadrature::get_normal_vectors() const + { + return normals; + } + + + + template class ImmersedSurfaceQuadrature<1>; + template class ImmersedSurfaceQuadrature<2>; + template class ImmersedSurfaceQuadrature<3>; + +} +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/non_matching/CMakeLists.txt b/tests/non_matching/CMakeLists.txt new file mode 100644 index 0000000000..f8ac3eb029 --- /dev/null +++ b/tests/non_matching/CMakeLists.txt @@ -0,0 +1,4 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.9) +INCLUDE(../setup_testsubproject.cmake) +PROJECT(testsuite CXX) +DEAL_II_PICKUP_TESTS() diff --git a/tests/non_matching/immersed_surface_quadrature.cc b/tests/non_matching/immersed_surface_quadrature.cc new file mode 100644 index 0000000000..51307b8811 --- /dev/null +++ b/tests/non_matching/immersed_surface_quadrature.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include "../tests.h" + +using namespace dealii; + +//Test that an ImmersedSurfaceQuadrature can be constructed for each dimension +//and that quadrature points can be added to it. + + + +template +void print_quadrature(const NonMatching::ImmersedSurfaceQuadrature &quadrature) +{ + for (unsigned int i = 0; i < quadrature.size(); ++i) + { + deallog << quadrature.point(i) << ", " << quadrature.weight(i) << ", " + << quadrature.normal_vector(i) << std::endl; + } +} + + + +//Check that get_normals() are callable and are of the same size as +//points and weights. +template +void check_get_normals( + const NonMatching::ImmersedSurfaceQuadrature &quadrature) +{ + const std::vector> &points= + quadrature.get_points(); + const std::vector> &normals= + quadrature.get_normal_vectors(); + AssertThrow(points.size()==normals.size(),ExcInternalError()) +} + + + +template +void test_non_default_constructor() +{ + deallog<<"Using constructor"<> points(1); + std::vector weights(1,1); + std::vector> normals; + normals.push_back(Point::unit_vector(dim-1)); + NonMatching::ImmersedSurfaceQuadrature quadrature(points,weights,normals); + + print_quadrature(quadrature); +} + + + +template +void test_push_back() +{ + deallog<<"Using push_back"< point; + const double weight=1; + const Tensor<1,dim> normal=Point::unit_vector(dim-1); + + NonMatching::ImmersedSurfaceQuadrature quadrature; + quadrature.push_back(point,weight,normal); + + print_quadrature(quadrature); +} + + + +template +void construct_quadrature_and_print_points() +{ + test_push_back(); + test_non_default_constructor(); +} + + + +int main() +{ + initlog(); + construct_quadrature_and_print_points<1>(); + construct_quadrature_and_print_points<2>(); + construct_quadrature_and_print_points<3>(); +} diff --git a/tests/non_matching/immersed_surface_quadrature.output b/tests/non_matching/immersed_surface_quadrature.output new file mode 100644 index 0000000000..7143d449ee --- /dev/null +++ b/tests/non_matching/immersed_surface_quadrature.output @@ -0,0 +1,13 @@ + +DEAL::Using push_back +DEAL::0.00000, 1.00000, 1.00000 +DEAL::Using constructor +DEAL::0.00000, 1.00000, 1.00000 +DEAL::Using push_back +DEAL::0.00000 0.00000, 1.00000, 0.00000 1.00000 +DEAL::Using constructor +DEAL::0.00000 0.00000, 1.00000, 0.00000 1.00000 +DEAL::Using push_back +DEAL::0.00000 0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 +DEAL::Using constructor +DEAL::0.00000 0.00000 0.00000, 1.00000, 0.00000 0.00000 1.00000 -- 2.39.5