From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Tue, 24 Sep 2013 09:04:36 +0000 (+0000)
Subject: add homogeneous boundary operator
X-Git-Tag: v8.1.0~738
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8da6ffaf0877632bb813ebe24d4dd6865cdba768;p=dealii.git

add homogeneous boundary operator

git-svn-id: https://svn.dealii.org/trunk@30907 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/integrators/elasticity.h b/deal.II/include/deal.II/integrators/elasticity.h
index d772cf695c..6623ba9390 100644
--- a/deal.II/include/deal.II/integrators/elasticity.h
+++ b/deal.II/include/deal.II/integrators/elasticity.h
@@ -214,6 +214,60 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * Homogeneous weak boundary condition for the elasticity operator by Nitsche,
+     * namely on the face <i>F</i> the vector
+     * @f[
+     * \int_F \Bigl(\gamma u \cdot v - n^T \epsilon(u) v - u \epsilon(v) n^T\Bigr)\;ds.
+     * @f]
+     *
+     * Here, <i>u</i> is the finite element function whose values and
+     * gradient are given in the arguments <tt>input</tt> and
+     * <tt>Dinput</tt>, respectively. $n$ is the outer
+     * normal vector and $\gamma$ is the usual penalty parameter.
+     *
+     * @author Guido Kanschat
+     * @date 2013
+     */
+    template <int dim, typename number>
+    void nitsche_residual_homogeneous (
+      Vector<number> &result,
+      const FEValuesBase<dim> &fe,
+      const VectorSlice<const std::vector<std::vector<double> > > &input,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput,
+      double penalty,
+      double factor = 1.)
+    {
+      const unsigned int n_dofs = fe.dofs_per_cell;
+      AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
+      AssertVectorVectorDimension(Dinput, dim, fe.n_quadrature_points);
+
+      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
+        {
+          const double dx = factor * fe.JxW(k);
+          const Point<dim> &n = fe.normal_vector(k);
+          for (unsigned int i=0; i<n_dofs; ++i)
+            for (unsigned int d1=0; d1<dim; ++d1)
+              {
+                const double u= input[d1][k];
+                const double v= fe.shape_value_component(i,k,d1);
+                result(i) += dx + 2.*penalty * u * v;
+
+                for (unsigned int d2=0; d2<dim; ++d2)
+                  {
+                    // v . nabla u n
+                    result(i) -= .5*dx* v * Dinput[d1][k][d2] * n(d2);
+                    // v . (nabla u)^T n
+                    result(i) -= .5*dx* v * Dinput[d2][k][d1] * n(d2);
+                    // u  nabla v n
+                    result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
+                    // u  (nabla v)^T n
+                    result(i) -= .5*dx * u * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
+                  }
+              }
+        }
+    }
+
     /**
      * The interior penalty flux
      * for symmetric gradients.