]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add residuals
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 4 Mar 2013 13:13:57 +0000 (13:13 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 4 Mar 2013 13:13:57 +0000 (13:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@28721 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/integrators/elasticity.h

index 29dfc868c7de80a9c23f6d68260164bc95c12498..8fe0766f8bd8aebfae4a17be733386cbc6fb8f35 100644 (file)
@@ -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);
+                 }
+              }
+        }
+    }
   }
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.