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
advection.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) 2012 - 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_advection_h
14#define dealii_integrators_advection_h
15
16
17#include <deal.II/base/config.h>
18
21
23#include <deal.II/fe/mapping.h>
24
26
28
30
32{
47 namespace Advection
48 {
70 template <int dim>
73 const FEValuesBase<dim> &fe,
74 const FEValuesBase<dim> &fetest,
75 const ArrayView<const std::vector<double>> &velocity,
76 const double factor = 1.)
77 {
78 const unsigned int n_dofs = fe.dofs_per_cell;
79 const unsigned int t_dofs = fetest.dofs_per_cell;
80 const unsigned int n_components = fe.get_fe().n_components();
81
82 AssertDimension(velocity.size(), dim);
83 // If the size of the
84 // velocity vectors is one,
85 // then do not increment
86 // between quadrature points.
87 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
88
89 if (v_increment == 1)
90 {
92 }
93
94 AssertDimension(M.n(), n_dofs);
95 AssertDimension(M.m(), t_dofs);
96
97 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
98 {
99 const double dx = factor * fe.JxW(k);
100 const unsigned int vindex = k * v_increment;
101
102 for (unsigned j = 0; j < n_dofs; ++j)
103 for (unsigned i = 0; i < t_dofs; ++i)
104 for (unsigned int c = 0; c < n_components; ++c)
105 {
106 double wgradv =
107 velocity[0][vindex] * fe.shape_grad_component(i, k, c)[0];
108 for (unsigned int d = 1; d < dim; ++d)
109 wgradv +=
110 velocity[d][vindex] * fe.shape_grad_component(i, k, c)[d];
111 M(i, j) -= dx * wgradv * fe.shape_value_component(j, k, c);
112 }
113 }
114 }
115
116
117
126 template <int dim>
127 DEAL_II_DEPRECATED inline void
129 const FEValuesBase<dim> &fe,
130 const std::vector<Tensor<1, dim>> &input,
131 const ArrayView<const std::vector<double>> &velocity,
132 double factor = 1.)
133 {
134 const unsigned int nq = fe.n_quadrature_points;
135 const unsigned int n_dofs = fe.dofs_per_cell;
136 Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
137 Assert(result.size() == n_dofs,
138 ExcDimensionMismatch(result.size(), n_dofs));
139
140 AssertDimension(velocity.size(), dim);
141 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
142 if (v_increment == 1)
143 {
145 }
146
147 for (unsigned k = 0; k < nq; ++k)
148 {
149 const double dx = factor * fe.JxW(k);
150 for (unsigned i = 0; i < n_dofs; ++i)
151 for (unsigned int d = 0; d < dim; ++d)
152 result(i) += dx * input[k][d] * fe.shape_value(i, k) *
153 velocity[d][k * v_increment];
154 }
155 }
156
157
158
169 template <int dim>
170 DEAL_II_DEPRECATED inline void
172 const FEValuesBase<dim> &fe,
173 const ArrayView<const std::vector<Tensor<1, dim>>> &input,
174 const ArrayView<const std::vector<double>> &velocity,
175 double factor = 1.)
176 {
177 const unsigned int nq = fe.n_quadrature_points;
178 const unsigned int n_dofs = fe.dofs_per_cell;
179 const unsigned int n_comp = fe.get_fe().n_components();
180
182 Assert(result.size() == n_dofs,
183 ExcDimensionMismatch(result.size(), n_dofs));
184
185 AssertDimension(velocity.size(), dim);
186 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
187 if (v_increment == 1)
188 {
190 }
191
192 for (unsigned k = 0; k < nq; ++k)
193 {
194 const double dx = factor * fe.JxW(k);
195 for (unsigned i = 0; i < n_dofs; ++i)
196 for (unsigned int c = 0; c < n_comp; ++c)
197 for (unsigned int d = 0; d < dim; ++d)
198 result(i) += dx * input[c][k][d] *
199 fe.shape_value_component(i, k, c) *
200 velocity[d][k * v_increment];
201 }
202 }
203
204
205
211 template <int dim>
212 DEAL_II_DEPRECATED inline void
214 const FEValuesBase<dim> &fe,
215 const std::vector<double> &input,
216 const ArrayView<const std::vector<double>> &velocity,
217 double factor = 1.)
218 {
219 const unsigned int nq = fe.n_quadrature_points;
220 const unsigned int n_dofs = fe.dofs_per_cell;
221 Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
222 Assert(result.size() == n_dofs,
223 ExcDimensionMismatch(result.size(), n_dofs));
224
225 AssertDimension(velocity.size(), dim);
226 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
227 if (v_increment == 1)
228 {
230 }
231
232 for (unsigned k = 0; k < nq; ++k)
233 {
234 const double dx = factor * fe.JxW(k);
235 for (unsigned i = 0; i < n_dofs; ++i)
236 for (unsigned int d = 0; d < dim; ++d)
237 result(i) -= dx * input[k] * fe.shape_grad(i, k)[d] *
238 velocity[d][k * v_increment];
239 }
240 }
241
242
243
251 template <int dim>
252 DEAL_II_DEPRECATED inline void
254 const FEValuesBase<dim> &fe,
255 const ArrayView<const std::vector<double>> &input,
256 const ArrayView<const std::vector<double>> &velocity,
257 double factor = 1.)
258 {
259 const unsigned int nq = fe.n_quadrature_points;
260 const unsigned int n_dofs = fe.dofs_per_cell;
261 const unsigned int n_comp = fe.get_fe().n_components();
262
264 Assert(result.size() == n_dofs,
265 ExcDimensionMismatch(result.size(), n_dofs));
266
267 AssertDimension(velocity.size(), dim);
268 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
269 if (v_increment == 1)
270 {
272 }
273
274 for (unsigned k = 0; k < nq; ++k)
275 {
276 const double dx = factor * fe.JxW(k);
277 for (unsigned i = 0; i < n_dofs; ++i)
278 for (unsigned int c = 0; c < n_comp; ++c)
279 for (unsigned int d = 0; d < dim; ++d)
280 result(i) -= dx * input[c][k] *
281 fe.shape_grad_component(i, k, c)[d] *
282 velocity[d][k * v_increment];
283 }
284 }
285
286
287
305 template <int dim>
308 const FEValuesBase<dim> &fe,
309 const FEValuesBase<dim> &fetest,
310 const ArrayView<const std::vector<double>> &velocity,
311 double factor = 1.)
312 {
313 const unsigned int n_dofs = fe.dofs_per_cell;
314 const unsigned int t_dofs = fetest.dofs_per_cell;
315 unsigned int n_components = fe.get_fe().n_components();
316 AssertDimension(M.m(), n_dofs);
317 AssertDimension(M.n(), n_dofs);
318
319 AssertDimension(velocity.size(), dim);
320 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
321 if (v_increment == 1)
322 {
324 }
325
326 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
327 {
328 const double dx = factor * fe.JxW(k);
329
330 double nv = 0.;
331 for (unsigned int d = 0; d < dim; ++d)
332 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
333
334 if (nv > 0)
335 {
336 for (unsigned i = 0; i < t_dofs; ++i)
337 for (unsigned j = 0; j < n_dofs; ++j)
338 {
339 if (fe.get_fe().is_primitive())
340 M(i, j) +=
341 dx * nv * fe.shape_value(i, k) * fe.shape_value(j, k);
342 else
343 for (unsigned int c = 0; c < n_components; ++c)
344 M(i, j) += dx * nv *
345 fetest.shape_value_component(i, k, c) *
346 fe.shape_value_component(j, k, c);
347 }
348 }
349 }
350 }
351
352
353
378 template <int dim>
379 DEAL_II_DEPRECATED inline void
381 const FEValuesBase<dim> &fe,
382 const std::vector<double> &input,
383 const std::vector<double> &data,
384 const ArrayView<const std::vector<double>> &velocity,
385 double factor = 1.)
386 {
387 const unsigned int n_dofs = fe.dofs_per_cell;
388
389 AssertDimension(input.size(), fe.n_quadrature_points);
391
392 AssertDimension(velocity.size(), dim);
393 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
394 if (v_increment == 1)
395 {
397 }
398
399
400 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
401 {
402 const double dx = factor * fe.JxW(k);
403
404 double nv = 0.;
405 for (unsigned int d = 0; d < dim; ++d)
406 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
407
408 // Always use the upwind value
409 const double val = (nv > 0.) ? input[k] : -data[k];
410
411 for (unsigned i = 0; i < n_dofs; ++i)
412 {
413 const double v = fe.shape_value(i, k);
414 result(i) += dx * nv * val * v;
415 }
416 }
417 }
418
419
420
445 template <int dim>
446 DEAL_II_DEPRECATED inline void
448 const FEValuesBase<dim> &fe,
449 const ArrayView<const std::vector<double>> &input,
450 const ArrayView<const std::vector<double>> &data,
451 const ArrayView<const std::vector<double>> &velocity,
452 double factor = 1.)
453 {
454 const unsigned int n_dofs = fe.dofs_per_cell;
455 const unsigned int n_comp = fe.get_fe().n_components();
456
459
460 AssertDimension(velocity.size(), dim);
461 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
462 if (v_increment == 1)
463 {
465 }
466
467
468 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
469 {
470 const double dx = factor * fe.JxW(k);
471
472 double nv = 0.;
473 for (unsigned int d = 0; d < dim; ++d)
474 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
475
476 std::vector<double> val(n_comp);
477
478 for (unsigned int d = 0; d < n_comp; ++d)
479 {
480 val[d] = (nv > 0.) ? input[d][k] : -data[d][k];
481 for (unsigned i = 0; i < n_dofs; ++i)
482 {
483 const double v = fe.shape_value_component(i, k, d);
484 result(i) += dx * nv * val[d] * v;
485 }
486 }
487 }
488 }
489
490
491
512 template <int dim>
518 const FEValuesBase<dim> &fe1,
519 const FEValuesBase<dim> &fe2,
520 const FEValuesBase<dim> &fetest1,
521 const FEValuesBase<dim> &fetest2,
522 const ArrayView<const std::vector<double>> &velocity,
523 const double factor = 1.)
524 {
525 const unsigned int n1 = fe1.dofs_per_cell;
526 // Multiply the quadrature point
527 // index below with this factor to
528 // have simpler data for constant
529 // velocities.
530 AssertDimension(velocity.size(), dim);
531 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
532 if (v_increment == 1)
533 {
535 }
536
537 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
538 {
539 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
540 for (unsigned int d = 1; d < dim; ++d)
541 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
542 const double dx_nbeta = factor * std::abs(nbeta) * fe1.JxW(k);
543 FullMatrix<double> &M1 = nbeta > 0. ? M11 : M22;
544 FullMatrix<double> &M2 = nbeta > 0. ? M21 : M12;
545 const FEValuesBase<dim> &fe = nbeta > 0. ? fe1 : fe2;
546 const FEValuesBase<dim> &fetest = nbeta > 0. ? fetest1 : fetest2;
547 const FEValuesBase<dim> &fetestn = nbeta > 0. ? fetest2 : fetest1;
548 for (unsigned i = 0; i < n1; ++i)
549 for (unsigned j = 0; j < n1; ++j)
550 {
551 if (fe1.get_fe().is_primitive())
552 {
553 M1(i, j) += dx_nbeta * fe.shape_value(j, k) *
554 fetest.shape_value(i, k);
555 M2(i, j) -= dx_nbeta * fe.shape_value(j, k) *
556 fetestn.shape_value(i, k);
557 }
558 else
559 {
560 for (unsigned int d = 0; d < fe1.get_fe().n_components();
561 ++d)
562 {
563 M1(i, j) += dx_nbeta *
564 fe.shape_value_component(j, k, d) *
565 fetest.shape_value_component(i, k, d);
566 M2(i, j) -= dx_nbeta *
567 fe.shape_value_component(j, k, d) *
568 fetestn.shape_value_component(i, k, d);
569 }
570 }
571 }
572 }
573 }
574
575
576
597 template <int dim>
600 Vector<double> &result2,
601 const FEValuesBase<dim> &fe1,
602 const FEValuesBase<dim> &fe2,
603 const std::vector<double> &input1,
604 const std::vector<double> &input2,
605 const ArrayView<const std::vector<double>> &velocity,
606 const double factor = 1.)
607 {
608 Assert(fe1.get_fe().n_components() == 1,
610 Assert(fe2.get_fe().n_components() == 1,
612
613 const unsigned int n1 = fe1.dofs_per_cell;
614 // Multiply the quadrature point
615 // index below with this factor to
616 // have simpler data for constant
617 // velocities.
618 AssertDimension(velocity.size(), dim);
619 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
620 if (v_increment == 1)
621 {
623 }
624
625 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
626 {
627 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
628 for (unsigned int d = 1; d < dim; ++d)
629 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
630 const double dx_nbeta = factor * nbeta * fe1.JxW(k);
631
632 for (unsigned i = 0; i < n1; ++i)
633 {
634 const double v1 = fe1.shape_value(i, k);
635 const double v2 = fe2.shape_value(i, k);
636 const double u1 = input1[k];
637 const double u2 = input2[k];
638 if (nbeta > 0)
639 {
640 result1(i) += dx_nbeta * u1 * v1;
641 result2(i) -= dx_nbeta * u1 * v2;
642 }
643 else
644 {
645 result1(i) += dx_nbeta * u2 * v1;
646 result2(i) -= dx_nbeta * u2 * v2;
647 }
648 }
649 }
650 }
651
652
653
674 template <int dim>
677 Vector<double> &result2,
678 const FEValuesBase<dim> &fe1,
679 const FEValuesBase<dim> &fe2,
680 const ArrayView<const std::vector<double>> &input1,
681 const ArrayView<const std::vector<double>> &input2,
682 const ArrayView<const std::vector<double>> &velocity,
683 const double factor = 1.)
684 {
685 const unsigned int n_comp = fe1.get_fe().n_components();
686 const unsigned int n1 = fe1.dofs_per_cell;
689
690 // Multiply the quadrature point
691 // index below with this factor to
692 // have simpler data for constant
693 // velocities.
694 AssertDimension(velocity.size(), dim);
695 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
696 if (v_increment == 1)
697 {
699 }
700
701 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
702 {
703 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
704 for (unsigned int d = 1; d < dim; ++d)
705 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
706 const double dx_nbeta = factor * nbeta * fe1.JxW(k);
707
708 for (unsigned i = 0; i < n1; ++i)
709 for (unsigned int d = 0; d < n_comp; ++d)
710 {
711 const double v1 = fe1.shape_value_component(i, k, d);
712 const double v2 = fe2.shape_value_component(i, k, d);
713 const double u1 = input1[d][k];
714 const double u2 = input2[d][k];
715 if (nbeta > 0)
716 {
717 result1(i) += dx_nbeta * u1 * v1;
718 result2(i) -= dx_nbeta * u1 * v2;
719 }
720 else
721 {
722 result1(i) += dx_nbeta * u2 * v1;
723 result2(i) -= dx_nbeta * u2 * v2;
724 }
725 }
726 }
727 }
728
729 } // namespace Advection
730} // namespace LocalIntegrators
731
732
734
735#endif
const unsigned int dofs_per_cell
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
double shape_value_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const unsigned int n_quadrature_points
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
unsigned int n_components() const
bool is_primitive() const
size_type n() const
size_type m() const
virtual size_type size() const override
#define DEAL_II_DEPRECATED
Definition config.h:294
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
const unsigned int v1
#define Assert(cond, exc)
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< index_type > data
Definition mpi.cc:734
void upwind_face_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const std::vector< double > &input1, const std::vector< double > &input2, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:599
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:72
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:128
void upwind_value_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const std::vector< double > &data, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:380
void upwind_value_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:307
Library of integrals over cells and faces.
Definition advection.h:32
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)