Reference documentation for deal.II version 9.4.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 - 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
26#include <deal.II/fe/mapping.h>
27
29
31
33
35
36namespace 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;
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);
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);
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
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
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
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
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);
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
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 {
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
const unsigned int dofs_per_cell
Definition: fe_values.h:2450
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
double JxW(const unsigned int quadrature_point) const
const FiniteElement< dim, spacedim > & get_fe() const
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
unsigned int n_components() const
size_type n() const
size_type m() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
const unsigned int v1
Definition: grid_tools.cc:1000
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1690
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type size() const
void u_times_n_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const std::vector< double > &data, double factor=1.)
Definition: divergence.h:324
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
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
void u_dot_n_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:259
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
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 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
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:157
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: divergence.h:54
Library of integrals over cells and faces.
Definition: advection.h:35
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)