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 ef99627c83c1ce4e5ea584e330c2779c94d2aba4..836936425bee7f6c151e93db461e4f8eb1e01175 100644 (file)
--- 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 b38480189ed080b80e774cabe91c053a65716fc3..a02e161f5c07d84438a3a39be518ee49f7364aa9 100644 (file)
--- 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);