From: Timo Heister <timo.heister@gmail.com>
Date: Tue, 2 Aug 2016 15:30:47 +0000 (+0200)
Subject: add divergence integrators
X-Git-Tag: v8.5.0-rc1~789^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=252aa54b0f99b510bac18db1fd76b897c41976d4;p=dealii.git

add divergence integrators
---

diff --git a/include/deal.II/integrators/divergence.h b/include/deal.II/integrators/divergence.h
index 4bd9c90eb7..fcbc72811b 100644
--- a/include/deal.II/integrators/divergence.h
+++ b/include/deal.II/integrators/divergence.h
@@ -529,6 +529,114 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * The matrix for the weak boundary condition of Nitsche type for linear elasticity:
+     * @f[
+     * \int_F \Bigl(\gamma (u \cdot n)(v \cdot n)  - \nabla\cdot u
+     * v\cdot n - u \cdot n \nabla \cdot v \Bigr)\;ds.
+     * @f]
+     */
+    template <int dim>
+    inline void nitsche_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);
+      AssertDimension(M.m(), n_dofs);
+      AssertDimension(M.n(), n_dofs);
+
+      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
+        {
+          const double dx = factor * fe.JxW(k);
+          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 un = 0., vn = 0., divu = 0., divv = 0.;
+                for (unsigned int d=0; d<dim; ++d)
+                  {
+                    un += fe.shape_value_component(j,k,d) * n[d];
+                    vn += fe.shape_value_component(i,k,d) * n[d];
+                    divu += fe.shape_grad_component(j,k,d)[d];
+                    divv += fe.shape_grad_component(i,k,d)[d];
+                  }
+
+                M(i,j) += dx * 2. * penalty * un * vn;
+                M(i,j) -= dx*(divu*vn+divv*un);
+              }
+        }
+    }
+
+    /**
+     * The interior penalty flux for 
+     */
+    template <int dim>
+    void ip_matrix (
+      FullMatrix<double> &M11,
+      FullMatrix<double> &M12,
+      FullMatrix<double> &M21,
+      FullMatrix<double> &M22,
+      const FEValuesBase<dim> &fe1,
+      const FEValuesBase<dim> &fe2,
+      double penalty,
+      double factor1 = 1.,
+      double factor2 = -1.)
+    {
+      const unsigned int n_dofs = fe1.dofs_per_cell;
+      AssertDimension(M11.n(), n_dofs);
+      AssertDimension(M11.m(), n_dofs);
+      AssertDimension(M12.n(), n_dofs);
+      AssertDimension(M12.m(), n_dofs);
+      AssertDimension(M21.n(), n_dofs);
+      AssertDimension(M21.m(), n_dofs);
+      AssertDimension(M22.n(), n_dofs);
+      AssertDimension(M22.m(), n_dofs);
+
+      const double nui = factor1;
+      const double nue = (factor2 < 0) ? factor1 : factor2;
+      const double nu = .5*(nui+nue);
+
+      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
+        {
+          const double dx = fe1.JxW(k);
+          const Tensor<1,dim> n = fe1.normal_vector(k);
+          for (unsigned int i=0; i<n_dofs; ++i)
+            {
+              for (unsigned int j=0; j<n_dofs; ++j)
+                {
+                  double uni = 0.;
+                  double une = 0.;
+                  double vni = 0.;
+                  double vne = 0.;
+                  double divui = 0.;
+                  double divue = 0.;
+                  double divvi = 0.;
+                  double divve = 0.;
+
+                  for (unsigned int d=0; d<dim; ++d)
+                    {
+                      uni += fe1.shape_value_component(j,k,d) * n[d];
+                      une += fe2.shape_value_component(j,k,d) * n[d];
+                      vni += fe1.shape_value_component(i,k,d) * n[d];
+                      vne += fe2.shape_value_component(i,k,d) * n[d];
+                      divui += fe1.shape_grad_component(j,k,d)[d];
+                      divue += fe2.shape_grad_component(j,k,d)[d];
+                      divvi += fe1.shape_grad_component(i,k,d)[d];
+                      divve += fe2.shape_grad_component(i,k,d)[d];
+                    }
+                  M11(i,j) += dx*(-.5*nui*divvi*uni-.5*nui*divui*vni+nu*penalty*uni*vni);
+                  M12(i,j) += dx*( .5*nui*divvi*une-.5*nue*divue*vni-nu*penalty*vni*une);
+                  M21(i,j) += dx*(-.5*nue*divve*uni+.5*nui*divui*vne-nu*penalty*uni*vne);
+                  M22(i,j) += dx*( .5*nue*divve*une+.5*nue*divue*vne+nu*penalty*une*vne);
+                }
+            }
+        }
+    }
+
     /**
      * The jump of the normal component
      * @f[