Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
divergence.h
1 // ---------------------------------------------------------------------
2 // @f$Id: divergence.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__integrators_divergence_h
18 #define __deal2__integrators_divergence_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
28 
30 
31 namespace LocalIntegrators
32 {
40  namespace Divergence
41  {
52  template <int dim>
55  const Tensor<2,dim> &h0,
56  const Tensor<2,dim> &h1,
57  const Tensor<2,dim> &h2)
58  {
59  Tensor<1,dim> result;
60  for (unsigned int d=0; d<dim; ++d)
61  {
62  result[d] += h0[d][0];
63  if (dim >=2) result[d] += h1[d][1];
64  if (dim >=3) result[d] += h2[d][2];
65  }
66  return result;
67  }
68 
69 
83  template <int dim>
84  void cell_matrix (
86  const FEValuesBase<dim> &fe,
87  const FEValuesBase<dim> &fetest,
88  double factor = 1.)
89  {
90  unsigned int fecomp = fe.get_fe().n_components();
91  const unsigned int n_dofs = fe.dofs_per_cell;
92  const unsigned int t_dofs = fetest.dofs_per_cell;
93  AssertDimension(fecomp, dim);
94  AssertDimension(fetest.get_fe().n_components(), 1);
95  AssertDimension(M.m(), t_dofs);
96  AssertDimension(M.n(), n_dofs);
97 
98  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
99  {
100  const double dx = fe.JxW(k) * factor;
101  for (unsigned int i=0; i<t_dofs; ++i)
102  {
103  const double vv = fetest.shape_value(i,k);
104  for (unsigned int d=0; d<dim; ++d)
105  for (unsigned int j=0; j<n_dofs; ++j)
106  {
107  const double du = fe.shape_grad_component(j,k,d)[d];
108  M(i,j) += dx * du * vv;
109  }
110  }
111  }
112  }
113 
129  template <int dim, typename number>
131  Vector<number> &result,
132  const FEValuesBase<dim> &fetest,
133  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
134  const double factor = 1.)
135  {
136  AssertDimension(fetest.get_fe().n_components(), 1);
138  const unsigned int t_dofs = fetest.dofs_per_cell;
139  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
140 
141  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
142  {
143  const double dx = factor * fetest.JxW(k);
144 
145  for (unsigned int i=0; i<t_dofs; ++i)
146  for (unsigned int d=0; d<dim; ++d)
147  result(i) += dx * input[d][k][d] * fetest.shape_value(i,k);
148  }
149  }
150 
151 
167  template <int dim, typename number>
169  Vector<number> &result,
170  const FEValuesBase<dim> &fetest,
171  const VectorSlice<const std::vector<std::vector<double> > > &input,
172  const double factor = 1.)
173  {
174  AssertDimension(fetest.get_fe().n_components(), 1);
176  const unsigned int t_dofs = fetest.dofs_per_cell;
177  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
178 
179  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
180  {
181  const double dx = factor * fetest.JxW(k);
182 
183  for (unsigned int i=0; i<t_dofs; ++i)
184  for (unsigned int d=0; d<dim; ++d)
185  result(i) -= dx * input[d][k] * fetest.shape_grad(i,k)[d];
186  }
187  }
188 
189 
202  template <int dim>
205  const FEValuesBase<dim> &fe,
206  const FEValuesBase<dim> &fetest,
207  double factor = 1.)
208  {
209  unsigned int fecomp = fetest.get_fe().n_components();
210  const unsigned int t_dofs = fetest.dofs_per_cell;
211  const unsigned int n_dofs = fe.dofs_per_cell;
212 
213  AssertDimension(fecomp, dim);
214  AssertDimension(fe.get_fe().n_components(), 1);
215  AssertDimension(M.m(), t_dofs);
216  AssertDimension(M.n(), n_dofs);
217 
218  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
219  {
220  const double dx = fe.JxW(k) * factor;
221  for (unsigned int d=0; d<dim; ++d)
222  for (unsigned int i=0; i<t_dofs; ++i)
223  {
224  const double vv = fetest.shape_value_component(i,k,d);
225  for (unsigned int j=0; j<n_dofs; ++j)
226  {
227  const Tensor<1,dim> &Du = fe.shape_grad(j,k);
228  M(i,j) += dx * vv * Du[d];
229  }
230  }
231  }
232  }
233 
248  template <int dim, typename number>
250  Vector<number> &result,
251  const FEValuesBase<dim> &fetest,
252  const std::vector<Tensor<1,dim> > &input,
253  const double factor = 1.)
254  {
255  AssertDimension(fetest.get_fe().n_components(), dim);
256  AssertDimension(input.size(), fetest.n_quadrature_points);
257  const unsigned int t_dofs = fetest.dofs_per_cell;
258  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
259 
260  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
261  {
262  const double dx = factor * fetest.JxW(k);
263 
264  for (unsigned int i=0; i<t_dofs; ++i)
265  for (unsigned int d=0; d<dim; ++d)
266  result(i) += dx * input[k][d] * fetest.shape_value_component(i,k,d);
267  }
268  }
269 
284  template <int dim, typename number>
286  Vector<number> &result,
287  const FEValuesBase<dim> &fetest,
288  const std::vector<double> &input,
289  const double factor = 1.)
290  {
291  AssertDimension(fetest.get_fe().n_components(), dim);
292  AssertDimension(input.size(), fetest.n_quadrature_points);
293  const unsigned int t_dofs = fetest.dofs_per_cell;
294  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
295 
296  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
297  {
298  const double dx = factor * fetest.JxW(k);
299 
300  for (unsigned int i=0; i<t_dofs; ++i)
301  for (unsigned int d=0; d<dim; ++d)
302  result(i) -= dx * input[k] * fetest.shape_grad_component(i,k,d)[d];
303  }
304  }
305 
315  template<int dim>
316  void
319  const FEValuesBase<dim> &fe,
320  const FEValuesBase<dim> &fetest,
321  double factor = 1.)
322  {
323  unsigned int fecomp = fe.get_fe().n_components();
324  const unsigned int n_dofs = fe.dofs_per_cell;
325  const unsigned int t_dofs = fetest.dofs_per_cell;
326 
327  AssertDimension(fecomp, dim);
328  AssertDimension(fetest.get_fe().n_components(), 1);
329  AssertDimension(M.m(), t_dofs);
330  AssertDimension(M.n(), n_dofs);
331 
332  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
333  {
334  const Tensor<1,dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
335  for (unsigned int i=0; i<t_dofs; ++i)
336  for (unsigned int j=0; j<n_dofs; ++j)
337  for (unsigned int d=0; d<dim; ++d)
338  M(i,j) += ndx[d] * fe.shape_value_component(j,k,d)
339  * fetest.shape_value(i,k);
340  }
341  }
342 
356  template<int dim, typename number>
357  void
359  Vector<number> &result,
360  const FEValuesBase<dim> &fe,
361  const FEValuesBase<dim> &fetest,
362  const VectorSlice<const std::vector<std::vector<double> > > &data,
363  double factor = 1.)
364  {
365  unsigned int fecomp = fe.get_fe().n_components();
366  const unsigned int t_dofs = fetest.dofs_per_cell;
367 
368  AssertDimension(fecomp, dim);
369  AssertDimension(fetest.get_fe().n_components(), 1);
370  AssertDimension(result.size(), t_dofs);
372 
373  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
374  {
375  const Tensor<1,dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
376 
377  for (unsigned int i=0; i<t_dofs; ++i)
378  for (unsigned int d=0; d<dim; ++d)
379  result(i) += ndx[d] * fetest.shape_value(i,k) * data[d][k];
380  }
381  }
382 
396  template<int dim, typename number>
397  void
399  Vector<number> &result,
400  const FEValuesBase<dim> &fetest,
401  const std::vector<double> &data,
402  double factor = 1.)
403  {
404  const unsigned int t_dofs = fetest.dofs_per_cell;
405 
406  AssertDimension(fetest.get_fe().n_components(), dim);
407  AssertDimension(result.size(), t_dofs);
408  AssertDimension(data.size(), fetest.n_quadrature_points);
409 
410  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
411  {
412  const Tensor<1,dim> ndx = factor * fetest.normal_vector(k) * fetest.JxW(k);
413 
414  for (unsigned int i=0; i<t_dofs; ++i)
415  for (unsigned int d=0; d<dim; ++d)
416  result(i) += ndx[d] * fetest.shape_value_component(i,k,d) * data[k];
417  }
418  }
419 
433  template<int dim>
434  void
436  FullMatrix<double> &M11,
437  FullMatrix<double> &M12,
438  FullMatrix<double> &M21,
439  FullMatrix<double> &M22,
440  const FEValuesBase<dim> &fe1,
441  const FEValuesBase<dim> &fe2,
442  const FEValuesBase<dim> &fetest1,
443  const FEValuesBase<dim> &fetest2,
444  double factor = 1.)
445  {
446  const unsigned int n_dofs = fe1.dofs_per_cell;
447  const unsigned int t_dofs = fetest1.dofs_per_cell;
448 
449  AssertDimension(fe1.get_fe().n_components(), dim);
450  AssertDimension(fe2.get_fe().n_components(), dim);
451  AssertDimension(fetest1.get_fe().n_components(), 1);
452  AssertDimension(fetest2.get_fe().n_components(), 1);
453  AssertDimension(M11.m(), t_dofs);
454  AssertDimension(M11.n(), n_dofs);
455  AssertDimension(M12.m(), t_dofs);
456  AssertDimension(M12.n(), n_dofs);
457  AssertDimension(M21.m(), t_dofs);
458  AssertDimension(M21.n(), n_dofs);
459  AssertDimension(M22.m(), t_dofs);
460  AssertDimension(M22.n(), n_dofs);
461 
462  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
463  {
464  const double dx = factor * fe1.JxW(k);
465  for (unsigned int i=0; i<t_dofs; ++i)
466  for (unsigned int j=0; j<n_dofs; ++j)
467  for (unsigned int d=0; d<dim; ++d)
468  {
469  const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
470  const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
471  const double v1 = fetest1.shape_value(i,k);
472  const double v2 = fetest2.shape_value(i,k);
473 
474  M11(i,j) += .5 * dx * un1 * v1;
475  M12(i,j) += .5 * dx * un2 * v1;
476  M21(i,j) += .5 * dx * un1 * v2;
477  M22(i,j) += .5 * dx * un2 * v2;
478  }
479  }
480  }
481 
491  template <int dim>
494  const FEValuesBase<dim> &fe,
495  double factor = 1.)
496  {
497  const unsigned int n_dofs = fe.dofs_per_cell;
498 
499  AssertDimension(fe.get_fe().n_components(), dim);
500  AssertDimension(M.m(), n_dofs);
501  AssertDimension(M.n(), n_dofs);
502 
503  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
504  {
505  const double dx = factor * fe.JxW(k);
506  for (unsigned int i=0; i<n_dofs; ++i)
507  for (unsigned int j=0; j<n_dofs; ++j)
508  {
509  double dv = 0.;
510  double du = 0.;
511  for (unsigned int d=0; d<dim; ++d)
512  {
513  dv += fe.shape_grad_component(i,k,d)[d];
514  du += fe.shape_grad_component(j,k,d)[d];
515  }
516 
517  M(i,j) += dx * du * dv;
518  }
519  }
520  }
521 
534  template<int dim>
535  void
537  FullMatrix<double> &M11,
538  FullMatrix<double> &M12,
539  FullMatrix<double> &M21,
540  FullMatrix<double> &M22,
541  const FEValuesBase<dim> &fe1,
542  const FEValuesBase<dim> &fe2,
543  double factor = 1.)
544  {
545  const unsigned int n_dofs = fe1.dofs_per_cell;
546 
547  AssertDimension(fe1.get_fe().n_components(), dim);
548  AssertDimension(fe2.get_fe().n_components(), dim);
549  AssertDimension(M11.m(), n_dofs);
550  AssertDimension(M11.n(), n_dofs);
551  AssertDimension(M12.m(), n_dofs);
552  AssertDimension(M12.n(), n_dofs);
553  AssertDimension(M21.m(), n_dofs);
554  AssertDimension(M21.n(), n_dofs);
555  AssertDimension(M22.m(), n_dofs);
556  AssertDimension(M22.n(), n_dofs);
557 
558  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
559  {
560  const double dx = factor * fe1.JxW(k);
561  for (unsigned int i=0; i<n_dofs; ++i)
562  for (unsigned int j=0; j<n_dofs; ++j)
563  for (unsigned int d=0; d<dim; ++d)
564  {
565  const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
566  const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
567  const double vn1 = fe1.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
568  const double vn2 =-fe2.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
569 
570  M11(i,j) += dx * un1 * vn1;
571  M12(i,j) += dx * un2 * vn1;
572  M21(i,j) += dx * un1 * vn2;
573  M22(i,j) += dx * un2 * vn2;
574  }
575  }
576  }
577 
590  template <int dim>
591  double norm(const FEValuesBase<dim> &fe,
592  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Du)
593  {
594  unsigned int fecomp = fe.get_fe().n_components();
595 
596  AssertDimension(fecomp, dim);
598 
599  double result = 0;
600  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
601  {
602  double div = Du[0][k][0];
603  for (unsigned int d=1; d<dim; ++d)
604  div += Du[d][k][d];
605  result += div*div*fe.JxW(k);
606  }
607  return result;
608  }
609 
610  }
611 }
612 
613 
614 DEAL_II_NAMESPACE_CLOSE
615 
616 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
size_type m() const
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:870
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &input, const double factor=1.)
Definition: divergence.h:130
size_type n() const
const Point< spacedim > & normal_vector(const unsigned int i) const
Library of integrals over cells and faces.
Definition: elasticity.h:31
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:84
double JxW(const unsigned int quadrature_point) const
#define Assert(cond, exc)
Definition: exceptions.h:299
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:203
void u_dot_n_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:317
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void u_dot_n_jump_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, double factor=1.)
Definition: divergence.h:536
void u_times_n_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< double > &data, double factor=1.)
Definition: divergence.h:398
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
Definition: divergence.h:591
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void u_dot_n_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const VectorSlice< const std::vector< std::vector< double > > > &data, double factor=1.)
Definition: divergence.h:358
const FiniteElement< dim, spacedim > & get_fe() const
void gradient_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< Tensor< 1, dim > > &input, const double factor=1.)
Definition: divergence.h:249
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
Tensor< 1, dim > grad_div(const Tensor< 2, dim > &h0, const Tensor< 2, dim > &h1, const Tensor< 2, dim > &h2)
Definition: divergence.h:54
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
void grad_div_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: divergence.h:492
std::size_t size() const