From 6b0cc0df8bbae6950027db61ac5f2c3db11e2169 Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 4 Mar 2013 13:13:57 +0000
Subject: [PATCH] add residuals

git-svn-id: https://svn.dealii.org/trunk@28721 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/deal.II/integrators/elasticity.h  | 181 +++++++++++++++++-
 1 file changed, 175 insertions(+), 6 deletions(-)

diff --git a/deal.II/include/deal.II/integrators/elasticity.h b/deal.II/include/deal.II/integrators/elasticity.h
index 29dfc868c7..8fe0766f8b 100644
--- a/deal.II/include/deal.II/integrators/elasticity.h
+++ b/deal.II/include/deal.II/integrators/elasticity.h
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -35,10 +35,11 @@ namespace LocalIntegrators
   namespace Elasticity
   {
     /**
-     * Scalar product of symmetric gradients.
+     * The linear elasticity operator in weak form, namely double
+     * contraction of symmetric gradients.
      *
      * \f[
-     * (\varepsilon(u), \varepsilon(v))
+     * \int_Z \varepsilon(u): \varepsilon(v)\,dx
      * \f]
      */
     template <int dim>
@@ -68,9 +69,48 @@ namespace LocalIntegrators
 
 
     /**
-     * The weak boundary condition
-     * of Nitsche type for
-     * symmetric gradients.
+     * Vector-valued residual operator for linear elasticity in weak form
+     *
+     * \f[
+     * - \int_Z \varepsilon(u): \varepsilon(v) \,dx
+     * \f]
+     */
+    template <int dim>
+    inline void
+    cell_residual  (
+      Vector<double> &result,
+      const FEValuesBase<dim> &fe,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
+      double factor = 1.)
+    {
+      const unsigned int nq = fe.n_quadrature_points;
+      const unsigned int n_dofs = fe.dofs_per_cell;
+      AssertDimension(fe.get_fe().n_components(), dim);
+
+      AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
+      Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+
+      for (unsigned int k=0; k<nq; ++k)
+        {
+          const double dx = factor * fe.JxW(k);
+          for (unsigned int i=0; i<n_dofs; ++i)
+            for (unsigned int d1=0; d1<dim; ++d1)
+	      for (unsigned int d2=0; d2<dim; ++d2)
+		{
+		  result(i) += dx * .25 *
+		    (input[d1][k][d2] + input[d2][k][d1]) *
+		    (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
+		}
+        }
+    }
+
+
+    /**
+     * The weak boundary condition of Nitsche type for linear
+     * elasticity:
+     * @f[
+     * \int_F \Bigl(\gamma (u-g) \cdot v - n^T \epsilon(u) v - (u-g) \epsilon(v) n^T\Bigr)\;ds.
+     * @f]
      */
     template <int dim>
     inline void nitsche_matrix (
@@ -111,6 +151,64 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * Weak boundary condition for the elasticity operator by Nitsche,
+     * namely on the face <i>F</i> the vector
+     * @f[
+     * \int_F \Bigl(\gamma (u-g) \cdot v - n^T \epsilon(u) v - (u-g) \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. <i>g</i> is the inhomogeneous
+     * boundary value in the argument <tt>data</tt>. $n$ is the outer
+     * normal vector and $\gamma$ is the usual penalty parameter.
+     *
+     * @author Guido Kanschat
+     * @date 2013
+     */
+    template <int dim>
+    void nitsche_residual (
+      Vector<double> &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,
+      const VectorSlice<const std::vector<std::vector<double> > > &data,
+      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);
+      AssertVectorVectorDimension(data, 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);
+		const double g= data[d1][k];
+		result(i) += dx + 2.*penalty * (u-g) * 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-g) * fe.shape_grad_component(i,k,d1)[d2] * n(d2);
+		    // u  (nabla v)^T n
+		    result(i) -= .5*dx * (u-g) * fe.shape_grad_component(i,k,d2)[d1] * n(d2);
+		  }
+	      }
+	}
+    }
+    
     /**
      * The interior penalty flux
      * for symmetric gradients.
@@ -188,6 +286,77 @@ namespace LocalIntegrators
                 }
         }
     }
+    /**
+     * Elasticity residual term for the symmetric interior penalty method.
+     *
+     * @author Guido Kanschat
+     * @date 2013
+     */
+    template<int dim>
+    void
+    ip_residual(
+      Vector<double> &result1,
+      Vector<double> &result2,
+      const FEValuesBase<dim> &fe1,
+      const FEValuesBase<dim> &fe2,
+      const VectorSlice<const std::vector<std::vector<double> > > &input1,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput1,
+      const VectorSlice<const std::vector<std::vector<double> > > &input2,
+      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput2,
+      double pen,
+      double int_factor = 1.,
+      double ext_factor = -1.)
+    {
+      const unsigned int n1 = fe1.dofs_per_cell;
+      
+      AssertDimension(fe1.get_fe().n_components(), dim);
+      AssertDimension(fe2.get_fe().n_components(), dim);
+      AssertVectorVectorDimension(input1, dim, fe1.n_quadrature_points);
+      AssertVectorVectorDimension(Dinput1, dim, fe1.n_quadrature_points);
+      AssertVectorVectorDimension(input2, dim, fe2.n_quadrature_points);
+      AssertVectorVectorDimension(Dinput2, dim, fe2.n_quadrature_points);
+
+      const double nu1 = int_factor;
+      const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
+      const double penalty = .5 * pen * (nu1 + nu2);
+
+
+      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
+        {
+          const double dx = fe1.JxW(k);
+          const Point<dim> &n = fe1.normal_vector(k);
+	  
+          for (unsigned int i=0; i<n1; ++i)
+            for (unsigned int d1=0; d1<dim; ++d1)
+              {
+                const double v1 = fe1.shape_value_component(i,k,d1);
+                const double v2 = fe2.shape_value_component(i,k,d1);
+                const double u1 = input1[d1][k];
+                const double u2 = input2[d1][k];
+		
+                result1(i) += dx * penalty * u1*v1;
+		result1(i) -= dx * penalty * u2*v1;
+		result2(i) -= dx * penalty * u1*v2;
+                result2(i) += dx * penalty * u2*v2;
+		
+		for (unsigned int d2=0; d2<dim; ++d2)
+		  {
+		    // v . nabla u n
+		    result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput1[d1][k][d2]) * n(d2) * v1;
+		    result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput1[d1][k][d2]) * n(d2) * v1;
+		    // v . (nabla u)^T n
+		    result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput1[d2][k][d1]) * n(d2) * v1;
+		    result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput1[d2][k][d1]) * n(d2) * v1;
+		    // u  nabla v n
+		    result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
+		    result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * (u1-u2);
+		    // u  (nabla v)^T n
+		    result1(i) -= .25*dx* nu1*fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
+		    result2(i) -= .25*dx* nu2*fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * (u1-u2);
+		  }
+              }
+        }
+    }
   }
 }
 
-- 
2.39.5