From 3b48b751ee4a72edc0ac166b6861f9f750a878bd Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Thu, 5 Jan 2017 18:07:23 +0100
Subject: [PATCH] tangential jumps for Laplacian

---
 include/deal.II/integrators/laplace.h |  61 ++++++++++-
 tests/integrators/laplacian_02.cc     | 149 ++++++++++++++++++++++++++
 tests/integrators/laplacian_02.output |  23 ++++
 3 files changed, 232 insertions(+), 1 deletion(-)
 create mode 100644 tests/integrators/laplacian_02.cc
 create mode 100644 tests/integrators/laplacian_02.output

diff --git a/include/deal.II/integrators/laplace.h b/include/deal.II/integrators/laplace.h
index 7379965b13..e7a156eac1 100644
--- a/include/deal.II/integrators/laplace.h
+++ b/include/deal.II/integrators/laplace.h
@@ -184,6 +184,65 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * Weak boundary condition of Nitsche type for the Laplacian applied to the
+     * tangential component only, namely on
+     * the face <i>F</i> the matrix
+     * @f[
+     * \int_F \Bigl(\gamma u_\tau v_\tau - \partial_n u_\tau v_\tau - u_\tau \partial_n v_\tau\Bigr)\;ds.
+     * @f]
+     *
+     * Here, $\gamma$ is the <tt>penalty</tt> parameter suitably computed with
+     * compute_penalty().
+     *
+     * @author Guido Kanschat
+     * @date 2017
+     */
+    template <int dim>
+    void nitsche_tangential_matrix (
+      FullMatrix<double> &M,
+      const FEValuesBase<dim> &fe,
+      double penalty,
+      double factor = 1.)
+    {
+      const unsigned int n_dofs = fe.dofs_per_cell;
+      AssertDimension(fe.get_fe().n_components(), dim);
+      Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
+      Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
+
+      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
+        {
+          const double dx = fe.JxW(k) * factor;
+          const Tensor<1,dim> n = fe.normal_vector(k);
+          for (unsigned int i=0; i<n_dofs; ++i)
+            for (unsigned int j=0; j<n_dofs; ++j)
+              {
+                double udotn = 0.;
+                double vdotn = 0.;
+                double ngradun = 0.;
+                double ngradvn = 0.;
+
+                for (unsigned int d=0; d<dim; ++d)
+                  {
+                    udotn += n[d]*fe.shape_value_component(j,k,d);
+                    vdotn += n[d]*fe.shape_value_component(i,k,d);
+                    ngradun += n*fe.shape_grad_component(j,k,d)*n[d];
+                    ngradvn += n*fe.shape_grad_component(i,k,d)*n[d];
+                  }
+
+                for (unsigned int d=0; d<dim; ++d)
+                  {
+                    const double v = fe.shape_value_component(i,k,d)-vdotn*n[d];
+                    const double dnv = n * fe.shape_grad_component(i,k,d)-ngradvn*n[d];
+                    const double u = fe.shape_value_component(j,k,d)-udotn*n[d];
+                    const double dnu = n * fe.shape_grad_component(j,k,d)-ngradun*n[d];
+
+                    M(i,j) += dx *(2.*penalty*u*v-dnu*v-dnv*u);
+                  }
+              }
+        }
+    }
+
     /**
      * Weak boundary condition for the Laplace operator by Nitsche, scalar
      * version, namely on the face <i>F</i> the vector
@@ -428,7 +487,7 @@ namespace LocalIntegrators
                       ngradv2n += n*fe2.shape_grad_component(i,k,d)*n[d];
                     }
 
-                  for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
+                  for (unsigned int d=0; d<dim; ++d)
                     {
                       const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n[d];
                       const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n[d];
diff --git a/tests/integrators/laplacian_02.cc b/tests/integrators/laplacian_02.cc
new file mode 100644
index 0000000000..cb58203449
--- /dev/null
+++ b/tests/integrators/laplacian_02.cc
@@ -0,0 +1,149 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test the tangential functions in integrators/laplace.h
+// Output matrices and assert consistency of residuals
+
+#include "../tests.h"
+#include "../test_grids.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+
+using namespace LocalIntegrators::Laplace;
+
+
+
+template <int dim>
+void test_boundary(const FiniteElement<dim> &fe, bool diff=false)
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  DoFTools::make_zero_boundary_constraints(dof, constraints);
+
+  QGauss<dim-1> quadrature(fe.tensor_degree()+1);
+  FEFaceValues<dim> fev(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values);
+
+  FullMatrix<double> M(fe.dofs_per_cell);
+  FullMatrix<double> Mglobal(dof.n_dofs());
+  std::vector<unsigned int> indices(fe.dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell = dof.begin_active();
+  cell->get_dof_indices(indices);
+  for (unsigned i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+    {
+      fev.reinit(cell, i);
+      nitsche_tangential_matrix(M, fev, 10.);
+      if (diff)
+        nitsche_matrix(M, fev, 10., -1.);
+      constraints.distribute_local_to_global(M, indices, indices, Mglobal);
+    }
+  deallog << fe.get_name() << ": bdry norm " << Mglobal.frobenius_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void test_face(const FiniteElement<dim> &fe, bool diff=false)
+{
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube(tr);
+  tr.refine_global();
+
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  DoFTools::make_zero_boundary_constraints(dof, constraints);
+
+  QGauss<dim-1> quadrature(fe.tensor_degree()+1);
+  FEFaceValues<dim> fev1(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values);
+  FEFaceValues<dim> fev2(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values);
+
+  FullMatrix<double> M11(fe.dofs_per_cell);
+  FullMatrix<double> M12(fe.dofs_per_cell);
+  FullMatrix<double> M21(fe.dofs_per_cell);
+  FullMatrix<double> M22(fe.dofs_per_cell);
+  FullMatrix<double> Mglobal(dof.n_dofs());
+  std::vector<unsigned int> indices1(fe.dofs_per_cell);
+  std::vector<unsigned int> indices2(fe.dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell1 = dof.begin_active();
+  typename DoFHandler<dim>::active_cell_iterator cell2 = ++dof.begin_active();
+
+  cell1->get_dof_indices(indices1);
+  cell2->get_dof_indices(indices2);
+
+  fev1.reinit(cell1, 1);
+  fev2.reinit(cell2, 0);
+  ip_tangential_matrix(M11, M12, M21, M22, fev1, fev2, 10);
+
+  if (diff)
+    ip_matrix(M11, M12, M21, M22, fev1, fev2, 10, -1.);
+
+  constraints.distribute_local_to_global(M11, indices1, indices1, Mglobal);
+  constraints.distribute_local_to_global(M21, indices2, indices1, Mglobal);
+  constraints.distribute_local_to_global(M12, indices1, indices2, Mglobal);
+  constraints.distribute_local_to_global(M22, indices2, indices2, Mglobal);
+  deallog << fe.get_name() << ": face norm " << Mglobal.frobenius_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+  // In each dimension, the first four outputs should be zero: the
+  // tangential jump of Nedelec elements is zero and for
+  // Raviart-Thomas, the tangential and the full jump are the same.
+  FE_Nedelec<dim> n1(2);
+  test_boundary(n1);
+  test_face(n1);
+
+  FE_RaviartThomas<dim> r1(2);
+  test_boundary(r1, true);
+  test_face(r1,true);
+  test_boundary(r1);
+  test_face(r1);
+
+  FE_DGQ<dim> q1(1);
+  FESystem<dim> sys1(q1,dim);
+  test_boundary(sys1);
+  test_boundary(sys1, true);
+  test_face(sys1);
+  test_face(sys1, true);
+  deallog << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/integrators/laplacian_02.output b/tests/integrators/laplacian_02.output
new file mode 100644
index 0000000000..0024ef3b7f
--- /dev/null
+++ b/tests/integrators/laplacian_02.output
@@ -0,0 +1,23 @@
+
+DEAL::FE_Nedelec<2>(2): bdry norm 0
+DEAL::FE_Nedelec<2>(2): face norm 0
+DEAL::FE_RaviartThomas<2>(2): bdry norm 0
+DEAL::FE_RaviartThomas<2>(2): face norm 0
+DEAL::FE_RaviartThomas<2>(2): bdry norm 1143.78
+DEAL::FE_RaviartThomas<2>(2): face norm 259.355
+DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: bdry norm 52.2707
+DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: bdry norm 52.2707
+DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: face norm 4.28174
+DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: face norm 4.28174
+DEAL::
+DEAL::FE_Nedelec<3>(2): bdry norm 0
+DEAL::FE_Nedelec<3>(2): face norm 0
+DEAL::FE_RaviartThomas<3>(2): bdry norm 0
+DEAL::FE_RaviartThomas<3>(2): face norm 0
+DEAL::FE_RaviartThomas<3>(2): bdry norm 5869.16
+DEAL::FE_RaviartThomas<3>(2): face norm 1270.58
+DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: bdry norm 86.7436
+DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: bdry norm 47.9857
+DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: face norm 1.59571
+DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: face norm 1.12834
+DEAL::
-- 
2.39.5