]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement jump_hessian and jump_3rd_derivative on FEInterfaceValues. 9259/head
authorSimon Sticko <simon@sticko.se>
Wed, 8 Jan 2020 13:37:47 +0000 (14:37 +0100)
committerSimon Sticko <simon@sticko.se>
Wed, 8 Jan 2020 13:45:15 +0000 (14:45 +0100)
In the same way as the other jump_* functions.

include/deal.II/fe/fe_interface_values.h
tests/feinterface/fe_interface_values_09.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_09.output [new file with mode: 0644]

index 4b3b173425649a89ae8f2ac7bb8f1978929b62be..b6923b0ddcf782ba4b828afad8f5457dd3201938 100644 (file)
@@ -371,6 +371,34 @@ public:
                 const unsigned int q_point,
                 const unsigned int component = 0) const;
 
+  /**
+   * Return the jump in the Hessian $[\nabla^2 u] = \nabla^2 u_{\text{cell0}} -
+   * \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>
+  jump_hessian(const unsigned int interface_dof_index,
+               const unsigned int q_point,
+               const unsigned int component = 0) const;
+
+  /**
+   * Return the jump in the third derivative $[\nabla^3 u] = \nabla^3
+   * u_{\text{cell0}} - \nabla^3 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^3 u] = \nabla^3 u_{\text{cell0}}$.
+   */
+  Tensor<3, dim>
+  jump_3rd_derivative(const unsigned int interface_dof_index,
+                      const unsigned int q_point,
+                      const unsigned int component = 0) const;
+
   /**
    * @}
    */
@@ -872,6 +900,67 @@ FEInterfaceValues<dim, spacedim>::jump_gradient(
 }
 
 
+
+template <int dim, int spacedim>
+Tensor<2, dim>
+FEInterfaceValues<dim, spacedim>::jump_hessian(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dof_pair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
+                                                         q_point,
+                                                         component);
+
+  Tensor<2, dim> value;
+
+  if (dof_pair[0] != numbers::invalid_unsigned_int)
+    value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
+                                                                 q_point,
+                                                                 component);
+  if (dof_pair[1] != numbers::invalid_unsigned_int)
+    value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
+                                                                  q_point,
+                                                                  component);
+
+  return value;
+}
+
+
+template <int dim, int spacedim>
+Tensor<3, dim>
+FEInterfaceValues<dim, spacedim>::jump_3rd_derivative(
+  const unsigned int interface_dof_index,
+  const unsigned int q_point,
+  const unsigned int component) const
+{
+  const auto dof_pair = dofmap[interface_dof_index];
+
+  if (at_boundary())
+    return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
+                                                                q_point,
+                                                                component);
+
+  Tensor<3, dim> value;
+
+  if (dof_pair[0] != numbers::invalid_unsigned_int)
+    value +=
+      1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
+                                                                 q_point,
+                                                                 component);
+  if (dof_pair[1] != numbers::invalid_unsigned_int)
+    value +=
+      -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1],
+                                                                  q_point,
+                                                                  component);
+
+  return value;
+}
+
+
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/feinterface/fe_interface_values_09.cc b/tests/feinterface/fe_interface_values_09.cc
new file mode 100644 (file)
index 0000000..55766fb
--- /dev/null
@@ -0,0 +1,150 @@
+// ---------------------------------------------------------------------
+//
+// 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 jump_hessian and jump_third_derivative 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
+print_norm_of_average_over_quadrature_points(const FEInterfaceValues<dim> &fiv)
+{
+  const unsigned int n_dofs = fiv.n_current_interface_dofs();
+  Vector<double>     cell_vector(n_dofs);
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < fiv.n_quadrature_points; ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) +=
+        fiv.jump_hessian(i, qpoint).norm() * fiv.get_JxW_values()[qpoint];
+  deallog << "jump_hessian.norm(): " << cell_vector;
+
+  cell_vector = 0.0;
+  for (unsigned int qpoint = 0; qpoint < fiv.n_quadrature_points; ++qpoint)
+    for (unsigned int i = 0; i < n_dofs; ++i)
+      cell_vector(i) += fiv.jump_3rd_derivative(i, qpoint).norm() *
+                        fiv.get_JxW_values()[qpoint];
+  deallog << "jump_3rd_derivative.norm(): " << cell_vector;
+}
+
+
+
+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);
+  const UpdateFlags update_flags =
+    update_hessians | update_3rd_derivatives | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+  auto cell = dofh.begin();
+
+  // Print jump in Hessian and third derivative for a face at the boundary.
+  {
+    unsigned int face = 0;
+    Assert(cell->at_boundary(face), ExcInternalError());
+
+    fiv.reinit(cell, face);
+
+    print_norm_of_average_over_quadrature_points(fiv);
+  }
+
+  // Print jump in Hessian and third derivative for a face between two cells.
+  {
+    const unsigned int face = 1;
+    Assert(!cell->at_boundary(face), ExcInternalError());
+
+    fiv.reinit(cell,
+               face,
+               numbers::invalid_unsigned_int,
+               cell->neighbor(face),
+               cell->neighbor_of_neighbor(face),
+               numbers::invalid_unsigned_int);
+
+    print_norm_of_average_over_quadrature_points(fiv);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  // Test the lowest order cg and dg elements which have a non-zero third
+  // derivative.
+  test<2>(FE_Q<2>(2));
+  test<2>(FE_DGQ<2>(2));
+  test<3>(FE_Q<3>(1));
+  test<3>(FE_DGQ<3>(1));
+}
diff --git a/tests/feinterface/fe_interface_values_09.output b/tests/feinterface/fe_interface_values_09.output
new file mode 100644 (file)
index 0000000..904fbfa
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::FE_Q<2>(2)
+DEAL::jump_hessian.norm(): 7.17398 2.12446 7.17398 2.12446 12.5704 4.36931 7.68564 7.68564 13.4538 
+DEAL::jump_3rd_derivative.norm(): 23.1831 11.8202 23.1831 11.8202 44.4667 20.3528 34.2248 34.2248 63.5828 
+DEAL::jump_hessian.norm(): 2.12446 11.0742 2.12446 11.0742 4.36931 14.6059 7.68564 7.68564 13.4538 2.12446 2.12446 4.36931 7.68564 7.68564 13.4538 
+DEAL::jump_3rd_derivative.norm(): 11.8202 41.5692 11.8202 41.5692 20.3528 83.1384 34.2248 34.2248 63.5828 11.8202 11.8202 20.3528 34.2248 34.2248 63.5828 
+DEAL::FE_DGQ<2>(2)
+DEAL::jump_hessian.norm(): 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 
+DEAL::jump_3rd_derivative.norm(): 23.1831 34.2248 11.8202 44.4667 63.5828 20.3528 23.1831 34.2248 11.8202 
+DEAL::jump_hessian.norm(): 2.12446 7.68564 7.17398 4.36931 13.4538 12.5704 2.12446 7.68564 7.17398 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 
+DEAL::jump_3rd_derivative.norm(): 11.8202 34.2248 23.1831 20.3528 63.5828 44.4667 11.8202 34.2248 23.1831 23.1831 34.2248 11.8202 44.4667 63.5828 20.3528 23.1831 34.2248 11.8202 
+DEAL::FE_Q<3>(1)
+DEAL::jump_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 
+DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 
+DEAL::jump_hessian.norm(): 1.07735 2.15470 1.07735 2.15470 1.07735 2.15470 1.07735 2.15470 1.07735 1.07735 1.07735 1.07735 
+DEAL::jump_3rd_derivative.norm(): 2.44949 4.89898 2.44949 4.89898 2.44949 4.89898 2.44949 4.89898 2.44949 2.44949 2.44949 2.44949 
+DEAL::FE_DGQ<3>(1)
+DEAL::jump_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 
+DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 
+DEAL::jump_hessian.norm(): 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 
+DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.