include/deal.II/integrators/l2.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: l2.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_l2_h
00013 #define __deal2__integrators_l2_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 L2
00036   {
00052     template <int dim>
00053     void mass_matrix (
00054       FullMatrix<double>& M,
00055       const FEValuesBase<dim>& fe,
00056       const double factor = 1.)
00057     {
00058       const unsigned int n_dofs = fe.dofs_per_cell;
00059       const unsigned int n_components = fe.get_fe().n_components();
00060 
00061       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00062         {
00063           const double dx = fe.JxW(k) * factor;
00064 
00065           for (unsigned i=0;i<n_dofs;++i)
00066             for (unsigned j=0;j<n_dofs;++j)
00067               for (unsigned int d=0;d<n_components;++d)
00068                 M(i,j) += dx
00069                           * fe.shape_value_component(j,k,d)
00070                           * fe.shape_value_component(i,k,d);
00071         }
00072     }
00073 
00081     template <int dim>
00082     void L2 (
00083       Vector<double>& result,
00084       const FEValuesBase<dim>& fe,
00085       const std::vector<double>& input,
00086       const double factor = 1.)
00087     {
00088       const unsigned int n_dofs = fe.dofs_per_cell;
00089       AssertDimension(result.size(), n_dofs);
00090       AssertDimension(fe.get_fe().n_components(), 1);
00091       AssertDimension(input.size(), fe.n_quadrature_points);
00092 
00093       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00094         for (unsigned i=0;i<n_dofs;++i)
00095           result(i) += fe.JxW(k) * factor * input[k] * fe.shape_value(i,k);
00096     }
00097 
00111     template <int dim>
00112     void L2 (
00113       Vector<double>& result,
00114       const FEValuesBase<dim>& fe,
00115       const VectorSlice<const std::vector<std::vector<double> > >& input,
00116       const double factor = 1.)
00117     {
00118       const unsigned int n_dofs = fe.dofs_per_cell;
00119       const unsigned int fe_components = fe.get_fe().n_components();
00120       const unsigned int n_components = input.size();
00121 
00122       AssertDimension(result.size(), n_dofs);
00123       AssertDimension(input.size(), fe_components);
00124 
00125       for (unsigned k=0;k<fe.n_quadrature_points;++k)
00126         for (unsigned i=0;i<n_dofs;++i)
00127           for (unsigned int d=0;d<n_components;++d)
00128             result(i) += fe.JxW(k) * factor * fe.shape_value_component(i,k,d) * input[d][k];
00129     }
00130 
00148     template <int dim>
00149     void jump_matrix (
00150       FullMatrix<double>& M11,
00151       FullMatrix<double>& M12,
00152       FullMatrix<double>& M21,
00153       FullMatrix<double>& M22,
00154       const FEValuesBase<dim>& fe1,
00155       const FEValuesBase<dim>& fe2,
00156       const double factor1 = 1.,
00157       const double factor2 = 1.)
00158     {
00159       const unsigned int n1_dofs = fe1.dofs_per_cell;
00160       const unsigned int n2_dofs = fe2.dofs_per_cell;
00161       const unsigned int n_components = fe1.get_fe().n_components();
00162 
00163       Assert(n1_dofs == n2_dofs, ExcNotImplemented());
00164       AssertDimension(n_components, fe2.get_fe().n_components());
00165       AssertDimension(M11.m(), n1_dofs);
00166       AssertDimension(M12.m(), n1_dofs);
00167       AssertDimension(M21.m(), n2_dofs);
00168       AssertDimension(M22.m(), n2_dofs);
00169       AssertDimension(M11.n(), n1_dofs);
00170       AssertDimension(M12.n(), n2_dofs);
00171       AssertDimension(M21.n(), n1_dofs);
00172       AssertDimension(M22.n(), n2_dofs);
00173 
00174       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00175         {
00176           const double dx = fe1.JxW(k);
00177 
00178           for (unsigned i=0;i<n1_dofs;++i)
00179             for (unsigned j=0;j<n1_dofs;++j)
00180               for (unsigned int d=0;d<n_components;++d)
00181                 {
00182                   const double u1 = factor1*fe1.shape_value_component(j,k,d);
00183                   const double u2 =-factor2*fe2.shape_value_component(j,k,d);
00184                   const double v1 = factor1*fe1.shape_value_component(i,k,d);
00185                   const double v2 =-factor2*fe2.shape_value_component(i,k,d);
00186 
00187                   M11(i,j) += dx * u1*v1;
00188                   M12(i,j) += dx * u2*v1;
00189                   M21(i,j) += dx * u1*v2;
00190                   M22(i,j) += dx * u2*v2;
00191                 }
00192         }
00193     }
00194   }
00195 }
00196 
00197 DEAL_II_NAMESPACE_CLOSE
00198 
00199 #endif
00200 
 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