Reference documentation for deal.II version 9.5.1
\(\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\}}\)
Loading...
Searching...
No Matches
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
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
34DEAL_II_NAMESPACE_OPEN
35
36namespace LocalIntegrators
37{
44 namespace Divergence
45 {
52 template <int dim>
53 void
54 cell_matrix(FullMatrix<double> & M,
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
93 cell_residual(Vector<number> & result,
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);
99 AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
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
126 cell_residual(Vector<number> & result,
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);
132 AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
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
157 gradient_matrix(FullMatrix<double> & M,
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
197 gradient_residual(Vector<number> & result,
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);
203 AssertDimension(input.size(), fetest.n_quadrature_points);
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
230 gradient_residual(Vector<number> & result,
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);
236 AssertDimension(input.size(), fetest.n_quadrature_points);
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
259 u_dot_n_matrix(FullMatrix<double> & M,
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
292 u_dot_n_residual(Vector<number> & result,
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);
303 AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
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
324 u_times_n_residual(Vector<number> & result,
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);
333 AssertDimension(data.size(), fetest.n_quadrature_points);
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
358 u_dot_n_matrix(FullMatrix<double> & M11,
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
417 u_dot_n_jump_matrix(FullMatrix<double> & M11,
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
472 norm(const FEValuesBase<dim> & fe,
473 const ArrayView<const std::vector<Tensor<1, dim>>> &Du)
474 {
475 AssertDimension(fe.get_fe().n_components(), dim);
476 AssertVectorVectorDimension(Du, dim, fe.n_quadrature_points);
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
493DEAL_II_NAMESPACE_CLOSE
494
495#endif
ArrayView::size
std::size_t size() const
Definition array_view.h:576
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
fe_values.h
full_matrix.h
grad_div.h
v1
const unsigned int v1
Definition grid_tools.cc:1067
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition exceptions.h:1810
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
dof_info.h
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:324
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
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:93
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:259
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:197
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:292
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:417
LocalIntegrators::Divergence::gradient_matrix
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition divergence.h:157
LocalIntegrators::Divergence::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition divergence.h:54
LocalIntegrators
Library of integrals over cells and faces.
Definition advection.h:35
Generated by doxygen 1.9.8