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