00001
00002
00003
00004
00005
00006
00007
00008
00009
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