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