Reference documentation for deal.II version 9.0.0
divergence.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2017 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_integrators_divergence_h
17 #define dealii_integrators_divergence_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/quadrature.h>
23 #include <deal.II/lac/full_matrix.h>
24 #include <deal.II/fe/mapping.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/meshworker/dof_info.h>
27 #include <deal.II/integrators/grad_div.h>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 namespace LocalIntegrators
32 {
41  namespace Divergence
42  {
52  template <int dim>
53  void cell_matrix (
55  const FEValuesBase<dim> &fe,
56  const FEValuesBase<dim> &fetest,
57  double factor = 1.)
58  {
59  const unsigned int n_dofs = fe.dofs_per_cell;
60  const unsigned int t_dofs = fetest.dofs_per_cell;
61  AssertDimension(fe.get_fe().n_components(), dim);
62  AssertDimension(fetest.get_fe().n_components(), 1);
63  AssertDimension(M.m(), t_dofs);
64  AssertDimension(M.n(), n_dofs);
65 
66  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
67  {
68  const double dx = fe.JxW(k) * factor;
69  for (unsigned int i=0; i<t_dofs; ++i)
70  {
71  const double vv = fetest.shape_value(i,k);
72  for (unsigned int d=0; d<dim; ++d)
73  for (unsigned int j=0; j<n_dofs; ++j)
74  {
75  const double du = fe.shape_grad_component(j,k,d)[d];
76  M(i,j) += dx * du * vv;
77  }
78  }
79  }
80  }
81 
94  template <int dim, typename number>
96  Vector<number> &result,
97  const FEValuesBase<dim> &fetest,
98  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
99  const double factor = 1.)
100  {
101  AssertDimension(fetest.get_fe().n_components(), 1);
103  const unsigned int t_dofs = fetest.dofs_per_cell;
104  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
105 
106  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
107  {
108  const double dx = factor * fetest.JxW(k);
109 
110  for (unsigned int i=0; i<t_dofs; ++i)
111  for (unsigned int d=0; d<dim; ++d)
112  result(i) += dx * input[d][k][d] * fetest.shape_value(i,k);
113  }
114  }
115 
116 
129  template <int dim, typename number>
131  Vector<number> &result,
132  const FEValuesBase<dim> &fetest,
133  const VectorSlice<const std::vector<std::vector<double> > > &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] * fetest.shape_grad(i,k)[d];
148  }
149  }
150 
151 
162  template <int dim>
165  const FEValuesBase<dim> &fe,
166  const FEValuesBase<dim> &fetest,
167  double factor = 1.)
168  {
169  const unsigned int t_dofs = fetest.dofs_per_cell;
170  const unsigned int n_dofs = fe.dofs_per_cell;
171 
172  AssertDimension(fetest.get_fe().n_components(), dim);
173  AssertDimension(fe.get_fe().n_components(), 1);
174  AssertDimension(M.m(), t_dofs);
175  AssertDimension(M.n(), n_dofs);
176 
177  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
178  {
179  const double dx = fe.JxW(k) * factor;
180  for (unsigned int d=0; d<dim; ++d)
181  for (unsigned int i=0; i<t_dofs; ++i)
182  {
183  const double vv = fetest.shape_value_component(i,k,d);
184  for (unsigned int j=0; j<n_dofs; ++j)
185  {
186  const Tensor<1,dim> &Du = fe.shape_grad(j,k);
187  M(i,j) += dx * vv * Du[d];
188  }
189  }
190  }
191  }
192 
205  template <int dim, typename number>
207  Vector<number> &result,
208  const FEValuesBase<dim> &fetest,
209  const std::vector<Tensor<1,dim> > &input,
210  const double factor = 1.)
211  {
212  AssertDimension(fetest.get_fe().n_components(), dim);
213  AssertDimension(input.size(), fetest.n_quadrature_points);
214  const unsigned int t_dofs = fetest.dofs_per_cell;
215  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
216 
217  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
218  {
219  const double dx = factor * fetest.JxW(k);
220 
221  for (unsigned int i=0; i<t_dofs; ++i)
222  for (unsigned int d=0; d<dim; ++d)
223  result(i) += dx * input[k][d] * fetest.shape_value_component(i,k,d);
224  }
225  }
226 
239  template <int dim, typename number>
241  Vector<number> &result,
242  const FEValuesBase<dim> &fetest,
243  const std::vector<double> &input,
244  const double factor = 1.)
245  {
246  AssertDimension(fetest.get_fe().n_components(), dim);
247  AssertDimension(input.size(), fetest.n_quadrature_points);
248  const unsigned int t_dofs = fetest.dofs_per_cell;
249  Assert (result.size() == t_dofs, ExcDimensionMismatch(result.size(), t_dofs));
250 
251  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
252  {
253  const double dx = factor * fetest.JxW(k);
254 
255  for (unsigned int i=0; i<t_dofs; ++i)
256  for (unsigned int d=0; d<dim; ++d)
257  result(i) -= dx * input[k] * fetest.shape_grad_component(i,k,d)[d];
258  }
259  }
260 
269  template <int dim>
270  void
273  const FEValuesBase<dim> &fe,
274  const FEValuesBase<dim> &fetest,
275  double factor = 1.)
276  {
277  const unsigned int n_dofs = fe.dofs_per_cell;
278  const unsigned int t_dofs = fetest.dofs_per_cell;
279 
280  AssertDimension(fe.get_fe().n_components(), dim);
281  AssertDimension(fetest.get_fe().n_components(), 1);
282  AssertDimension(M.m(), t_dofs);
283  AssertDimension(M.n(), n_dofs);
284 
285  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
286  {
287  const Tensor<1,dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
288  for (unsigned int i=0; i<t_dofs; ++i)
289  for (unsigned int j=0; j<n_dofs; ++j)
290  for (unsigned int d=0; d<dim; ++d)
291  M(i,j) += ndx[d] * fe.shape_value_component(j,k,d)
292  * fetest.shape_value(i,k);
293  }
294  }
295 
306  template <int dim, typename number>
307  void
309  Vector<number> &result,
310  const FEValuesBase<dim> &fe,
311  const FEValuesBase<dim> &fetest,
312  const VectorSlice<const std::vector<std::vector<double> > > &data,
313  double factor = 1.)
314  {
315  const unsigned int t_dofs = fetest.dofs_per_cell;
316 
317  AssertDimension(fe.get_fe().n_components(), dim);
318  AssertDimension(fetest.get_fe().n_components(), 1);
319  AssertDimension(result.size(), t_dofs);
321 
322  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
323  {
324  const Tensor<1,dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
325 
326  for (unsigned int i=0; i<t_dofs; ++i)
327  for (unsigned int d=0; d<dim; ++d)
328  result(i) += ndx[d] * fetest.shape_value(i,k) * data[d][k];
329  }
330  }
331 
342  template <int dim, typename number>
343  void
345  Vector<number> &result,
346  const FEValuesBase<dim> &fetest,
347  const std::vector<double> &data,
348  double factor = 1.)
349  {
350  const unsigned int t_dofs = fetest.dofs_per_cell;
351 
352  AssertDimension(fetest.get_fe().n_components(), dim);
353  AssertDimension(result.size(), t_dofs);
354  AssertDimension(data.size(), fetest.n_quadrature_points);
355 
356  for (unsigned int k=0; k<fetest.n_quadrature_points; ++k)
357  {
358  const Tensor<1,dim> ndx = factor * fetest.normal_vector(k) * fetest.JxW(k);
359 
360  for (unsigned int i=0; i<t_dofs; ++i)
361  for (unsigned int d=0; d<dim; ++d)
362  result(i) += ndx[d] * fetest.shape_value_component(i,k,d) * data[k];
363  }
364  }
365 
377  template <int dim>
378  void
380  FullMatrix<double> &M11,
381  FullMatrix<double> &M12,
382  FullMatrix<double> &M21,
383  FullMatrix<double> &M22,
384  const FEValuesBase<dim> &fe1,
385  const FEValuesBase<dim> &fe2,
386  const FEValuesBase<dim> &fetest1,
387  const FEValuesBase<dim> &fetest2,
388  double factor = 1.)
389  {
390  const unsigned int n_dofs = fe1.dofs_per_cell;
391  const unsigned int t_dofs = fetest1.dofs_per_cell;
392 
393  AssertDimension(fe1.get_fe().n_components(), dim);
394  AssertDimension(fe2.get_fe().n_components(), dim);
395  AssertDimension(fetest1.get_fe().n_components(), 1);
396  AssertDimension(fetest2.get_fe().n_components(), 1);
397  AssertDimension(M11.m(), t_dofs);
398  AssertDimension(M11.n(), n_dofs);
399  AssertDimension(M12.m(), t_dofs);
400  AssertDimension(M12.n(), n_dofs);
401  AssertDimension(M21.m(), t_dofs);
402  AssertDimension(M21.n(), n_dofs);
403  AssertDimension(M22.m(), t_dofs);
404  AssertDimension(M22.n(), n_dofs);
405 
406  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
407  {
408  const double dx = factor * fe1.JxW(k);
409  for (unsigned int i=0; i<t_dofs; ++i)
410  for (unsigned int j=0; j<n_dofs; ++j)
411  for (unsigned int d=0; d<dim; ++d)
412  {
413  const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
414  const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
415  const double v1 = fetest1.shape_value(i,k);
416  const double v2 = fetest2.shape_value(i,k);
417 
418  M11(i,j) += .5 * dx * un1 * v1;
419  M12(i,j) += .5 * dx * un2 * v1;
420  M21(i,j) += .5 * dx * un1 * v2;
421  M22(i,j) += .5 * dx * un2 * v2;
422  }
423  }
424  }
425 
429  template <int dim>
430  DEAL_II_DEPRECATED
431  void grad_div_matrix (
433  const FEValuesBase<dim> &fe,
434  const double factor = 1.);
435 
436  template <int dim>
439  const FEValuesBase<dim> &fe,
440  const double factor)
441  {
442  GradDiv::cell_matrix(M, fe, factor);
443  }
444 
448  template <int dim, typename number>
449  DEAL_II_DEPRECATED
450  void grad_div_residual (
451  Vector<number> &result,
452  const FEValuesBase<dim> &fetest,
453  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
454  const double factor = 1.);
455 
456  template <int dim, typename number>
458  Vector<number> &result,
459  const FEValuesBase<dim> &fetest,
460  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
461  const double factor)
462  {
463  GradDiv::cell_residual(result, fetest, input, factor);
464  }
465 
478  template <int dim>
479  void
481  FullMatrix<double> &M11,
482  FullMatrix<double> &M12,
483  FullMatrix<double> &M21,
484  FullMatrix<double> &M22,
485  const FEValuesBase<dim> &fe1,
486  const FEValuesBase<dim> &fe2,
487  double factor = 1.)
488  {
489  const unsigned int n_dofs = fe1.dofs_per_cell;
490 
491  AssertDimension(fe1.get_fe().n_components(), dim);
492  AssertDimension(fe2.get_fe().n_components(), dim);
493  AssertDimension(M11.m(), n_dofs);
494  AssertDimension(M11.n(), n_dofs);
495  AssertDimension(M12.m(), n_dofs);
496  AssertDimension(M12.n(), n_dofs);
497  AssertDimension(M21.m(), n_dofs);
498  AssertDimension(M21.n(), n_dofs);
499  AssertDimension(M22.m(), n_dofs);
500  AssertDimension(M22.n(), n_dofs);
501 
502  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
503  {
504  const double dx = factor * fe1.JxW(k);
505  for (unsigned int i=0; i<n_dofs; ++i)
506  for (unsigned int j=0; j<n_dofs; ++j)
507  for (unsigned int d=0; d<dim; ++d)
508  {
509  const double un1 = fe1.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
510  const double un2 =-fe2.shape_value_component(j,k,d) * fe1.normal_vector(k)[d];
511  const double vn1 = fe1.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
512  const double vn2 =-fe2.shape_value_component(i,k,d) * fe1.normal_vector(k)[d];
513 
514  M11(i,j) += dx * un1 * vn1;
515  M12(i,j) += dx * un2 * vn1;
516  M21(i,j) += dx * un1 * vn2;
517  M22(i,j) += dx * un2 * vn2;
518  }
519  }
520  }
521 
533  template <int dim>
534  double norm(const FEValuesBase<dim> &fe,
535  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Du)
536  {
537  AssertDimension(fe.get_fe().n_components(), dim);
539 
540  double result = 0;
541  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
542  {
543  double div = Du[0][k][0];
544  for (unsigned int d=1; d<dim; ++d)
545  div += Du[d][k][d];
546  result += div*div*fe.JxW(k);
547  }
548  return result;
549  }
550 
551  }
552 }
553 
554 
555 DEAL_II_NAMESPACE_CLOSE
556 
557 #endif
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1248
const unsigned int dofs_per_cell
Definition: fe_values.h:1847
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: grad_div.h:52
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:1259
const FiniteElement< dim, spacedim > & get_fe() const
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:95
void grad_div_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:457
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
Library of integrals over cells and faces.
Definition: advection.h:30
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:53
size_type n() const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
#define Assert(cond, exc)
Definition: exceptions.h:1142
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:163
void u_dot_n_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:271
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:480
const unsigned int n_quadrature_points
Definition: fe_values.h:1840
void u_times_n_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< double > &data, double factor=1.)
Definition: divergence.h:344
void grad_div_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: divergence.h:437
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
Definition: divergence.h:534
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:308
void gradient_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< Tensor< 1, dim > > &input, const double factor=1.)
Definition: divergence.h:206
double JxW(const unsigned int quadrature_point) const
size_type size() const
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: grad_div.h:87
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const