include/deal.II/integrators/laplace.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: laplace.h 25370 2012-04-02 20:53:57Z 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_laplace_h
00013 #define __deal2__integrators_laplace_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 Laplace
00036   {
00051     template<int dim>
00052     void cell_matrix (
00053       FullMatrix<double>& M,
00054       const FEValuesBase<dim>& fe,
00055       const double factor = 1.)
00056     {
00057       const unsigned int n_dofs = fe.dofs_per_cell;
00058       const unsigned int n_components = fe.get_fe().n_components();
00059 
00060       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00061         {
00062           const double dx = fe.JxW(k) * factor;
00063           for (unsigned i=0;i<n_dofs;++i)
00064             {
00065               for (unsigned j=0;j<n_dofs;++j)
00066                 for (unsigned int d=0;d<n_components;++d)
00067                   M(i,j) += dx *
00068                             (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
00069             }
00070         }
00071     }
00085     template <int dim>
00086     void nitsche_matrix (
00087       FullMatrix<double>& M,
00088       const FEValuesBase<dim>& fe,
00089       double penalty,
00090       double factor = 1.)
00091     {
00092       const unsigned int n_dofs = fe.dofs_per_cell;
00093       const unsigned int n_comp = fe.get_fe().n_components();
00094 
00095       Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
00096       Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
00097 
00098       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00099         {
00100           const double dx = fe.JxW(k) * factor;
00101           const Point<dim>& n = fe.normal_vector(k);
00102           for (unsigned i=0;i<n_dofs;++i)
00103             for (unsigned j=0;j<n_dofs;++j)
00104               for (unsigned int d=0;d<n_comp;++d)
00105                 M(i,j) += dx *
00106                           (2. * fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
00107                            - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
00108                            - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
00109         }
00110     }
00111 
00130       template <int dim>
00131       void nitsche_residual (
00132         Vector<double>& result,
00133         const FEValuesBase<dim>& fe,
00134         const VectorSlice<const std::vector<std::vector<double> > >& input,
00135         const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput,
00136         const VectorSlice<const std::vector<std::vector<double> > >& data,
00137         double penalty,
00138         double factor = 1.)
00139       {
00140         const unsigned int n_dofs = fe.dofs_per_cell;
00141 
00142         const unsigned int n_comp = fe.get_fe().n_components();
00143         AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
00144         AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points);
00145         AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);
00146 
00147         for (unsigned k=0;k<fe.n_quadrature_points;++k)
00148           {
00149             const double dx = factor * fe.JxW(k);
00150             const Point<dim>& n = fe.normal_vector(k);
00151             for (unsigned i=0;i<n_dofs;++i)
00152               for (unsigned int d=0;d<n_comp;++d)
00153                 {
00154                   const double dnv = fe.shape_grad_component(i,k,d) * n;
00155                   const double dnu = Dinput[d][k] * n;
00156                   const double v= fe.shape_value_component(i,k,d);
00157                   const double u= input[d][k];
00158                   const double g= data[d][k];
00159 
00160                   result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
00161                 }
00162           }
00163       }
00164 
00178       template <int dim>
00179       void nitsche_residual (
00180         Vector<double>& result,
00181         const FEValuesBase<dim>& fe,
00182         const std::vector<double>& input,
00183         const std::vector<Tensor<1,dim> >& Dinput,
00184         const std::vector<double>& data,
00185         double penalty,
00186         double factor = 1.)
00187       {
00188         const unsigned int n_dofs = fe.dofs_per_cell;
00189         AssertDimension(input.size(), fe.n_quadrature_points);
00190         AssertDimension(Dinput.size(), fe.n_quadrature_points);
00191         AssertDimension(data.size(), fe.n_quadrature_points);
00192 
00193         for (unsigned k=0;k<fe.n_quadrature_points;++k)
00194           {
00195             const double dx = factor * fe.JxW(k);
00196             const Point<dim>& n = fe.normal_vector(k);
00197             for (unsigned i=0;i<n_dofs;++i)
00198               {
00199                 const double dnv = fe.shape_grad(i,k) * n;
00200                 const double dnu = Dinput[k] * n;
00201                 const double v= fe.shape_value(i,k);
00202                 const double u= input[k];
00203                 const double g= data[k];
00204 
00205                 result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
00206               }
00207           }
00208       }
00209 
00230     template <int dim>
00231     void ip_matrix (
00232       FullMatrix<double>& M11,
00233       FullMatrix<double>& M12,
00234       FullMatrix<double>& M21,
00235       FullMatrix<double>& M22,
00236       const FEValuesBase<dim>& fe1,
00237       const FEValuesBase<dim>& fe2,
00238       double penalty,
00239       double factor1 = 1.,
00240       double factor2 = -1.)
00241     {
00242       const unsigned int n_dofs = fe1.dofs_per_cell;
00243       AssertDimension(M11.n(), n_dofs);
00244       AssertDimension(M11.m(), n_dofs);
00245       AssertDimension(M12.n(), n_dofs);
00246       AssertDimension(M12.m(), n_dofs);
00247       AssertDimension(M21.n(), n_dofs);
00248       AssertDimension(M21.m(), n_dofs);
00249       AssertDimension(M22.n(), n_dofs);
00250       AssertDimension(M22.m(), n_dofs);
00251 
00252       const double nui = factor1;
00253       const double nue = (factor2 < 0) ? factor1 : factor2;
00254 
00255       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00256         {
00257           const double dx = fe1.JxW(k);
00258           const Point<dim>& n = fe1.normal_vector(k);
00259           for (unsigned int d=0;d<fe1.get_fe().n_components();++d)
00260             {
00261               for (unsigned i=0;i<n_dofs;++i)
00262                 {
00263                   for (unsigned j=0;j<n_dofs;++j)
00264                     {
00265                       const double vi = fe1.shape_value_component(i,k,d);
00266                       const double dnvi = n * fe1.shape_grad_component(i,k,d);
00267                       const double ve = fe2.shape_value_component(i,k,d);
00268                       const double dnve = n * fe2.shape_grad_component(i,k,d);
00269                       const double ui = fe1.shape_value_component(j,k,d);
00270                       const double dnui = n * fe1.shape_grad_component(j,k,d);
00271                       const double ue = fe2.shape_value_component(j,k,d);
00272                       const double dnue = n * fe2.shape_grad_component(j,k,d);
00273                       M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
00274                       M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
00275                       M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
00276                       M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
00277                     }
00278                 }
00279             }
00280         }
00281     }
00282     
00297     template <int dim>
00298     double compute_penalty(
00299       const MeshWorker::DoFInfo<dim>& dinfo1,
00300       const MeshWorker::DoFInfo<dim>& dinfo2,
00301       unsigned int deg1,
00302       unsigned int deg2)
00303     {
00304       const unsigned int normal1 = GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
00305       const unsigned int normal2 = GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
00306       const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1+1);
00307       const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2+1);
00308 
00309       double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
00310       double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
00311       if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
00312         {
00313           Assert (dinfo1.face == dinfo2.face, ExcInternalError());
00314           Assert (dinfo1.face->has_children(), ExcInternalError());
00315           penalty1 *= 2;
00316         }
00317       const double penalty = 0.5*(penalty1 + penalty2);
00318       return penalty;
00319     }
00320   }
00321 }
00322 
00323 
00324 DEAL_II_NAMESPACE_CLOSE
00325 
00326 #endif
00327 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

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