include/deal.II/integrators/elasticity.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: elasticity.h 25383 2012-04-04 19:37:47Z kanschat @f$
00003 //
00004 //    Copyright (C) 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__integrators_elasticity_h
00013 #define __deal2__integrators_elasticity_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/quadrature.h>
00019 #include <deal.II/lac/full_matrix.h>
00020 #include <deal.II/fe/mapping.h>
00021 #include <deal.II/fe/fe_values.h>
00022 #include <deal.II/meshworker/dof_info.h>
00023 
00024 DEAL_II_NAMESPACE_OPEN
00025 
00026 namespace LocalIntegrators
00027 {
00035   namespace Elasticity
00036   {
00044     template <int dim>
00045     inline void cell_matrix (
00046       FullMatrix<double>& M,
00047       const FEValuesBase<dim>& fe,
00048       const double factor = 1.)
00049     {
00050       const unsigned int n_dofs = fe.dofs_per_cell;
00051       
00052       AssertDimension(fe.get_fe().n_components(), dim);
00053       AssertDimension(M.m(), n_dofs);
00054       AssertDimension(M.n(), n_dofs);
00055       
00056       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00057         {
00058           const double dx = factor * fe.JxW(k);
00059           for (unsigned i=0;i<n_dofs;++i)
00060             for (unsigned j=0;j<n_dofs;++j)
00061               for (unsigned int d1=0;d1<dim;++d1)
00062                 for (unsigned int d2=0;d2<dim;++d2)
00063                   M(i,j) += dx * .25 *
00064                             (fe.shape_grad_component(j,k,d1)[d2] + fe.shape_grad_component(j,k,d2)[d1]) *
00065                             (fe.shape_grad_component(i,k,d1)[d2] + fe.shape_grad_component(i,k,d2)[d1]);
00066         }
00067     }
00068 
00069 
00075     template <int dim>
00076     inline void nitsche_matrix (
00077       FullMatrix<double>& M,
00078       const FEValuesBase<dim>& fe,
00079       double penalty,
00080       double factor = 1.)
00081     {
00082       const unsigned int n_dofs = fe.dofs_per_cell;
00083       
00084       AssertDimension(fe.get_fe().n_components(), dim);
00085       AssertDimension(M.m(), n_dofs);
00086       AssertDimension(M.n(), n_dofs);
00087       
00088       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00089         {
00090           const double dx = factor * fe.JxW(k);
00091           const Point<dim>& n = fe.normal_vector(k);
00092           for (unsigned i=0;i<n_dofs;++i)
00093             for (unsigned j=0;j<n_dofs;++j)
00094               for (unsigned int d1=0;d1<dim;++d1)
00095                 {
00096                   const double u = fe.shape_value_component(j,k,d1);
00097                   const double v = fe.shape_value_component(i,k,d1);
00098                   M(i,j) += dx * penalty * u * v;
00099                   for (unsigned int d2=0;d2<dim;++d2)
00100                     {
00101                                                        // v . nabla u n
00102                       M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d1)[d2] *n(d2)* v;
00103                                                        // v (nabla u)^T n                     
00104                       M(i,j) -= .5*dx* fe.shape_grad_component(j,k,d2)[d1] *n(d2)* v;
00105                                                        // u  nabla v n
00106                       M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d1)[d2] *n(d2)* u;
00107                                                        // u (nabla v)^T n                     
00108                       M(i,j) -= .5*dx* fe.shape_grad_component(i,k,d2)[d1] *n(d2)* u;
00109                     }
00110                 }
00111         }
00112     }
00113     
00118     template <int dim>
00119     inline void ip_matrix (
00120       FullMatrix<double>& M11,
00121       FullMatrix<double>& M12,
00122       FullMatrix<double>& M21,
00123       FullMatrix<double>& M22,
00124       const FEValuesBase<dim>& fe1,
00125       const FEValuesBase<dim>& fe2,
00126       const double pen,
00127       const double int_factor = 1.,
00128       const double ext_factor = -1.)
00129     {
00130       const unsigned int n_dofs = fe1.dofs_per_cell;
00131         
00132       AssertDimension(fe1.get_fe().n_components(), dim);
00133       AssertDimension(fe2.get_fe().n_components(), dim);
00134       AssertDimension(M11.m(), n_dofs);
00135       AssertDimension(M11.n(), n_dofs);
00136       AssertDimension(M12.m(), n_dofs);
00137       AssertDimension(M12.n(), n_dofs);
00138       AssertDimension(M21.m(), n_dofs);
00139       AssertDimension(M21.n(), n_dofs);
00140       AssertDimension(M22.m(), n_dofs);
00141       AssertDimension(M22.n(), n_dofs);
00142         
00143       const double nu1 = int_factor;
00144       const double nu2 = (ext_factor < 0) ? int_factor : ext_factor;
00145       const double penalty = .5 * pen * (nu1 + nu2);
00146         
00147       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00148         {
00149           const double dx = fe1.JxW(k);
00150           const Point<dim>& n = fe1.normal_vector(k);
00151           for (unsigned i=0;i<n_dofs;++i)
00152             for (unsigned j=0;j<n_dofs;++j)
00153               for (unsigned int d1=0;d1<dim;++d1)
00154                 {
00155                   const double u1 = fe1.shape_value_component(j,k,d1);
00156                   const double u2 = fe2.shape_value_component(j,k,d1);
00157                   const double v1 = fe1.shape_value_component(i,k,d1);
00158                   const double v2 = fe2.shape_value_component(i,k,d1);
00159                     
00160                   M11(i,j) += dx * penalty * u1*v1;
00161                   M12(i,j) -= dx * penalty * u2*v1;
00162                   M21(i,j) -= dx * penalty * u1*v2;
00163                   M22(i,j) += dx * penalty * u2*v2;
00164                     
00165                   for (unsigned int d2=0;d2<dim;++d2)
00166                     {
00167                                                        // v . nabla u n
00168                       M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
00169                       M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v1;
00170                       M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
00171                       M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d1)[d2] * n(d2) * v2;
00172                                                        // v (nabla u)^T n                     
00173                       M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
00174                       M12(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v1;
00175                       M21(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
00176                       M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(j,k,d2)[d1] * n(d2) * v2;
00177                                                        // u  nabla v n
00178                       M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
00179                       M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
00180                       M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u1;
00181                       M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d1)[d2] * n(d2) * u2;
00182                                                        // u (nabla v)^T n                     
00183                       M11(i,j) -= .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
00184                       M12(i,j) += .25 * dx * nu1 * fe1.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
00185                       M21(i,j) -= .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u1;
00186                       M22(i,j) += .25 * dx * nu2 * fe2.shape_grad_component(i,k,d2)[d1] * n(d2) * u2;
00187                     }
00188                 }
00189         }
00190     }
00191   }
00192 }
00193 
00194 DEAL_II_NAMESPACE_CLOSE
00195 
00196 #endif
00197 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:09 by doxygen 1.7.3