namespace NonMatching
{
/**
- * This class 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$.
+ * This class defines a quadrature formula to integrate over the intersection
+ * between an oriented surface, $\hat{S}$, and a cell or face. The word
+ * "immersed" in the class name reflects that the surface may intersect the
+ * cell/face in an arbitrary way.
*
- * 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
+ * The spacedim template parameter of this class is the dimension that the
+ * (spacedim-1)-dimensional surface is embedded in:
+ * $\hat{S} \subset \mathbb{R}^{\text{spacedim}}$. The dim parameter
+ * describes the dimension of the "object" that the surface intersects. That
+ * is, dim = spacedim corresponds to the surface intersecting a cell and
+ * dim = spacedim - 1 corresponds to the surface intersecting a face. The
+ * quadrature formula is described by a set of quadrature points,
+ * $\hat{x}_q \in \mathbb{R}^{\text{dim}}$, weights, $w_q$, and normalized
+ * surface normals, $\hat{n}_q \in \mathbb{R}^{\text{spacedim}}$.
+ *
+ * Consider first the case dim = spacedim. We typically want to compute
+ * 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 an integral over $S\cap K$ in real space can be
+ * transformed to an 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
+ * \int_{S\cap K} f dS =
+ * \int_{S\cap K} f |d\bar{S}| =
+ * \int_{\hat{S}\cap\hat{K}} f \circ F_{K} \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,
+ * Jacobian matrix. This transformation is possible since the continuous
+ * surface elements are vectors:
+ * $d\bar{S}, d\hat{S} \in \mathbb{R}^{spacedim}$, which are parallel to the
+ * normals of $S$ and $\hat{S}$. That is, the normal is needed 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 \dealcoloneq 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
+ * \int_{S\cap K} f 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
+ *
+ * When dim = spacedim - 1, this class represents a (spacedim-2)-dimensional
+ * integral. That is, if spacedim = 3 we have a line integral immersed in a
+ * face. Here, the transformation between the face, $F$, and reference face,
+ * $\hat{F}$, reads
+ * @f[
+ * \int_{S\cap F} f dr =
+ * \int_{S\cap F} f |d\bar{r}| =
+ * \int_{\hat{S}\cap\hat{F}} f \circ F_{K} | J d\hat{r}|
+ * \approx \sum_{q} f \left(F_{K}(\hat{x}_{q}) \right) |J_q \hat{t}_q| w_q,
+ * @f]
+ * where $\hat{t}_q = \hat{n}_q \times \hat{n}_F$ is the tangent to the curve
+ * at $\hat{x}_q$ and $\hat{n}_F$ is the face normal. It would be possible to
+ * compute the tangent by only knowing the normal to the curve in the face
+ * plane (i.e. the dim-dimensional normal). However, when these quadratures
+ * are used, the weak form typically involves the so-called conormal, which
+ * can not be computed without knowing the surface normal in
+ * $\mathbb{R}^{\text{spacedim}}$. The conormal is the unit vector parallel to
+ * the projection of the face normal into the surface plane. This is
+ * essentially the same thing as the normalized @ref GlossBoundaryForm
+ * "boundary form".
*/
- template <int dim>
+ template <int dim, int spacedim = dim>
class ImmersedSurfaceQuadrature : public Quadrature<dim>
{
public:
* surface normals. The points, weights and normals should be with respect
* to reference space, and the normals should be normalized.
*/
- ImmersedSurfaceQuadrature(const std::vector<Point<dim>> & points,
- const std::vector<double> & weights,
- const std::vector<Tensor<1, dim>> &normals);
+ ImmersedSurfaceQuadrature(const std::vector<Point<dim>> & points,
+ const std::vector<double> & weights,
+ const std::vector<Tensor<1, spacedim>> &normals);
/**
* Extend the given formula by an additional quadrature point.
* quadrature formula.
*/
void
- push_back(const Point<dim> & point,
- const double weight,
- const Tensor<1, dim> &normal);
+ push_back(const Point<dim> & point,
+ const double weight,
+ const Tensor<1, spacedim> &normal);
/**
* Return a reference to the <tt>i</tt>th surface normal.
*/
- const Tensor<1, dim> &
+ const Tensor<1, spacedim> &
normal_vector(const unsigned int i) const;
/**
* Return a reference to the whole %vector of normals.
*/
- const std::vector<Tensor<1, dim>> &
+ const std::vector<Tensor<1, spacedim>> &
get_normal_vectors() const;
protected:
/**
* %Vector of surface normals at each quadrature point.
*/
- std::vector<Tensor<1, dim>> normals;
+ std::vector<Tensor<1, spacedim>> normals;
};
} // namespace NonMatching
DEAL_II_NAMESPACE_OPEN
namespace NonMatching
{
- template <int dim>
- ImmersedSurfaceQuadrature<dim>::ImmersedSurfaceQuadrature(
- const std::vector<Point<dim>> & points,
- const std::vector<double> & weights,
- const std::vector<Tensor<1, dim>> &normals)
+ template <int dim, int spacedim>
+ ImmersedSurfaceQuadrature<dim, spacedim>::ImmersedSurfaceQuadrature(
+ const std::vector<Point<dim>> & points,
+ const std::vector<double> & weights,
+ const std::vector<Tensor<1, spacedim>> &normals)
: Quadrature<dim>(points, weights)
, normals(normals)
{
- template <int dim>
+ template <int dim, int spacedim>
void
- ImmersedSurfaceQuadrature<dim>::push_back(const Point<dim> & point,
- const double weight,
- const Tensor<1, dim> &normal)
+ ImmersedSurfaceQuadrature<dim, spacedim>::push_back(
+ const Point<dim> & point,
+ const double weight,
+ const Tensor<1, spacedim> &normal)
{
this->quadrature_points.push_back(point);
this->weights.push_back(weight);
- template <int dim>
- const Tensor<1, dim> &
- ImmersedSurfaceQuadrature<dim>::normal_vector(const unsigned int i) const
+ template <int dim, int spacedim>
+ const Tensor<1, spacedim> &
+ ImmersedSurfaceQuadrature<dim, spacedim>::normal_vector(
+ const unsigned int i) const
{
AssertIndexRange(i, this->size());
return normals[i];
- template <int dim>
- const std::vector<Tensor<1, dim>> &
- ImmersedSurfaceQuadrature<dim>::get_normal_vectors() const
+ template <int dim, int spacedim>
+ const std::vector<Tensor<1, spacedim>> &
+ ImmersedSurfaceQuadrature<dim, spacedim>::get_normal_vectors() const
{
return normals;
}
- template class ImmersedSurfaceQuadrature<1>;
- template class ImmersedSurfaceQuadrature<2>;
- template class ImmersedSurfaceQuadrature<3>;
+ template class ImmersedSurfaceQuadrature<1, 1>;
+ template class ImmersedSurfaceQuadrature<2, 2>;
+ template class ImmersedSurfaceQuadrature<3, 3>;
+ template class ImmersedSurfaceQuadrature<1, 2>;
+ template class ImmersedSurfaceQuadrature<2, 3>;
} // namespace NonMatching
DEAL_II_NAMESPACE_CLOSE
-template <int dim>
+template <int dim, int spacedim>
void
-print_quadrature(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
+print_quadrature(
+ const NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> &quadrature)
{
for (unsigned int i = 0; i < quadrature.size(); ++i)
{
// Check that get_normals() are callable and are of the same size as
// points and weights.
-template <int dim>
+template <int dim, int spacedim>
void
-check_get_normals(const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature)
+check_get_normals(
+ const NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> &quadrature)
{
const std::vector<Point<dim>> & points = quadrature.get_points();
const std::vector<Tensor<1, dim>> &normals = quadrature.get_normal_vectors();
-template <int dim>
+template <int dim, int spacedim>
void
test_non_default_constructor()
{
deallog << "Using constructor" << std::endl;
- std::vector<Point<dim>> points(1);
- std::vector<double> weights(1, 1);
- std::vector<Tensor<1, dim>> normals;
- normals.push_back(Point<dim>::unit_vector(dim - 1));
- NonMatching::ImmersedSurfaceQuadrature<dim> quadrature(points,
- weights,
- normals);
+ std::vector<Point<dim>> points(1);
+ std::vector<double> weights(1, 1);
+ std::vector<Tensor<1, spacedim>> normals;
+ normals.push_back(Point<spacedim>::unit_vector(spacedim - 1));
+ NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> quadrature(points,
+ weights,
+ normals);
print_quadrature(quadrature);
}
-template <int dim>
+template <int dim, int spacedim>
void
test_push_back()
{
deallog << "Using push_back" << std::endl;
- const Point<dim> point;
- const double weight = 1;
- const Tensor<1, dim> normal = Point<dim>::unit_vector(dim - 1);
+ const Point<dim> point;
+ const double weight = 1;
+ const Tensor<1, spacedim> normal = Point<spacedim>::unit_vector(spacedim - 1);
- NonMatching::ImmersedSurfaceQuadrature<dim> quadrature;
+ NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> quadrature;
quadrature.push_back(point, weight, normal);
print_quadrature(quadrature);
-template <int dim>
+template <int dim, int spacedim>
void
construct_quadrature_and_print_points()
{
- test_push_back<dim>();
- test_non_default_constructor<dim>();
+ test_push_back<dim, spacedim>();
+ test_non_default_constructor<dim, spacedim>();
}
main()
{
initlog();
- construct_quadrature_and_print_points<1>();
- construct_quadrature_and_print_points<2>();
- construct_quadrature_and_print_points<3>();
+ construct_quadrature_and_print_points<1, 1>();
+ construct_quadrature_and_print_points<2, 2>();
+ construct_quadrature_and_print_points<3, 3>();
+
+ // Face quadrature
+ construct_quadrature_and_print_points<1, 2>();
+ construct_quadrature_and_print_points<2, 3>();
}
projects / dealii.git / commitdiff