Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
divergence.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 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.md at
12 // the top level directory of deal.II.
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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/quadrature.h>
24 
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
28 #include <deal.II/integrators/grad_div.h>
29 
30 #include <deal.II/lac/full_matrix.h>
31 
32 #include <deal.II/meshworker/dof_info.h>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 namespace LocalIntegrators
37 {
46  namespace Divergence
47  {
57  template <int dim>
58  void
59  cell_matrix(FullMatrix<double> & M,
60  const FEValuesBase<dim> &fe,
61  const FEValuesBase<dim> &fetest,
62  double factor = 1.)
63  {
64  const unsigned int n_dofs = fe.dofs_per_cell;
65  const unsigned int t_dofs = fetest.dofs_per_cell;
66  AssertDimension(fe.get_fe().n_components(), dim);
67  AssertDimension(fetest.get_fe().n_components(), 1);
68  AssertDimension(M.m(), t_dofs);
69  AssertDimension(M.n(), n_dofs);
70 
71  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
72  {
73  const double dx = fe.JxW(k) * factor;
74  for (unsigned int i = 0; i < t_dofs; ++i)
75  {
76  const double vv = fetest.shape_value(i, k);
77  for (unsigned int d = 0; d < dim; ++d)
78  for (unsigned int j = 0; j < n_dofs; ++j)
79  {
80  const double du = fe.shape_grad_component(j, k, d)[d];
81  M(i, j) += dx * du * vv;
82  }
83  }
84  }
85  }
86 
99  template <int dim, typename number>
100  void
101  cell_residual(Vector<number> & result,
102  const FEValuesBase<dim> & fetest,
103  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
104  const double factor = 1.)
105  {
106  AssertDimension(fetest.get_fe().n_components(), 1);
107  AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
108  const unsigned int t_dofs = fetest.dofs_per_cell;
109  Assert(result.size() == t_dofs,
110  ExcDimensionMismatch(result.size(), t_dofs));
111 
112  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
113  {
114  const double dx = factor * fetest.JxW(k);
115 
116  for (unsigned int i = 0; i < t_dofs; ++i)
117  for (unsigned int d = 0; d < dim; ++d)
118  result(i) += dx * input[d][k][d] * fetest.shape_value(i, k);
119  }
120  }
121 
122 
135  template <int dim, typename number>
136  void
137  cell_residual(Vector<number> & result,
138  const FEValuesBase<dim> & fetest,
139  const ArrayView<const std::vector<double>> &input,
140  const double factor = 1.)
141  {
142  AssertDimension(fetest.get_fe().n_components(), 1);
143  AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
144  const unsigned int t_dofs = fetest.dofs_per_cell;
145  Assert(result.size() == t_dofs,
146  ExcDimensionMismatch(result.size(), t_dofs));
147 
148  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
149  {
150  const double dx = factor * fetest.JxW(k);
151 
152  for (unsigned int i = 0; i < t_dofs; ++i)
153  for (unsigned int d = 0; d < dim; ++d)
154  result(i) -= dx * input[d][k] * fetest.shape_grad(i, k)[d];
155  }
156  }
157 
158 
169  template <int dim>
170  void
171  gradient_matrix(FullMatrix<double> & M,
172  const FEValuesBase<dim> &fe,
173  const FEValuesBase<dim> &fetest,
174  double factor = 1.)
175  {
176  const unsigned int t_dofs = fetest.dofs_per_cell;
177  const unsigned int n_dofs = fe.dofs_per_cell;
178 
179  AssertDimension(fetest.get_fe().n_components(), dim);
180  AssertDimension(fe.get_fe().n_components(), 1);
181  AssertDimension(M.m(), t_dofs);
182  AssertDimension(M.n(), n_dofs);
183 
184  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
185  {
186  const double dx = fe.JxW(k) * factor;
187  for (unsigned int d = 0; d < dim; ++d)
188  for (unsigned int i = 0; i < t_dofs; ++i)
189  {
190  const double vv = fetest.shape_value_component(i, k, d);
191  for (unsigned int j = 0; j < n_dofs; ++j)
192  {
193  const Tensor<1, dim> &Du = fe.shape_grad(j, k);
194  M(i, j) += dx * vv * Du[d];
195  }
196  }
197  }
198  }
199 
212  template <int dim, typename number>
213  void
214  gradient_residual(Vector<number> & result,
215  const FEValuesBase<dim> & fetest,
216  const std::vector<Tensor<1, dim>> &input,
217  const double factor = 1.)
218  {
219  AssertDimension(fetest.get_fe().n_components(), dim);
220  AssertDimension(input.size(), fetest.n_quadrature_points);
221  const unsigned int t_dofs = fetest.dofs_per_cell;
222  Assert(result.size() == t_dofs,
223  ExcDimensionMismatch(result.size(), t_dofs));
224 
225  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
226  {
227  const double dx = factor * fetest.JxW(k);
228 
229  for (unsigned int i = 0; i < t_dofs; ++i)
230  for (unsigned int d = 0; d < dim; ++d)
231  result(i) +=
232  dx * input[k][d] * fetest.shape_value_component(i, k, d);
233  }
234  }
235 
248  template <int dim, typename number>
249  void
250  gradient_residual(Vector<number> & result,
251  const FEValuesBase<dim> & fetest,
252  const std::vector<double> &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,
259  ExcDimensionMismatch(result.size(), t_dofs));
260 
261  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
262  {
263  const double dx = factor * fetest.JxW(k);
264 
265  for (unsigned int i = 0; i < t_dofs; ++i)
266  for (unsigned int d = 0; d < dim; ++d)
267  result(i) -=
268  dx * input[k] * fetest.shape_grad_component(i, k, d)[d];
269  }
270  }
271 
280  template <int dim>
281  void
282  u_dot_n_matrix(FullMatrix<double> & M,
283  const FEValuesBase<dim> &fe,
284  const FEValuesBase<dim> &fetest,
285  double factor = 1.)
286  {
287  const unsigned int n_dofs = fe.dofs_per_cell;
288  const unsigned int t_dofs = fetest.dofs_per_cell;
289 
290  AssertDimension(fe.get_fe().n_components(), dim);
291  AssertDimension(fetest.get_fe().n_components(), 1);
292  AssertDimension(M.m(), t_dofs);
293  AssertDimension(M.n(), n_dofs);
294 
295  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
296  {
297  const Tensor<1, dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
298  for (unsigned int i = 0; i < t_dofs; ++i)
299  for (unsigned int j = 0; j < n_dofs; ++j)
300  for (unsigned int d = 0; d < dim; ++d)
301  M(i, j) += ndx[d] * fe.shape_value_component(j, k, d) *
302  fetest.shape_value(i, k);
303  }
304  }
305 
316  template <int dim, typename number>
317  void
318  u_dot_n_residual(Vector<number> & result,
319  const FEValuesBase<dim> & fe,
320  const FEValuesBase<dim> & fetest,
321  const ArrayView<const std::vector<double>> &data,
322  double factor = 1.)
323  {
324  const unsigned int t_dofs = fetest.dofs_per_cell;
325 
326  AssertDimension(fe.get_fe().n_components(), dim);
327  AssertDimension(fetest.get_fe().n_components(), 1);
328  AssertDimension(result.size(), t_dofs);
329  AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
330 
331  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
332  {
333  const Tensor<1, dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
334 
335  for (unsigned int i = 0; i < t_dofs; ++i)
336  for (unsigned int d = 0; d < dim; ++d)
337  result(i) += ndx[d] * fetest.shape_value(i, k) * data[d][k];
338  }
339  }
340 
351  template <int dim, typename number>
352  void
353  u_times_n_residual(Vector<number> & result,
354  const FEValuesBase<dim> & fetest,
355  const std::vector<double> &data,
356  double factor = 1.)
357  {
358  const unsigned int t_dofs = fetest.dofs_per_cell;
359 
360  AssertDimension(fetest.get_fe().n_components(), dim);
361  AssertDimension(result.size(), t_dofs);
362  AssertDimension(data.size(), fetest.n_quadrature_points);
363 
364  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
365  {
366  const Tensor<1, dim> ndx =
367  factor * fetest.normal_vector(k) * fetest.JxW(k);
368 
369  for (unsigned int i = 0; i < t_dofs; ++i)
370  for (unsigned int d = 0; d < dim; ++d)
371  result(i) +=
372  ndx[d] * fetest.shape_value_component(i, k, d) * data[k];
373  }
374  }
375 
388  template <int dim>
389  void
390  u_dot_n_matrix(FullMatrix<double> & M11,
391  FullMatrix<double> & M12,
392  FullMatrix<double> & M21,
393  FullMatrix<double> & M22,
394  const FEValuesBase<dim> &fe1,
395  const FEValuesBase<dim> &fe2,
396  const FEValuesBase<dim> &fetest1,
397  const FEValuesBase<dim> &fetest2,
398  double factor = 1.)
399  {
400  const unsigned int n_dofs = fe1.dofs_per_cell;
401  const unsigned int t_dofs = fetest1.dofs_per_cell;
402 
403  AssertDimension(fe1.get_fe().n_components(), dim);
404  AssertDimension(fe2.get_fe().n_components(), dim);
405  AssertDimension(fetest1.get_fe().n_components(), 1);
406  AssertDimension(fetest2.get_fe().n_components(), 1);
407  AssertDimension(M11.m(), t_dofs);
408  AssertDimension(M11.n(), n_dofs);
409  AssertDimension(M12.m(), t_dofs);
410  AssertDimension(M12.n(), n_dofs);
411  AssertDimension(M21.m(), t_dofs);
412  AssertDimension(M21.n(), n_dofs);
413  AssertDimension(M22.m(), t_dofs);
414  AssertDimension(M22.n(), n_dofs);
415 
416  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
417  {
418  const double dx = factor * fe1.JxW(k);
419  for (unsigned int i = 0; i < t_dofs; ++i)
420  for (unsigned int j = 0; j < n_dofs; ++j)
421  for (unsigned int d = 0; d < dim; ++d)
422  {
423  const double un1 = fe1.shape_value_component(j, k, d) *
424  fe1.normal_vector(k)[d];
425  const double un2 = -fe2.shape_value_component(j, k, d) *
426  fe1.normal_vector(k)[d];
427  const double v1 = fetest1.shape_value(i, k);
428  const double v2 = fetest2.shape_value(i, k);
429 
430  M11(i, j) += .5 * dx * un1 * v1;
431  M12(i, j) += .5 * dx * un2 * v1;
432  M21(i, j) += .5 * dx * un1 * v2;
433  M22(i, j) += .5 * dx * un2 * v2;
434  }
435  }
436  }
437 
441  template <int dim>
442  DEAL_II_DEPRECATED void
443  grad_div_matrix(FullMatrix<double> & M,
444  const FEValuesBase<dim> &fe,
445  const double factor = 1.);
446 
447  template <int dim>
448  void
449  grad_div_matrix(FullMatrix<double> & M,
450  const FEValuesBase<dim> &fe,
451  const double factor)
452  {
453  GradDiv::cell_matrix(M, fe, factor);
454  }
455 
459  template <int dim, typename number>
460  DEAL_II_DEPRECATED void
461  grad_div_residual(Vector<number> & result,
462  const FEValuesBase<dim> &fetest,
463  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
464  const double factor = 1.);
465 
466  template <int dim, typename number>
467  void
468  grad_div_residual(Vector<number> & result,
469  const FEValuesBase<dim> &fetest,
470  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
471  const double factor)
472  {
473  GradDiv::cell_residual(result, fetest, input, factor);
474  }
475 
488  template <int dim>
489  void
490  u_dot_n_jump_matrix(FullMatrix<double> & M11,
491  FullMatrix<double> & M12,
492  FullMatrix<double> & M21,
493  FullMatrix<double> & M22,
494  const FEValuesBase<dim> &fe1,
495  const FEValuesBase<dim> &fe2,
496  double factor = 1.)
497  {
498  const unsigned int n_dofs = fe1.dofs_per_cell;
499 
500  AssertDimension(fe1.get_fe().n_components(), dim);
501  AssertDimension(fe2.get_fe().n_components(), dim);
502  AssertDimension(M11.m(), n_dofs);
503  AssertDimension(M11.n(), n_dofs);
504  AssertDimension(M12.m(), n_dofs);
505  AssertDimension(M12.n(), n_dofs);
506  AssertDimension(M21.m(), n_dofs);
507  AssertDimension(M21.n(), n_dofs);
508  AssertDimension(M22.m(), n_dofs);
509  AssertDimension(M22.n(), n_dofs);
510 
511  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
512  {
513  const double dx = factor * fe1.JxW(k);
514  for (unsigned int i = 0; i < n_dofs; ++i)
515  for (unsigned int j = 0; j < n_dofs; ++j)
516  for (unsigned int d = 0; d < dim; ++d)
517  {
518  const double un1 = fe1.shape_value_component(j, k, d) *
519  fe1.normal_vector(k)[d];
520  const double un2 = -fe2.shape_value_component(j, k, d) *
521  fe1.normal_vector(k)[d];
522  const double vn1 = fe1.shape_value_component(i, k, d) *
523  fe1.normal_vector(k)[d];
524  const double vn2 = -fe2.shape_value_component(i, k, d) *
525  fe1.normal_vector(k)[d];
526 
527  M11(i, j) += dx * un1 * vn1;
528  M12(i, j) += dx * un2 * vn1;
529  M21(i, j) += dx * un1 * vn2;
530  M22(i, j) += dx * un2 * vn2;
531  }
532  }
533  }
534 
546  template <int dim>
547  double
548  norm(const FEValuesBase<dim> & fe,
549  const ArrayView<const std::vector<Tensor<1, dim>>> &Du)
550  {
551  AssertDimension(fe.get_fe().n_components(), dim);
552  AssertVectorVectorDimension(Du, dim, fe.n_quadrature_points);
553 
554  double result = 0;
555  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
556  {
557  double div = Du[0][k][0];
558  for (unsigned int d = 1; d < dim; ++d)
559  div += Du[d][k][d];
560  result += div * div * fe.JxW(k);
561  }
562  return result;
563  }
564 
565  } // namespace Divergence
566 } // namespace LocalIntegrators
567 
568 
569 DEAL_II_NAMESPACE_CLOSE
570 
571 #endif
FEValuesBase::shape_value
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
LocalIntegrators::Divergence::grad_div_residual
void grad_div_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, const double factor=1.)
Definition: divergence.h:468
fe_values.h
LocalIntegrators::Divergence::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:59
LocalIntegrators::GradDiv::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: grad_div.h:57
LocalIntegrators::Divergence::u_dot_n_residual
void u_dot_n_residual(Vector< number > &result, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &data, double factor=1.)
Definition: divergence.h:318
FEValuesBase::shape_value_component
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
ArrayView
Definition: array_view.h:77
LocalIntegrators::Divergence::u_dot_n_matrix
void u_dot_n_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:282
FullMatrix::m
size_type m() const
FEValuesBase::normal_vector
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
DataOutBase::dx
@ dx
Definition: data_out_base.h:1562
FullMatrix::n
size_type n() const
Vector::size
size_type size() const
LocalIntegrators::Divergence::gradient_matrix
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:171
LocalIntegrators::Divergence::u_dot_n_jump_matrix
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:490
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
LocalIntegrators::Divergence::u_times_n_residual
void u_times_n_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< double > &data, double factor=1.)
Definition: divergence.h:353
LocalIntegrators::Divergence::cell_residual
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, const double factor=1.)
Definition: divergence.h:101
FEValuesBase::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_values.h:2108
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1607
LocalIntegrators::Divergence::grad_div_matrix
void grad_div_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: divergence.h:449
FEValuesBase::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
Tensor< 1, dim >
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
FEValuesBase::shape_grad_component
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:98
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
LocalIntegrators::GradDiv::cell_residual
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, const double factor=1.)
Definition: grad_div.h:94
FEValuesBase::JxW
double JxW(const unsigned int quadrature_point) const
LocalIntegrators
Library of integrals over cells and faces.
Definition: advection.h:34
dof_info.h
FEValuesBase
Definition: fe.h:36
exceptions.h
mapping.h
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
FEValuesBase::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
grad_div.h
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
quadrature.h
LocalIntegrators::Divergence::gradient_residual
void gradient_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< Tensor< 1, dim >> &input, const double factor=1.)
Definition: divergence.h:214
FEValuesBase::n_quadrature_points
const unsigned int n_quadrature_points
Definition: fe_values.h:2101
config.h
FullMatrix< double >
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Vector< number >
full_matrix.h
Generated by   doxygen 1.8.17