Reference documentation for deal.II version 9.3.0
divergence.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2020 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
24
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27
29
31
33
35
36 namespace LocalIntegrators
37 {
44  namespace Divergence
45  {
52  template <int dim>
53  void
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
91  template <int dim, typename number>
92  void
94  const FEValuesBase<dim> & fetest,
95  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
96  const double factor = 1.)
97  {
98  AssertDimension(fetest.get_fe().n_components(), 1);
100  const unsigned int t_dofs = fetest.dofs_per_cell;
101  Assert(result.size() == t_dofs,
102  ExcDimensionMismatch(result.size(), t_dofs));
103
104  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
105  {
106  const double dx = factor * fetest.JxW(k);
107
108  for (unsigned int i = 0; i < t_dofs; ++i)
109  for (unsigned int d = 0; d < dim; ++d)
110  result(i) += dx * input[d][k][d] * fetest.shape_value(i, k);
111  }
112  }
113
114
124  template <int dim, typename number>
125  void
127  const FEValuesBase<dim> & fetest,
128  const ArrayView<const std::vector<double>> &input,
129  const double factor = 1.)
130  {
131  AssertDimension(fetest.get_fe().n_components(), 1);
133  const unsigned int t_dofs = fetest.dofs_per_cell;
134  Assert(result.size() == t_dofs,
135  ExcDimensionMismatch(result.size(), t_dofs));
136
137  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
138  {
139  const double dx = factor * fetest.JxW(k);
140
141  for (unsigned int i = 0; i < t_dofs; ++i)
142  for (unsigned int d = 0; d < dim; ++d)
143  result(i) -= dx * input[d][k] * fetest.shape_grad(i, k)[d];
144  }
145  }
146
147
155  template <int dim>
156  void
158  const FEValuesBase<dim> &fe,
159  const FEValuesBase<dim> &fetest,
160  double factor = 1.)
161  {
162  const unsigned int t_dofs = fetest.dofs_per_cell;
163  const unsigned int n_dofs = fe.dofs_per_cell;
164
165  AssertDimension(fetest.get_fe().n_components(), dim);
166  AssertDimension(fe.get_fe().n_components(), 1);
167  AssertDimension(M.m(), t_dofs);
168  AssertDimension(M.n(), n_dofs);
169
170  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
171  {
172  const double dx = fe.JxW(k) * factor;
173  for (unsigned int d = 0; d < dim; ++d)
174  for (unsigned int i = 0; i < t_dofs; ++i)
175  {
176  const double vv = fetest.shape_value_component(i, k, d);
177  for (unsigned int j = 0; j < n_dofs; ++j)
178  {
179  const Tensor<1, dim> &Du = fe.shape_grad(j, k);
180  M(i, j) += dx * vv * Du[d];
181  }
182  }
183  }
184  }
185
195  template <int dim, typename number>
196  void
198  const FEValuesBase<dim> & fetest,
199  const std::vector<Tensor<1, dim>> &input,
200  const double factor = 1.)
201  {
202  AssertDimension(fetest.get_fe().n_components(), dim);
204  const unsigned int t_dofs = fetest.dofs_per_cell;
205  Assert(result.size() == t_dofs,
206  ExcDimensionMismatch(result.size(), t_dofs));
207
208  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
209  {
210  const double dx = factor * fetest.JxW(k);
211
212  for (unsigned int i = 0; i < t_dofs; ++i)
213  for (unsigned int d = 0; d < dim; ++d)
214  result(i) +=
215  dx * input[k][d] * fetest.shape_value_component(i, k, d);
216  }
217  }
218
228  template <int dim, typename number>
229  void
231  const FEValuesBase<dim> & fetest,
232  const std::vector<double> &input,
233  const double factor = 1.)
234  {
235  AssertDimension(fetest.get_fe().n_components(), dim);
237  const unsigned int t_dofs = fetest.dofs_per_cell;
238  Assert(result.size() == t_dofs,
239  ExcDimensionMismatch(result.size(), t_dofs));
240
241  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
242  {
243  const double dx = factor * fetest.JxW(k);
244
245  for (unsigned int i = 0; i < t_dofs; ++i)
246  for (unsigned int d = 0; d < dim; ++d)
247  result(i) -=
248  dx * input[k] * fetest.shape_grad_component(i, k, d)[d];
249  }
250  }
251
257  template <int dim>
258  void
260  const FEValuesBase<dim> &fe,
261  const FEValuesBase<dim> &fetest,
262  double factor = 1.)
263  {
264  const unsigned int n_dofs = fe.dofs_per_cell;
265  const unsigned int t_dofs = fetest.dofs_per_cell;
266
267  AssertDimension(fe.get_fe().n_components(), dim);
268  AssertDimension(fetest.get_fe().n_components(), 1);
269  AssertDimension(M.m(), t_dofs);
270  AssertDimension(M.n(), n_dofs);
271
272  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
273  {
274  const Tensor<1, dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
275  for (unsigned int i = 0; i < t_dofs; ++i)
276  for (unsigned int j = 0; j < n_dofs; ++j)
277  for (unsigned int d = 0; d < dim; ++d)
278  M(i, j) += ndx[d] * fe.shape_value_component(j, k, d) *
279  fetest.shape_value(i, k);
280  }
281  }
282
290  template <int dim, typename number>
291  void
293  const FEValuesBase<dim> & fe,
294  const FEValuesBase<dim> & fetest,
295  const ArrayView<const std::vector<double>> &data,
296  double factor = 1.)
297  {
298  const unsigned int t_dofs = fetest.dofs_per_cell;
299
300  AssertDimension(fe.get_fe().n_components(), dim);
301  AssertDimension(fetest.get_fe().n_components(), 1);
302  AssertDimension(result.size(), t_dofs);
304
305  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
306  {
307  const Tensor<1, dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
308
309  for (unsigned int i = 0; i < t_dofs; ++i)
310  for (unsigned int d = 0; d < dim; ++d)
311  result(i) += ndx[d] * fetest.shape_value(i, k) * data[d][k];
312  }
313  }
314
322  template <int dim, typename number>
323  void
325  const FEValuesBase<dim> & fetest,
326  const std::vector<double> &data,
327  double factor = 1.)
328  {
329  const unsigned int t_dofs = fetest.dofs_per_cell;
330
331  AssertDimension(fetest.get_fe().n_components(), dim);
332  AssertDimension(result.size(), t_dofs);
334
335  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
336  {
337  const Tensor<1, dim> ndx =
338  factor * fetest.normal_vector(k) * fetest.JxW(k);
339
340  for (unsigned int i = 0; i < t_dofs; ++i)
341  for (unsigned int d = 0; d < dim; ++d)
342  result(i) +=
343  ndx[d] * fetest.shape_value_component(i, k, d) * data[k];
344  }
345  }
346
356  template <int dim>
357  void
359  FullMatrix<double> & M12,
360  FullMatrix<double> & M21,
361  FullMatrix<double> & M22,
362  const FEValuesBase<dim> &fe1,
363  const FEValuesBase<dim> &fe2,
364  const FEValuesBase<dim> &fetest1,
365  const FEValuesBase<dim> &fetest2,
366  double factor = 1.)
367  {
368  const unsigned int n_dofs = fe1.dofs_per_cell;
369  const unsigned int t_dofs = fetest1.dofs_per_cell;
370
371  AssertDimension(fe1.get_fe().n_components(), dim);
372  AssertDimension(fe2.get_fe().n_components(), dim);
373  AssertDimension(fetest1.get_fe().n_components(), 1);
374  AssertDimension(fetest2.get_fe().n_components(), 1);
375  AssertDimension(M11.m(), t_dofs);
376  AssertDimension(M11.n(), n_dofs);
377  AssertDimension(M12.m(), t_dofs);
378  AssertDimension(M12.n(), n_dofs);
379  AssertDimension(M21.m(), t_dofs);
380  AssertDimension(M21.n(), n_dofs);
381  AssertDimension(M22.m(), t_dofs);
382  AssertDimension(M22.n(), n_dofs);
383
384  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
385  {
386  const double dx = factor * fe1.JxW(k);
387  for (unsigned int i = 0; i < t_dofs; ++i)
388  for (unsigned int j = 0; j < n_dofs; ++j)
389  for (unsigned int d = 0; d < dim; ++d)
390  {
391  const double un1 = fe1.shape_value_component(j, k, d) *
392  fe1.normal_vector(k)[d];
393  const double un2 = -fe2.shape_value_component(j, k, d) *
394  fe1.normal_vector(k)[d];
395  const double v1 = fetest1.shape_value(i, k);
396  const double v2 = fetest2.shape_value(i, k);
397
398  M11(i, j) += .5 * dx * un1 * v1;
399  M12(i, j) += .5 * dx * un2 * v1;
400  M21(i, j) += .5 * dx * un1 * v2;
401  M22(i, j) += .5 * dx * un2 * v2;
402  }
403  }
404  }
405
415  template <int dim>
416  void
418  FullMatrix<double> & M12,
419  FullMatrix<double> & M21,
420  FullMatrix<double> & M22,
421  const FEValuesBase<dim> &fe1,
422  const FEValuesBase<dim> &fe2,
423  double factor = 1.)
424  {
425  const unsigned int n_dofs = fe1.dofs_per_cell;
426
427  AssertDimension(fe1.get_fe().n_components(), dim);
428  AssertDimension(fe2.get_fe().n_components(), dim);
429  AssertDimension(M11.m(), n_dofs);
430  AssertDimension(M11.n(), n_dofs);
431  AssertDimension(M12.m(), n_dofs);
432  AssertDimension(M12.n(), n_dofs);
433  AssertDimension(M21.m(), n_dofs);
434  AssertDimension(M21.n(), n_dofs);
435  AssertDimension(M22.m(), n_dofs);
436  AssertDimension(M22.n(), n_dofs);
437
438  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
439  {
440  const double dx = factor * fe1.JxW(k);
441  for (unsigned int i = 0; i < n_dofs; ++i)
442  for (unsigned int j = 0; j < n_dofs; ++j)
443  for (unsigned int d = 0; d < dim; ++d)
444  {
445  const double un1 = fe1.shape_value_component(j, k, d) *
446  fe1.normal_vector(k)[d];
447  const double un2 = -fe2.shape_value_component(j, k, d) *
448  fe1.normal_vector(k)[d];
449  const double vn1 = fe1.shape_value_component(i, k, d) *
450  fe1.normal_vector(k)[d];
451  const double vn2 = -fe2.shape_value_component(i, k, d) *
452  fe1.normal_vector(k)[d];
453
454  M11(i, j) += dx * un1 * vn1;
455  M12(i, j) += dx * un2 * vn1;
456  M21(i, j) += dx * un1 * vn2;
457  M22(i, j) += dx * un2 * vn2;
458  }
459  }
460  }
461
470  template <int dim>
471  double
473  const ArrayView<const std::vector<Tensor<1, dim>>> &Du)
474  {
475  AssertDimension(fe.get_fe().n_components(), dim);
477
478  double result = 0;
479  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
480  {
481  double div = Du[0][k][0];
482  for (unsigned int d = 1; d < dim; ++d)
483  div += Du[d][k][d];
484  result += div * div * fe.JxW(k);
485  }
486  return result;
487  }
488
489  } // namespace Divergence
490 } // namespace LocalIntegrators
491
492
494
495 #endif
size_type m() const
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:93
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
const unsigned int dofs_per_cell
Definition: fe_values.h:2450
const FiniteElement< dim, spacedim > & get_fe() const
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1649
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
const unsigned int v1
Definition: grid_tools.cc:963
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.
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:292
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:54
size_type n() const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) 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:197
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:157
void u_dot_n_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:259
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:417
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
void u_times_n_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< double > &data, double factor=1.)
Definition: divergence.h:324
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
double JxW(const unsigned int quadrature_point) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
size_type size() const
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