add matrix integrators

git-svn-id: https://svn.dealii.org/trunk@25383 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 e742d031010a0a8d136a6c8eb6f2ac919ba3a483..69783504b6807ff1f40a2f2d7e1c4df9d1140f19 100644 (file)
--- a/deal.II/include/deal.II/integrators/elasticity.h
+++ b/deal.II/include/deal.II/integrators/elasticity.h
@@ -34,10 +34,163 @@ namespace LocalIntegrators
  */
   namespace Elasticity
   {
+/**
+ * Scalar product of symmetric gradients.
+ *
+ * \f[
+ * (\varepsilon(u), \varepsilon(v))
+ * \f]                               
+ */
+    template <int dim>
+    inline void cell_matrix (
+      FullMatrix<double>& M,
+      const FEValuesBase<dim>& fe,
+      const 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 k=0;k<fe.n_quadrature_points;++k)
+       {
+         const double dx = factor * fe.JxW(k);
+         for (unsigned i=0;i<n_dofs;++i)
+           for (unsigned j=0;j<n_dofs;++j)
+             for (unsigned int d1=0;d1<dim;++d1)
+               for (unsigned int d2=0;d2<dim;++d2)
+                 M(i,j) += dx * .25 *
+                           (fe.shape_grad_component(j,k,d1)[d2] + fe.shape_grad_component(j,k,d2)[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
+ * symmetric gradients.
+ */
+    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 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 i=0;i<n_dofs;++i)
+           for (unsigned j=0;j<n_dofs;++j)
+             for (unsigned int d1=0;d1<dim;++d1)
+               {
+                 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;
+                 for (unsigned int d2=0;d2<dim;++d2)
+                   {
+                                                      // v . nabla u n
+                     M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n(d2)* v;
+                                                      // v (nabla u)^T n                     
+                     M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n(d2)* v;
+                                                      // u  nabla v n
+                     M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n(d2)* u;
+                                                      // u (nabla v)^T n                     
+                     M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n(d2)* u;
+                   }
+               }
+       }
+    }
+    
+                                    /**
+                                     * The interior penalty flux
+                                     * for symmetric gradients.
+                                     */
+    template <int dim>
+    inline void ip_matrix (
+      FullMatrix<double>& M11,
+      FullMatrix<double>& M12,
+      FullMatrix<double>& M21,
+      FullMatrix<double>& M22,
+      const FEValuesBase<dim>& fe1,
+      const FEValuesBase<dim>& fe2,
+      const double pen,
+      const double int_factor = 1.,
+      const double ext_factor = -1.)
+    {
+      const unsigned int n_dofs = fe1.dofs_per_cell;
+       
+      AssertDimension(fe1.get_fe().n_components(), dim);
+      AssertDimension(fe2.get_fe().n_components(), dim);
+      AssertDimension(M11.m(), n_dofs);
+      AssertDimension(M11.n(), n_dofs);
+      AssertDimension(M12.m(), n_dofs);
+      AssertDimension(M12.n(), n_dofs);
+      AssertDimension(M21.m(), n_dofs);
+      AssertDimension(M21.n(), n_dofs);
+      AssertDimension(M22.m(), n_dofs);
+      AssertDimension(M22.n(), n_dofs);
+       
+      const double nu1 = int_factor;
+      const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
+      const double penalty = .5 * pen * (nu1 + nu2);
+       
+      for (unsigned k=0;k<fe1.n_quadrature_points;++k)
+       {
+         const double dx = fe1.JxW(k);
+         const Point<dim>& n = fe1.normal_vector(k);
+         for (unsigned i=0;i<n_dofs;++i)
+           for (unsigned j=0;j<n_dofs;++j)
+             for (unsigned int d1=0;d1<dim;++d1)
+               {
+                 const double u1 = fe1.shape_value_component(j,k,d1);
+                 const double u2 = fe2.shape_value_component(j,k,d1);
+                 const double v1 = fe1.shape_value_component(i,k,d1);
+                 const double v2 = fe2.shape_value_component(i,k,d1);
+                   
+                 M11(i,j) += dx * penalty * u1*v1;
+                 M12(i,j) -= dx * penalty * u2*v1;
+                 M21(i,j) -= dx * penalty * u1*v2;
+                 M22(i,j) += dx * penalty * u2*v2;
+                   
+                 for (unsigned int d2=0;d2<dim;++d2)
+                   {
+                                                      // v . nabla u n
+                     M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
+                     M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
+                     M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
+                     M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
+                                                      // v (nabla u)^T n                     
+                     M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
+                     M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
+                     M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
+                     M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
+                                                      // u  nabla v n
+                     M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
+                     M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
+                     M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
+                     M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
+                                                      // u (nabla v)^T n                     
+                     M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
+                     M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
+                     M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
+                     M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
+                   }
+               }
+       }
+    }
   }
 }
 
-
 DEAL_II_NAMESPACE_CLOSE
 
 #endif