include/deal.II/integrators/divergence.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: divergence.h 25345 2012-03-31 08:37:04Z bangerth @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_divergence_h
00013 #define __deal2__integrators_divergence_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 Divergence
00036   {
00047     template <int dim>
00048     Tensor<1,dim>
00049     grad_div(
00050       const Tensor<2,dim>& h0,
00051       const Tensor<2,dim>& h1,
00052       const Tensor<2,dim>& h2)
00053     {
00054       Tensor<1,dim> result;
00055       for (unsigned int d=0;d<dim;++d)
00056         {
00057           result[d] += h0[d][0];
00058           if (dim >=2) result[d] += h1[d][1];
00059           if (dim >=3) result[d] += h2[d][2];
00060         }
00061       return result;
00062     }
00063 
00064 
00075     template <int dim>
00076     void cell_matrix (
00077       FullMatrix<double>& M,
00078       const FEValuesBase<dim>& fe,
00079       const FEValuesBase<dim>& fetest,
00080       double factor = 1.)
00081     {
00082       unsigned int fecomp = fe.get_fe().n_components();
00083       const unsigned int n_dofs = fe.dofs_per_cell;
00084       const unsigned int t_dofs = fetest.dofs_per_cell;
00085       AssertDimension(fecomp, dim);
00086       AssertDimension(fetest.get_fe().n_components(), 1);
00087       AssertDimension(M.m(), t_dofs);
00088       AssertDimension(M.n(), n_dofs);
00089 
00090       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00091         {
00092           const double dx = fe.JxW(k) * factor;
00093           for (unsigned i=0;i<t_dofs;++i)
00094             {
00095               const double vv = fetest.shape_value(i,k);
00096               for (unsigned int d=0;d<dim;++d)
00097                 for (unsigned j=0;j<n_dofs;++j)
00098                   {
00099                     const double du = fe.shape_grad_component(j,k,d)[d];
00100                     M(i,j) += dx * du * vv;
00101                   }
00102             }
00103         }
00104     }
00115     template <int dim>
00116     void gradient_matrix(
00117       FullMatrix<double>& M,
00118       const FEValuesBase<dim>& fe,
00119       const FEValuesBase<dim>& fetest,
00120       double factor = 1.)
00121     {
00122       unsigned int fecomp = fetest.get_fe().n_components();
00123       const unsigned int t_dofs = fetest.dofs_per_cell;
00124       const unsigned int n_dofs = fe.dofs_per_cell;
00125 
00126       AssertDimension(fecomp, dim);
00127       AssertDimension(fe.get_fe().n_components(), 1);
00128       AssertDimension(M.m(), t_dofs);
00129       AssertDimension(M.n(), n_dofs);
00130 
00131       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00132         {
00133           const double dx = fe.JxW(k) * factor;
00134           for (unsigned int d=0;d<dim;++d)
00135             for (unsigned i=0;i<t_dofs;++i)
00136               {
00137                 const double vv = fetest.shape_value_component(i,k,d);
00138                 for (unsigned j=0;j<n_dofs;++j)
00139                   {
00140                     const Tensor<1,dim>& Du = fe.shape_grad(j,k);
00141                     M(i,j) += dx * vv * Du[d];
00142                   }
00143               }
00144         }
00145     }
00146 
00157     template<int dim>
00158     void
00159     u_dot_n_matrix (
00160       FullMatrix<double>& M,
00161       const FEValuesBase<dim>& fe,
00162       const FEValuesBase<dim>& fetest,
00163       double factor = 1.)
00164     {
00165       unsigned int fecomp = fe.get_fe().n_components();
00166       const unsigned int n_dofs = fe.dofs_per_cell;
00167       const unsigned int t_dofs = fetest.dofs_per_cell;
00168 
00169       AssertDimension(fecomp, dim);
00170       AssertDimension(fetest.get_fe().n_components(), 1);
00171       AssertDimension(M.m(), t_dofs);
00172       AssertDimension(M.n(), n_dofs);
00173 
00174       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00175         {
00176           const Tensor<1,dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
00177           for (unsigned i=0;i<t_dofs;++i)
00178             for (unsigned j=0;j<n_dofs;++j)
00179               for (unsigned int d=0;d<dim;++d)
00180                 M(i,j) += ndx[d] * fe.shape_value_component(j,k,d)
00181                           * fetest.shape_value(i,k);
00182         }
00183     }
00184 
00195     template<int dim>
00196     void
00197     u_dot_n_residual (
00198       Vector<double>& result,
00199       const FEValuesBase<dim>& fe,
00200       const FEValuesBase<dim>& fetest,
00201       const VectorSlice<const std::vector<std::vector<double> > >& data,
00202       double factor = 1.)
00203     {
00204       unsigned int fecomp = fe.get_fe().n_components();
00205       const unsigned int t_dofs = fetest.dofs_per_cell;
00206 
00207       AssertDimension(fecomp, dim);
00208       AssertDimension(fetest.get_fe().n_components(), 1);
00209       AssertDimension(result.size(), t_dofs);
00210       AssertVectorVectorDimension (data, dim, fe.n_quadrature_points);
00211 
00212       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00213         {
00214           const Tensor<1,dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
00215 
00216           for (unsigned i=0;i<t_dofs;++i)
00217             for (unsigned int d=0;d<dim;++d)
00218               result(i) += ndx[d] * fetest.shape_value(i,k) * data[d][k];
00219         }
00220     }
00231     template<int dim>
00232     void
00233     u_dot_n_matrix (
00234       FullMatrix<double>& M11,
00235       FullMatrix<double>& M12,
00236       FullMatrix<double>& M21,
00237       FullMatrix<double>& M22,
00238       const FEValuesBase<dim>& fe1,
00239       const FEValuesBase<dim>& fe2,
00240       const FEValuesBase<dim>& fetest1,
00241       const FEValuesBase<dim>& fetest2,
00242       double factor = 1.)
00243     {
00244       const unsigned int n_dofs = fe1.dofs_per_cell;
00245       const unsigned int t_dofs = fetest1.dofs_per_cell;
00246 
00247       AssertDimension(fe1.get_fe().n_components(), dim);
00248       AssertDimension(fe2.get_fe().n_components(), dim);
00249       AssertDimension(fetest1.get_fe().n_components(), 1);
00250       AssertDimension(fetest2.get_fe().n_components(), 1);
00251       AssertDimension(M11.m(), t_dofs);
00252       AssertDimension(M11.n(), n_dofs);
00253       AssertDimension(M12.m(), t_dofs);
00254       AssertDimension(M12.n(), n_dofs);
00255       AssertDimension(M21.m(), t_dofs);
00256       AssertDimension(M21.n(), n_dofs);
00257       AssertDimension(M22.m(), t_dofs);
00258       AssertDimension(M22.n(), n_dofs);
00259 
00260       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00261         {
00262           const double dx = factor * fe1.JxW(k);
00263           for (unsigned i=0;i<t_dofs;++i)
00264             for (unsigned j=0;j<n_dofs;++j)
00265               for (unsigned int d=0;d<dim;++d)
00266                 {
00267                   const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
00268                   const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
00269                   const double v1 = fetest1.shape_value(i,k);
00270                   const double v2 = fetest2.shape_value(i,k);
00271 
00272                   M11(i,j) += .5 * dx * un1 * v1;
00273                   M12(i,j) += .5 * dx * un2 * v1;
00274                   M21(i,j) += .5 * dx * un1 * v2;
00275                   M22(i,j) += .5 * dx * un2 * v2;
00276                 }
00277         }
00278     }
00279 
00286     template <int dim>
00287     void grad_div_matrix (
00288       FullMatrix<double>& M,
00289       const FEValuesBase<dim>& fe,
00290       double factor = 1.)
00291     {
00292       const unsigned int n_dofs = fe.dofs_per_cell;
00293 
00294       AssertDimension(fe.get_fe().n_components(), dim);
00295       AssertDimension(M.m(), n_dofs);
00296       AssertDimension(M.n(), n_dofs);
00297 
00298       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00299         {
00300           const double dx = factor * fe.JxW(k);
00301           for (unsigned i=0;i<n_dofs;++i)
00302             for (unsigned j=0;j<n_dofs;++j)
00303               {
00304                 double dv = 0.;
00305                 double du = 0.;
00306                 for (unsigned int d=0;d<dim;++d)
00307                   {
00308                     dv += fe.shape_grad_component(i,k,d)[d];
00309                     du += fe.shape_grad_component(j,k,d)[d];
00310                   }
00311 
00312                 M(i,j) += dx * du * dv;
00313               }
00314         }
00315     }
00316 
00326     template<int dim>
00327     void
00328     u_dot_n_jump_matrix (
00329       FullMatrix<double>& M11,
00330       FullMatrix<double>& M12,
00331       FullMatrix<double>& M21,
00332       FullMatrix<double>& M22,
00333       const FEValuesBase<dim>& fe1,
00334       const FEValuesBase<dim>& fe2,
00335       double factor = 1.)
00336     {
00337       const unsigned int n_dofs = fe1.dofs_per_cell;
00338 
00339       AssertDimension(fe1.get_fe().n_components(), dim);
00340       AssertDimension(fe2.get_fe().n_components(), dim);
00341       AssertDimension(M11.m(), n_dofs);
00342       AssertDimension(M11.n(), n_dofs);
00343       AssertDimension(M12.m(), n_dofs);
00344       AssertDimension(M12.n(), n_dofs);
00345       AssertDimension(M21.m(), n_dofs);
00346       AssertDimension(M21.n(), n_dofs);
00347       AssertDimension(M22.m(), n_dofs);
00348       AssertDimension(M22.n(), n_dofs);
00349 
00350       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00351         {
00352           const double dx = factor * fe1.JxW(k);
00353           for (unsigned i=0;i<n_dofs;++i)
00354             for (unsigned j=0;j<n_dofs;++j)
00355               for (unsigned int d=0;d<dim;++d)
00356                 {
00357                   const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
00358                   const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
00359                   const double vn1 = fe1.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
00360                   const double vn2 =-fe2.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
00361 
00362                   M11(i,j) += dx * un1 * vn1;
00363                   M12(i,j) += dx * un2 * vn1;
00364                   M21(i,j) += dx * un1 * vn2;
00365                   M22(i,j) += dx * un2 * vn2;
00366                 }
00367         }
00368     }
00369 
00370     template <int dim>
00371     double norm(const FEValuesBase<dim>& fe,
00372                 const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Du)
00373     {
00374       unsigned int fecomp = fe.get_fe().n_components();
00375 
00376       AssertDimension(fecomp, dim);
00377       AssertVectorVectorDimension (Du, dim, fe.n_quadrature_points);
00378 
00379       double result = 0;
00380       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00381         {
00382           double div = Du[0][k][0];
00383           for (unsigned int d=1;d<dim;++d)
00384             div += Du[d][k][d];
00385           result += div*div;
00386         }
00387       return result;
00388     }
00389 
00390   }
00391 }
00392 
00393 
00394 DEAL_II_NAMESPACE_CLOSE
00395 
00396 #endif
00397 
 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