From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Thu, 26 Sep 2013 18:42:28 +0000 (+0000)
Subject: new residual Divergence::u_times_n\
X-Git-Tag: v8.1.0~716
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6c92e62df8580aa43fe7b4605443339acd6bfdd5;p=dealii.git

new residual Divergence::u_times_n\
fixed elasticity residuals

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

diff --git a/deal.II/include/deal.II/integrators/divergence.h b/deal.II/include/deal.II/integrators/divergence.h
index ef99627c83..836936425b 100644
--- a/deal.II/include/deal.II/integrators/divergence.h
+++ b/deal.II/include/deal.II/integrators/divergence.h
@@ -380,6 +380,43 @@ namespace LocalIntegrators
         }
     }
 
+    /**
+     * The trace of the gradient
+     * operator, namely the product
+     * of the normal component of the
+     * vector valued test space and
+     * the trial space.
+     * @f[
+     * \int_F u (\mathbf v\cdot \mathbf n) \,ds
+     * @f]
+     *
+     * @author Guido Kanschat
+     * @date 2013
+     */
+    template<int dim, typename number>
+    void
+    u_times_n_residual (
+      Vector<number> &result,
+      const FEValuesBase<dim> &fetest,
+      const std::vector<double> &data,
+      double factor = 1.)
+    {
+      const unsigned int t_dofs = fetest.dofs_per_cell;
+
+      AssertDimension(fetest.get_fe().n_components(), dim);
+      AssertDimension(result.size(), t_dofs);
+      AssertDimension(data.size(), fetest.n_quadrature_points);
+
+      for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
+        {
+          const Tensor<1,dim> ndx = factor * fetest.normal_vector(k) * fetest.JxW(k);
+
+          for (unsigned int i=0; i<t_dofs; ++i)
+            for (unsigned int d=0; d<dim; ++d)
+              result(i) += ndx[d] * fetest.shape_value_component(i,k,d) * data[k];
+        }
+    }
+
     /**
      * The trace of the divergence
      * operator, namely the product
diff --git a/deal.II/include/deal.II/integrators/elasticity.h b/deal.II/include/deal.II/integrators/elasticity.h
index b38480189e..a02e161f5c 100644
--- a/deal.II/include/deal.II/integrators/elasticity.h
+++ b/deal.II/include/deal.II/integrators/elasticity.h
@@ -140,7 +140,7 @@ namespace LocalIntegrators
                 {
                   const double u = fe.shape_value_component(j,k,d1);
                   const double v = fe.shape_value_component(i,k,d1);
-                  M(i,j) += dx * penalty * u * v;
+                  M(i,j) += dx * 2. * penalty * u * v;
                   for (unsigned int d2=0; d2<dim; ++d2)
                     {
                       // v . nabla u n
@@ -197,7 +197,7 @@ namespace LocalIntegrators
                 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;
+                result(i) += dx * 2.*penalty * (u-g) * v;
 
                 for (unsigned int d2=0; d2<dim; ++d2)
                   {
@@ -401,11 +401,11 @@ namespace LocalIntegrators
                 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;
+                    result1(i) -= .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v1;
+                    result2(i) += .25*dx* (nu1*Dinput1[d1][k][d2]+nu2*Dinput2[d1][k][d2]) * n(d2) * v2;
                     // 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;
+                    result1(i) -= .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v1;
+                    result2(i) += .25*dx* (nu1*Dinput1[d2][k][d1]+nu2*Dinput2[d2][k][d1]) * n(d2) * v2;
                     // 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);