deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\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// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2010 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_integrators_divergence_h
14#define dealii_integrators_divergence_h
15
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/exceptions.h>
20#include <deal.II/base/quadrature.h>
21
22#include <deal.II/fe/fe_values.h>
23#include <deal.II/fe/mapping.h>
24
25#include <deal.II/integrators/grad_div.h>
26
27#include <deal.II/lac/full_matrix.h>
28
29#include <deal.II/meshworker/dof_info.h>
30
31DEAL_II_NAMESPACE_OPEN
32
33namespace LocalIntegrators
34{
41 namespace Divergence
42 {
49 template <int dim>
50 DEAL_II_DEPRECATED void
51 cell_matrix(FullMatrix<double> &M,
52 const FEValuesBase<dim> &fe,
53 const FEValuesBase<dim> &fetest,
54 double factor = 1.)
55 {
56 const unsigned int n_dofs = fe.dofs_per_cell;
57 const unsigned int t_dofs = fetest.dofs_per_cell;
58 AssertDimension(fe.get_fe().n_components(), dim);
59 AssertDimension(fetest.get_fe().n_components(), 1);
60 AssertDimension(M.m(), t_dofs);
61 AssertDimension(M.n(), n_dofs);
62
63 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
64 {
65 const double dx = fe.JxW(k) * factor;
66 for (unsigned int i = 0; i < t_dofs; ++i)
67 {
68 const double vv = fetest.shape_value(i, k);
69 for (unsigned int d = 0; d < dim; ++d)
70 for (unsigned int j = 0; j < n_dofs; ++j)
71 {
72 const double du = fe.shape_grad_component(j, k, d)[d];
73 M(i, j) += dx * du * vv;
74 }
75 }
76 }
77 }
78
88 template <int dim, typename number>
89 DEAL_II_DEPRECATED void
90 cell_residual(Vector<number> &result,
91 const FEValuesBase<dim> &fetest,
92 const ArrayView<const std::vector<Tensor<1, dim>>> &input,
93 const double factor = 1.)
94 {
95 AssertDimension(fetest.get_fe().n_components(), 1);
96 AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
97 const unsigned int t_dofs = fetest.dofs_per_cell;
98 Assert(result.size() == t_dofs,
99 ExcDimensionMismatch(result.size(), t_dofs));
100
101 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
102 {
103 const double dx = factor * fetest.JxW(k);
104
105 for (unsigned int i = 0; i < t_dofs; ++i)
106 for (unsigned int d = 0; d < dim; ++d)
107 result(i) += dx * input[d][k][d] * fetest.shape_value(i, k);
108 }
109 }
110
111
121 template <int dim, typename number>
122 DEAL_II_DEPRECATED void
123 cell_residual(Vector<number> &result,
124 const FEValuesBase<dim> &fetest,
125 const ArrayView<const std::vector<double>> &input,
126 const double factor = 1.)
127 {
128 AssertDimension(fetest.get_fe().n_components(), 1);
129 AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
130 const unsigned int t_dofs = fetest.dofs_per_cell;
131 Assert(result.size() == t_dofs,
132 ExcDimensionMismatch(result.size(), t_dofs));
133
134 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
135 {
136 const double dx = factor * fetest.JxW(k);
137
138 for (unsigned int i = 0; i < t_dofs; ++i)
139 for (unsigned int d = 0; d < dim; ++d)
140 result(i) -= dx * input[d][k] * fetest.shape_grad(i, k)[d];
141 }
142 }
143
144
152 template <int dim>
153 DEAL_II_DEPRECATED void
154 gradient_matrix(FullMatrix<double> &M,
155 const FEValuesBase<dim> &fe,
156 const FEValuesBase<dim> &fetest,
157 double factor = 1.)
158 {
159 const unsigned int t_dofs = fetest.dofs_per_cell;
160 const unsigned int n_dofs = fe.dofs_per_cell;
161
162 AssertDimension(fetest.get_fe().n_components(), dim);
163 AssertDimension(fe.get_fe().n_components(), 1);
164 AssertDimension(M.m(), t_dofs);
165 AssertDimension(M.n(), n_dofs);
166
167 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
168 {
169 const double dx = fe.JxW(k) * factor;
170 for (unsigned int d = 0; d < dim; ++d)
171 for (unsigned int i = 0; i < t_dofs; ++i)
172 {
173 const double vv = fetest.shape_value_component(i, k, d);
174 for (unsigned int j = 0; j < n_dofs; ++j)
175 {
176 const Tensor<1, dim> &Du = fe.shape_grad(j, k);
177 M(i, j) += dx * vv * Du[d];
178 }
179 }
180 }
181 }
182
192 template <int dim, typename number>
193 DEAL_II_DEPRECATED void
194 gradient_residual(Vector<number> &result,
195 const FEValuesBase<dim> &fetest,
196 const std::vector<Tensor<1, dim>> &input,
197 const double factor = 1.)
198 {
199 AssertDimension(fetest.get_fe().n_components(), dim);
200 AssertDimension(input.size(), fetest.n_quadrature_points);
201 const unsigned int t_dofs = fetest.dofs_per_cell;
202 Assert(result.size() == t_dofs,
203 ExcDimensionMismatch(result.size(), t_dofs));
204
205 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
206 {
207 const double dx = factor * fetest.JxW(k);
208
209 for (unsigned int i = 0; i < t_dofs; ++i)
210 for (unsigned int d = 0; d < dim; ++d)
211 result(i) +=
212 dx * input[k][d] * fetest.shape_value_component(i, k, d);
213 }
214 }
215
225 template <int dim, typename number>
226 DEAL_II_DEPRECATED void
227 gradient_residual(Vector<number> &result,
228 const FEValuesBase<dim> &fetest,
229 const std::vector<double> &input,
230 const double factor = 1.)
231 {
232 AssertDimension(fetest.get_fe().n_components(), dim);
233 AssertDimension(input.size(), fetest.n_quadrature_points);
234 const unsigned int t_dofs = fetest.dofs_per_cell;
235 Assert(result.size() == t_dofs,
236 ExcDimensionMismatch(result.size(), t_dofs));
237
238 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
239 {
240 const double dx = factor * fetest.JxW(k);
241
242 for (unsigned int i = 0; i < t_dofs; ++i)
243 for (unsigned int d = 0; d < dim; ++d)
244 result(i) -=
245 dx * input[k] * fetest.shape_grad_component(i, k, d)[d];
246 }
247 }
248
254 template <int dim>
255 DEAL_II_DEPRECATED void
256 u_dot_n_matrix(FullMatrix<double> &M,
257 const FEValuesBase<dim> &fe,
258 const FEValuesBase<dim> &fetest,
259 double factor = 1.)
260 {
261 const unsigned int n_dofs = fe.dofs_per_cell;
262 const unsigned int t_dofs = fetest.dofs_per_cell;
263
264 AssertDimension(fe.get_fe().n_components(), dim);
265 AssertDimension(fetest.get_fe().n_components(), 1);
266 AssertDimension(M.m(), t_dofs);
267 AssertDimension(M.n(), n_dofs);
268
269 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
270 {
271 const Tensor<1, dim> ndx = factor * fe.JxW(k) * fe.normal_vector(k);
272 for (unsigned int i = 0; i < t_dofs; ++i)
273 for (unsigned int j = 0; j < n_dofs; ++j)
274 for (unsigned int d = 0; d < dim; ++d)
275 M(i, j) += ndx[d] * fe.shape_value_component(j, k, d) *
276 fetest.shape_value(i, k);
277 }
278 }
279
287 template <int dim, typename number>
288 DEAL_II_DEPRECATED void
289 u_dot_n_residual(Vector<number> &result,
290 const FEValuesBase<dim> &fe,
291 const FEValuesBase<dim> &fetest,
292 const ArrayView<const std::vector<double>> &data,
293 double factor = 1.)
294 {
295 const unsigned int t_dofs = fetest.dofs_per_cell;
296
297 AssertDimension(fe.get_fe().n_components(), dim);
298 AssertDimension(fetest.get_fe().n_components(), 1);
299 AssertDimension(result.size(), t_dofs);
300 AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
301
302 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
303 {
304 const Tensor<1, dim> ndx = factor * fe.normal_vector(k) * fe.JxW(k);
305
306 for (unsigned int i = 0; i < t_dofs; ++i)
307 for (unsigned int d = 0; d < dim; ++d)
308 result(i) += ndx[d] * fetest.shape_value(i, k) * data[d][k];
309 }
310 }
311
319 template <int dim, typename number>
320 DEAL_II_DEPRECATED void
321 u_times_n_residual(Vector<number> &result,
322 const FEValuesBase<dim> &fetest,
323 const std::vector<double> &data,
324 double factor = 1.)
325 {
326 const unsigned int t_dofs = fetest.dofs_per_cell;
327
328 AssertDimension(fetest.get_fe().n_components(), dim);
329 AssertDimension(result.size(), t_dofs);
330 AssertDimension(data.size(), fetest.n_quadrature_points);
331
332 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
333 {
334 const Tensor<1, dim> ndx =
335 factor * fetest.normal_vector(k) * fetest.JxW(k);
336
337 for (unsigned int i = 0; i < t_dofs; ++i)
338 for (unsigned int d = 0; d < dim; ++d)
339 result(i) +=
340 ndx[d] * fetest.shape_value_component(i, k, d) * data[k];
341 }
342 }
343
353 template <int dim>
354 DEAL_II_DEPRECATED void
355 u_dot_n_matrix(FullMatrix<double> &M11,
356 FullMatrix<double> &M12,
357 FullMatrix<double> &M21,
358 FullMatrix<double> &M22,
359 const FEValuesBase<dim> &fe1,
360 const FEValuesBase<dim> &fe2,
361 const FEValuesBase<dim> &fetest1,
362 const FEValuesBase<dim> &fetest2,
363 double factor = 1.)
364 {
365 const unsigned int n_dofs = fe1.dofs_per_cell;
366 const unsigned int t_dofs = fetest1.dofs_per_cell;
367
368 AssertDimension(fe1.get_fe().n_components(), dim);
369 AssertDimension(fe2.get_fe().n_components(), dim);
370 AssertDimension(fetest1.get_fe().n_components(), 1);
371 AssertDimension(fetest2.get_fe().n_components(), 1);
372 AssertDimension(M11.m(), t_dofs);
373 AssertDimension(M11.n(), n_dofs);
374 AssertDimension(M12.m(), t_dofs);
375 AssertDimension(M12.n(), n_dofs);
376 AssertDimension(M21.m(), t_dofs);
377 AssertDimension(M21.n(), n_dofs);
378 AssertDimension(M22.m(), t_dofs);
379 AssertDimension(M22.n(), n_dofs);
380
381 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
382 {
383 const double dx = factor * fe1.JxW(k);
384 for (unsigned int i = 0; i < t_dofs; ++i)
385 for (unsigned int j = 0; j < n_dofs; ++j)
386 for (unsigned int d = 0; d < dim; ++d)
387 {
388 const double un1 = fe1.shape_value_component(j, k, d) *
389 fe1.normal_vector(k)[d];
390 const double un2 = -fe2.shape_value_component(j, k, d) *
391 fe1.normal_vector(k)[d];
392 const double v1 = fetest1.shape_value(i, k);
393 const double v2 = fetest2.shape_value(i, k);
394
395 M11(i, j) += .5 * dx * un1 * v1;
396 M12(i, j) += .5 * dx * un2 * v1;
397 M21(i, j) += .5 * dx * un1 * v2;
398 M22(i, j) += .5 * dx * un2 * v2;
399 }
400 }
401 }
402
412 template <int dim>
413 DEAL_II_DEPRECATED void
414 u_dot_n_jump_matrix(FullMatrix<double> &M11,
415 FullMatrix<double> &M12,
416 FullMatrix<double> &M21,
417 FullMatrix<double> &M22,
418 const FEValuesBase<dim> &fe1,
419 const FEValuesBase<dim> &fe2,
420 double factor = 1.)
421 {
422 const unsigned int n_dofs = fe1.dofs_per_cell;
423
424 AssertDimension(fe1.get_fe().n_components(), dim);
425 AssertDimension(fe2.get_fe().n_components(), dim);
426 AssertDimension(M11.m(), n_dofs);
427 AssertDimension(M11.n(), n_dofs);
428 AssertDimension(M12.m(), n_dofs);
429 AssertDimension(M12.n(), n_dofs);
430 AssertDimension(M21.m(), n_dofs);
431 AssertDimension(M21.n(), n_dofs);
432 AssertDimension(M22.m(), n_dofs);
433 AssertDimension(M22.n(), n_dofs);
434
435 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
436 {
437 const double dx = factor * fe1.JxW(k);
438 for (unsigned int i = 0; i < n_dofs; ++i)
439 for (unsigned int j = 0; j < n_dofs; ++j)
440 for (unsigned int d = 0; d < dim; ++d)
441 {
442 const double un1 = fe1.shape_value_component(j, k, d) *
443 fe1.normal_vector(k)[d];
444 const double un2 = -fe2.shape_value_component(j, k, d) *
445 fe1.normal_vector(k)[d];
446 const double vn1 = fe1.shape_value_component(i, k, d) *
447 fe1.normal_vector(k)[d];
448 const double vn2 = -fe2.shape_value_component(i, k, d) *
449 fe1.normal_vector(k)[d];
450
451 M11(i, j) += dx * un1 * vn1;
452 M12(i, j) += dx * un2 * vn1;
453 M21(i, j) += dx * un1 * vn2;
454 M22(i, j) += dx * un2 * vn2;
455 }
456 }
457 }
458
467 template <int dim>
468 DEAL_II_DEPRECATED double
469 norm(const FEValuesBase<dim> &fe,
470 const ArrayView<const std::vector<Tensor<1, dim>>> &Du)
471 {
472 AssertDimension(fe.get_fe().n_components(), dim);
473 AssertVectorVectorDimension(Du, dim, fe.n_quadrature_points);
474
475 double result = 0;
476 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
477 {
478 double div = Du[0][k][0];
479 for (unsigned int d = 1; d < dim; ++d)
480 div += Du[d][k][d];
481 result += div * div * fe.JxW(k);
482 }
483 return result;
484 }
485
486 } // namespace Divergence
487} // namespace LocalIntegrators
488
489
490DEAL_II_NAMESPACE_CLOSE
491
492#endif
FEValuesBase::dofs_per_cell
const unsigned int dofs_per_cell
Definition fe_values_base.h:192
FEValuesBase::normal_vector
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
FEValuesBase::shape_value_component
double shape_value_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
FEValuesBase::n_quadrature_points
const unsigned int n_quadrature_points
Definition fe_values_base.h:174
FEValuesBase::shape_grad_component
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
FEValuesBase::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
FEValuesBase::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
FEValuesBase::JxW
double JxW(const unsigned int q_point) const
FEValuesBase::shape_value
const double & shape_value(const unsigned int i, const unsigned int q_point) const
FiniteElementData::n_components
unsigned int n_components() const
FullMatrix::n
size_type n() const
FullMatrix::m
size_type m() const
Vector::size
virtual size_type size() const override
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition config.h:294
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
fe_values.h
full_matrix.h
grad_div.h
v1
const unsigned int v1
Definition grid_tools_topology.cc:648
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition exception_macros.h:873
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
dof_info.h
data
std::vector< index_type > data
Definition mpi.cc:734
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:321
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:469
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:90
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:256
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:194
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:289
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:414
LocalIntegrators::Divergence::gradient_matrix
void gradient_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition divergence.h:154
LocalIntegrators::Divergence::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition divergence.h:51
LocalIntegrators
Library of integrals over cells and faces.
Definition advection.h:32
Generated by doxygen 1.9.8