From: Timo Heister <timo.heister@gmail.com>
Date: Wed, 2 Oct 2019 17:20:15 +0000 (-0400)
Subject: FEInterfaceValues: gradients and hessians
X-Git-Tag: v9.2.0-rc1~943^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bb0735cd65243d2f41fffcd556a3a311a3c6ed8c;p=dealii.git

FEInterfaceValues: gradients and hessians

add average_gradient(), average_hessian(), jump_gradient()
---

diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h
index 9e611039cc..50100a27c8 100644
--- a/include/deal.II/fe/fe_interface_values.h
+++ b/include/deal.II/fe/fe_interface_values.h
@@ -307,6 +307,47 @@ public:
           const unsigned int q_point,
           const unsigned int component = 0) const;
 
+  /**
+   * Return the average of the gradient $\{\nabla u \} = \frac{1}{2}\nabla
+   * u_{\text{cell0}} + \frac{1}{2} \nabla u_{\text{cell1}}$ on the interface
+   * for the shape function @p interface_dof_index at the quadrature point @p
+   * q_point of component @p component.
+   *
+   * If this is a boundary face (at_boundary() returns true), then
+   * $\{\nabla u\}=\nabla u_{\text{cell0}}$.
+   */
+  Tensor<1, dim>
+  average_gradient(const unsigned int interface_dof_index,
+                   const unsigned int q_point,
+                   const unsigned int component = 0) const;
+
+  /**
+   * Return the average of the hessian $\{\nabla^2 u \} = \frac{1}{2}\nabla^2
+   * u_{\text{cell0}} + \frac{1}{2} \nabla^2 u_{\text{cell1}}$ on the interface
+   * for the shape function @p interface_dof_index at the quadrature point @p
+   * q_point of component @p component.
+   *
+   * If this is a boundary face (at_boundary() returns true), then
+   * $\{\nabla^2 u\}=\nabla^2 u_{\text{cell0}}$.
+   */
+  Tensor<2, dim>
+  average_hessian(const unsigned int interface_dof_index,
+                  const unsigned int q_point,
+                  const unsigned int component = 0) const;
+
+  /**
+   * Return the jump if the gradient $[\nabla u]=\nabla u_{\text{cell0}} -
+   * \nabla u_{\text{cell1}}$ on the interface for the shape function @p
+   * interface_dof_index at the quadrature point @p q_point of component @p
+   * component.
+   *
+   * If this is a boundary face (at_boundary() returns true), then
+   * $[\nabla u]=\nabla u_{\text{cell0}}$.
+   */
+  Tensor<1, dim>
+  jump_gradient(const unsigned int interface_dof_index,
+                const unsigned int q_point,
+                const unsigned int component = 0) const;
 
   /**
    * @}
@@ -705,6 +746,95 @@ FEInterfaceValues<dim, spacedim>::average(
   return value;
 }
 
+
+
+template <int dim, int spacedim>
+Tensor<1, dim>
+FEInterfaceValues<dim, spacedim>::average_gradient(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return get_fe_face_values(0).shape_grad_component(dofpair[0],
+                                                      q_point,
+                                                      component);
+
+  Tensor<1, dim> value;
+
+  if (dofpair[0] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(0).shape_grad_component(dofpair[0],
+                                                              q_point,
+                                                              component);
+  if (dofpair[1] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(1).shape_grad_component(dofpair[1],
+                                                              q_point,
+                                                              component);
+
+  return value;
+}
+
+
+
+template <int dim, int spacedim>
+Tensor<2, dim>
+FEInterfaceValues<dim, spacedim>::average_hessian(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return get_fe_face_values(0).shape_hessian_component(dofpair[0],
+                                                         q_point,
+                                                         component);
+
+  Tensor<2, dim> value;
+
+  if (dofpair[0] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(0).shape_hessian_component(dofpair[0],
+                                                                 q_point,
+                                                                 component);
+  if (dofpair[1] != numbers::invalid_unsigned_int)
+    value += 0.5 * get_fe_face_values(1).shape_hessian_component(dofpair[1],
+                                                                 q_point,
+                                                                 component);
+
+  return value;
+}
+
+template <int dim, int spacedim>
+Tensor<1, dim>
+FEInterfaceValues<dim, spacedim>::jump_gradient(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dofpair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return get_fe_face_values(0).shape_grad_component(dofpair[0],
+                                                      q_point,
+                                                      component);
+
+  Tensor<1, dim> value;
+
+  if (dofpair[0] != numbers::invalid_unsigned_int)
+    value += 1.0 * get_fe_face_values(0).shape_grad_component(dofpair[0],
+                                                              q_point,
+                                                              component);
+  if (dofpair[1] != numbers::invalid_unsigned_int)
+    value += -1.0 * get_fe_face_values(1).shape_grad_component(dofpair[1],
+                                                               q_point,
+                                                               component);
+
+  return value;
+}
+
+
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/feinterface/fe_interface_values_08.cc b/tests/feinterface/fe_interface_values_08.cc
new file mode 100644
index 0000000000..8568b49c23
--- /dev/null
+++ b/tests/feinterface/fe_interface_values_08.cc
@@ -0,0 +1,170 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate average_gradient average_hessian jump_gradient of FEInterfaceValues
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(const FiniteElement<dim> &fe)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim> dofh(tria);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_hessians | update_quadrature_points |
+                             update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+
+  auto cell = dofh.begin();
+
+  {
+    unsigned int f = 0;
+    Assert(cell->at_boundary(f), ExcInternalError());
+
+    fiv.reinit(cell, f);
+
+    const unsigned int n_dofs = fiv.n_current_interface_dofs();
+    Vector<double>     cell_vector(n_dofs);
+
+    const auto &q_points = fiv.get_quadrature_points();
+
+    cell_vector = 0.0;
+    for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        cell_vector(i) +=
+          fiv.average_gradient(i, qpoint).norm() * fiv.get_JxW_values()[qpoint];
+    deallog << "average_gradient.norm(): " << cell_vector;
+
+    cell_vector = 0.0;
+    for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        cell_vector(i) +=
+          fiv.average_hessian(i, qpoint).norm() * fiv.get_JxW_values()[qpoint];
+    deallog << "average_hessian.norm(): " << cell_vector;
+
+    cell_vector = 0.0;
+    for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+      for (unsigned int i = 0; i < n_dofs; ++i)
+        cell_vector(i) +=
+          fiv.jump_gradient(i, qpoint).norm() * fiv.get_JxW_values()[qpoint];
+    deallog << "jump_gradient.norm(): " << cell_vector;
+  }
+
+
+  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+    if (!cell->at_boundary(f))
+      {
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   cell->neighbor_of_neighbor(f),
+                   numbers::invalid_unsigned_int);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) += fiv.average_gradient(i, qpoint).norm() *
+                              fiv.get_JxW_values()[qpoint];
+        deallog << "average_gradient.norm(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) += fiv.average_hessian(i, qpoint).norm() *
+                              fiv.get_JxW_values()[qpoint];
+        deallog << "average_hessian.norm(): " << cell_vector;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) += fiv.jump_gradient(i, qpoint).norm() *
+                              fiv.get_JxW_values()[qpoint];
+        deallog << "jump_gradient.norm(): " << cell_vector;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(FE_Q<2>(1));
+  test<2>(FE_DGQ<2>(0));
+  test<2>(FE_DGQ<2>(1));
+  test<2>(FE_DGQ<2>(2));
+  test<3>(FE_Q<3>(1));
+  test<3>(FE_DGQ<3>(0));
+  test<3>(FE_DGQ<3>(1));
+}
diff --git a/tests/feinterface/fe_interface_values_08.output b/tests/feinterface/fe_interface_values_08.output
new file mode 100644
index 0000000000..6487816da7
--- /dev/null
+++ b/tests/feinterface/fe_interface_values_08.output
@@ -0,0 +1,50 @@
+
+DEAL::FE_Q<2>(1)
+DEAL::average_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 
+DEAL::average_hessian.norm(): 1.41421 1.41421 1.41421 1.41421 
+DEAL::jump_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 
+DEAL::average_gradient.norm(): 0.250000 1.00000 0.250000 1.00000 0.250000 0.250000 
+DEAL::average_hessian.norm(): 0.707107 0.00000 0.707107 0.00000 0.707107 0.707107 
+DEAL::jump_gradient.norm(): 0.500000 1.00000 0.500000 1.00000 0.500000 0.500000 
+DEAL::FE_DGQ<2>(0)
+DEAL::average_gradient.norm(): 0.00000 
+DEAL::average_hessian.norm(): 0.00000 
+DEAL::jump_gradient.norm(): 0.00000 
+DEAL::average_gradient.norm(): 0.00000 0.00000 
+DEAL::average_hessian.norm(): 0.00000 0.00000 
+DEAL::jump_gradient.norm(): 0.00000 0.00000 
+DEAL::FE_DGQ<2>(1)
+DEAL::average_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 
+DEAL::average_hessian.norm(): 1.41421 1.41421 1.41421 1.41421 
+DEAL::jump_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 
+DEAL::average_gradient.norm(): 0.250000 0.573917 0.250000 0.573917 0.573917 0.250000 0.573917 0.250000 
+DEAL::average_hessian.norm(): 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 
+DEAL::jump_gradient.norm(): 0.500000 1.14783 0.500000 1.14783 1.14783 0.500000 1.14783 0.500000 
+DEAL::FE_DGQ<2>(2)
+DEAL::average_gradient.norm(): 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 
+DEAL::average_hessian.norm(): 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 
+DEAL::jump_gradient.norm(): 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 
+DEAL::average_gradient.norm(): 0.107583 0.430331 0.762102 0.333333 1.33333 1.58962 0.107583 0.430331 0.762102 0.762102 0.430331 0.107583 1.58962 1.33333 0.333333 0.762102 0.430331 0.107583 
+DEAL::average_hessian.norm(): 1.06223 3.84282 3.58699 2.18466 6.72690 6.28519 1.06223 3.84282 3.58699 3.58699 3.84282 1.06223 6.28519 6.72690 2.18466 3.58699 3.84282 1.06223 
+DEAL::jump_gradient.norm(): 0.215166 0.860663 1.52420 0.666667 2.66667 3.17925 0.215166 0.860663 1.52420 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 
+DEAL::FE_Q<3>(1)
+DEAL::average_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 
+DEAL::average_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 
+DEAL::jump_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 
+DEAL::average_gradient.norm(): 0.125000 0.761802 0.125000 0.761802 0.125000 0.761802 0.125000 0.761802 0.125000 0.125000 0.125000 0.125000 
+DEAL::average_hessian.norm(): 0.538675 1.41421 0.538675 1.41421 0.538675 1.41421 0.538675 1.41421 0.538675 0.538675 0.538675 0.538675 
+DEAL::jump_gradient.norm(): 0.250000 0.500000 0.250000 0.500000 0.250000 0.500000 0.250000 0.500000 0.250000 0.250000 0.250000 0.250000 
+DEAL::FE_DGQ<3>(0)
+DEAL::average_gradient.norm(): 0.00000 
+DEAL::average_hessian.norm(): 0.00000 
+DEAL::jump_gradient.norm(): 0.00000 
+DEAL::average_gradient.norm(): 0.00000 0.00000 
+DEAL::average_hessian.norm(): 0.00000 0.00000 
+DEAL::jump_gradient.norm(): 0.00000 0.00000 
+DEAL::FE_DGQ<3>(1)
+DEAL::average_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 
+DEAL::average_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 
+DEAL::jump_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 
+DEAL::average_gradient.norm(): 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 
+DEAL::average_hessian.norm(): 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 
+DEAL::jump_gradient.norm(): 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000