Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+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: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2023 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_advection_h
16#define dealii_integrators_advection_h
17
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/fe/mapping.h>
26
28
30
32
34{
49 namespace Advection
50 {
72 template <int dim>
73 void
75 const FEValuesBase<dim> &fe,
76 const FEValuesBase<dim> &fetest,
77 const ArrayView<const std::vector<double>> &velocity,
78 const double factor = 1.)
79 {
80 const unsigned int n_dofs = fe.dofs_per_cell;
81 const unsigned int t_dofs = fetest.dofs_per_cell;
82 const unsigned int n_components = fe.get_fe().n_components();
83
84 AssertDimension(velocity.size(), dim);
85 // If the size of the
86 // velocity vectors is one,
87 // then do not increment
88 // between quadrature points.
89 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
90
91 if (v_increment == 1)
92 {
94 }
95
96 AssertDimension(M.n(), n_dofs);
97 AssertDimension(M.m(), t_dofs);
98
99 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
100 {
101 const double dx = factor * fe.JxW(k);
102 const unsigned int vindex = k * v_increment;
103
104 for (unsigned j = 0; j < n_dofs; ++j)
105 for (unsigned i = 0; i < t_dofs; ++i)
106 for (unsigned int c = 0; c < n_components; ++c)
107 {
108 double wgradv =
109 velocity[0][vindex] * fe.shape_grad_component(i, k, c)[0];
110 for (unsigned int d = 1; d < dim; ++d)
111 wgradv +=
112 velocity[d][vindex] * fe.shape_grad_component(i, k, c)[d];
113 M(i, j) -= dx * wgradv * fe.shape_value_component(j, k, c);
114 }
115 }
116 }
117
118
119
128 template <int dim>
129 inline void
131 const FEValuesBase<dim> &fe,
132 const std::vector<Tensor<1, dim>> &input,
133 const ArrayView<const std::vector<double>> &velocity,
134 double factor = 1.)
135 {
136 const unsigned int nq = fe.n_quadrature_points;
137 const unsigned int n_dofs = fe.dofs_per_cell;
138 Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
139 Assert(result.size() == n_dofs,
140 ExcDimensionMismatch(result.size(), n_dofs));
141
142 AssertDimension(velocity.size(), dim);
143 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
144 if (v_increment == 1)
145 {
147 }
148
149 for (unsigned k = 0; k < nq; ++k)
150 {
151 const double dx = factor * fe.JxW(k);
152 for (unsigned i = 0; i < n_dofs; ++i)
153 for (unsigned int d = 0; d < dim; ++d)
154 result(i) += dx * input[k][d] * fe.shape_value(i, k) *
155 velocity[d][k * v_increment];
156 }
157 }
158
159
160
171 template <int dim>
172 inline void
174 const FEValuesBase<dim> &fe,
175 const ArrayView<const std::vector<Tensor<1, dim>>> &input,
176 const ArrayView<const std::vector<double>> &velocity,
177 double factor = 1.)
178 {
179 const unsigned int nq = fe.n_quadrature_points;
180 const unsigned int n_dofs = fe.dofs_per_cell;
181 const unsigned int n_comp = fe.get_fe().n_components();
182
184 Assert(result.size() == n_dofs,
185 ExcDimensionMismatch(result.size(), n_dofs));
186
187 AssertDimension(velocity.size(), dim);
188 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
189 if (v_increment == 1)
190 {
192 }
193
194 for (unsigned k = 0; k < nq; ++k)
195 {
196 const double dx = factor * fe.JxW(k);
197 for (unsigned i = 0; i < n_dofs; ++i)
198 for (unsigned int c = 0; c < n_comp; ++c)
199 for (unsigned int d = 0; d < dim; ++d)
200 result(i) += dx * input[c][k][d] *
201 fe.shape_value_component(i, k, c) *
202 velocity[d][k * v_increment];
203 }
204 }
205
206
207
213 template <int dim>
214 inline void
216 const FEValuesBase<dim> &fe,
217 const std::vector<double> &input,
218 const ArrayView<const std::vector<double>> &velocity,
219 double factor = 1.)
220 {
221 const unsigned int nq = fe.n_quadrature_points;
222 const unsigned int n_dofs = fe.dofs_per_cell;
223 Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
224 Assert(result.size() == n_dofs,
225 ExcDimensionMismatch(result.size(), n_dofs));
226
227 AssertDimension(velocity.size(), dim);
228 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
229 if (v_increment == 1)
230 {
232 }
233
234 for (unsigned k = 0; k < nq; ++k)
235 {
236 const double dx = factor * fe.JxW(k);
237 for (unsigned i = 0; i < n_dofs; ++i)
238 for (unsigned int d = 0; d < dim; ++d)
239 result(i) -= dx * input[k] * fe.shape_grad(i, k)[d] *
240 velocity[d][k * v_increment];
241 }
242 }
243
244
245
253 template <int dim>
254 inline void
256 const FEValuesBase<dim> &fe,
257 const ArrayView<const std::vector<double>> &input,
258 const ArrayView<const std::vector<double>> &velocity,
259 double factor = 1.)
260 {
261 const unsigned int nq = fe.n_quadrature_points;
262 const unsigned int n_dofs = fe.dofs_per_cell;
263 const unsigned int n_comp = fe.get_fe().n_components();
264
266 Assert(result.size() == n_dofs,
267 ExcDimensionMismatch(result.size(), n_dofs));
268
269 AssertDimension(velocity.size(), dim);
270 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
271 if (v_increment == 1)
272 {
274 }
275
276 for (unsigned k = 0; k < nq; ++k)
277 {
278 const double dx = factor * fe.JxW(k);
279 for (unsigned i = 0; i < n_dofs; ++i)
280 for (unsigned int c = 0; c < n_comp; ++c)
281 for (unsigned int d = 0; d < dim; ++d)
282 result(i) -= dx * input[c][k] *
283 fe.shape_grad_component(i, k, c)[d] *
284 velocity[d][k * v_increment];
285 }
286 }
287
288
289
307 template <int dim>
308 void
310 const FEValuesBase<dim> &fe,
311 const FEValuesBase<dim> &fetest,
312 const ArrayView<const std::vector<double>> &velocity,
313 double factor = 1.)
314 {
315 const unsigned int n_dofs = fe.dofs_per_cell;
316 const unsigned int t_dofs = fetest.dofs_per_cell;
317 unsigned int n_components = fe.get_fe().n_components();
318 AssertDimension(M.m(), n_dofs);
319 AssertDimension(M.n(), n_dofs);
320
321 AssertDimension(velocity.size(), dim);
322 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
323 if (v_increment == 1)
324 {
326 }
327
328 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
329 {
330 const double dx = factor * fe.JxW(k);
331
332 double nv = 0.;
333 for (unsigned int d = 0; d < dim; ++d)
334 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
335
336 if (nv > 0)
337 {
338 for (unsigned i = 0; i < t_dofs; ++i)
339 for (unsigned j = 0; j < n_dofs; ++j)
340 {
341 if (fe.get_fe().is_primitive())
342 M(i, j) +=
343 dx * nv * fe.shape_value(i, k) * fe.shape_value(j, k);
344 else
345 for (unsigned int c = 0; c < n_components; ++c)
346 M(i, j) += dx * nv *
347 fetest.shape_value_component(i, k, c) *
348 fe.shape_value_component(j, k, c);
349 }
350 }
351 }
352 }
353
354
355
380 template <int dim>
381 inline void
383 const FEValuesBase<dim> &fe,
384 const std::vector<double> &input,
385 const std::vector<double> &data,
386 const ArrayView<const std::vector<double>> &velocity,
387 double factor = 1.)
388 {
389 const unsigned int n_dofs = fe.dofs_per_cell;
390
391 AssertDimension(input.size(), fe.n_quadrature_points);
392 AssertDimension(data.size(), fe.n_quadrature_points);
393
394 AssertDimension(velocity.size(), dim);
395 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
396 if (v_increment == 1)
397 {
399 }
400
401
402 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
403 {
404 const double dx = factor * fe.JxW(k);
405
406 double nv = 0.;
407 for (unsigned int d = 0; d < dim; ++d)
408 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
409
410 // Always use the upwind value
411 const double val = (nv > 0.) ? input[k] : -data[k];
412
413 for (unsigned i = 0; i < n_dofs; ++i)
414 {
415 const double v = fe.shape_value(i, k);
416 result(i) += dx * nv * val * v;
417 }
418 }
419 }
420
421
422
447 template <int dim>
448 inline void
450 const FEValuesBase<dim> &fe,
451 const ArrayView<const std::vector<double>> &input,
452 const ArrayView<const std::vector<double>> &data,
453 const ArrayView<const std::vector<double>> &velocity,
454 double factor = 1.)
455 {
456 const unsigned int n_dofs = fe.dofs_per_cell;
457 const unsigned int n_comp = fe.get_fe().n_components();
458
461
462 AssertDimension(velocity.size(), dim);
463 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
464 if (v_increment == 1)
465 {
467 }
468
469
470 for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
471 {
472 const double dx = factor * fe.JxW(k);
473
474 double nv = 0.;
475 for (unsigned int d = 0; d < dim; ++d)
476 nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
477
478 std::vector<double> val(n_comp);
479
480 for (unsigned int d = 0; d < n_comp; ++d)
481 {
482 val[d] = (nv > 0.) ? input[d][k] : -data[d][k];
483 for (unsigned i = 0; i < n_dofs; ++i)
484 {
485 const double v = fe.shape_value_component(i, k, d);
486 result(i) += dx * nv * val[d] * v;
487 }
488 }
489 }
490 }
491
492
493
514 template <int dim>
515 void
520 const FEValuesBase<dim> &fe1,
521 const FEValuesBase<dim> &fe2,
522 const FEValuesBase<dim> &fetest1,
523 const FEValuesBase<dim> &fetest2,
524 const ArrayView<const std::vector<double>> &velocity,
525 const double factor = 1.)
526 {
527 const unsigned int n1 = fe1.dofs_per_cell;
528 // Multiply the quadrature point
529 // index below with this factor to
530 // have simpler data for constant
531 // velocities.
532 AssertDimension(velocity.size(), dim);
533 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
534 if (v_increment == 1)
535 {
537 }
538
539 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
540 {
541 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
542 for (unsigned int d = 1; d < dim; ++d)
543 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
544 const double dx_nbeta = factor * std::abs(nbeta) * fe1.JxW(k);
545 FullMatrix<double> &M1 = nbeta > 0. ? M11 : M22;
546 FullMatrix<double> &M2 = nbeta > 0. ? M21 : M12;
547 const FEValuesBase<dim> &fe = nbeta > 0. ? fe1 : fe2;
548 const FEValuesBase<dim> &fetest = nbeta > 0. ? fetest1 : fetest2;
549 const FEValuesBase<dim> &fetestn = nbeta > 0. ? fetest2 : fetest1;
550 for (unsigned i = 0; i < n1; ++i)
551 for (unsigned j = 0; j < n1; ++j)
552 {
553 if (fe1.get_fe().is_primitive())
554 {
555 M1(i, j) += dx_nbeta * fe.shape_value(j, k) *
556 fetest.shape_value(i, k);
557 M2(i, j) -= dx_nbeta * fe.shape_value(j, k) *
558 fetestn.shape_value(i, k);
559 }
560 else
561 {
562 for (unsigned int d = 0; d < fe1.get_fe().n_components();
563 ++d)
564 {
565 M1(i, j) += dx_nbeta *
566 fe.shape_value_component(j, k, d) *
567 fetest.shape_value_component(i, k, d);
568 M2(i, j) -= dx_nbeta *
569 fe.shape_value_component(j, k, d) *
570 fetestn.shape_value_component(i, k, d);
571 }
572 }
573 }
574 }
575 }
576
577
578
599 template <int dim>
600 void
602 Vector<double> &result2,
603 const FEValuesBase<dim> &fe1,
604 const FEValuesBase<dim> &fe2,
605 const std::vector<double> &input1,
606 const std::vector<double> &input2,
607 const ArrayView<const std::vector<double>> &velocity,
608 const double factor = 1.)
609 {
610 Assert(fe1.get_fe().n_components() == 1,
612 Assert(fe2.get_fe().n_components() == 1,
614
615 const unsigned int n1 = fe1.dofs_per_cell;
616 // Multiply the quadrature point
617 // index below with this factor to
618 // have simpler data for constant
619 // velocities.
620 AssertDimension(velocity.size(), dim);
621 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
622 if (v_increment == 1)
623 {
625 }
626
627 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
628 {
629 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
630 for (unsigned int d = 1; d < dim; ++d)
631 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
632 const double dx_nbeta = factor * nbeta * fe1.JxW(k);
633
634 for (unsigned i = 0; i < n1; ++i)
635 {
636 const double v1 = fe1.shape_value(i, k);
637 const double v2 = fe2.shape_value(i, k);
638 const double u1 = input1[k];
639 const double u2 = input2[k];
640 if (nbeta > 0)
641 {
642 result1(i) += dx_nbeta * u1 * v1;
643 result2(i) -= dx_nbeta * u1 * v2;
644 }
645 else
646 {
647 result1(i) += dx_nbeta * u2 * v1;
648 result2(i) -= dx_nbeta * u2 * v2;
649 }
650 }
651 }
652 }
653
654
655
676 template <int dim>
677 void
679 Vector<double> &result2,
680 const FEValuesBase<dim> &fe1,
681 const FEValuesBase<dim> &fe2,
682 const ArrayView<const std::vector<double>> &input1,
683 const ArrayView<const std::vector<double>> &input2,
684 const ArrayView<const std::vector<double>> &velocity,
685 const double factor = 1.)
686 {
687 const unsigned int n_comp = fe1.get_fe().n_components();
688 const unsigned int n1 = fe1.dofs_per_cell;
691
692 // Multiply the quadrature point
693 // index below with this factor to
694 // have simpler data for constant
695 // velocities.
696 AssertDimension(velocity.size(), dim);
697 const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
698 if (v_increment == 1)
699 {
701 }
702
703 for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
704 {
705 double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
706 for (unsigned int d = 1; d < dim; ++d)
707 nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
708 const double dx_nbeta = factor * nbeta * fe1.JxW(k);
709
710 for (unsigned i = 0; i < n1; ++i)
711 for (unsigned int d = 0; d < n_comp; ++d)
712 {
713 const double v1 = fe1.shape_value_component(i, k, d);
714 const double v2 = fe2.shape_value_component(i, k, d);
715 const double u1 = input1[d][k];
716 const double u2 = input2[d][k];
717 if (nbeta > 0)
718 {
719 result1(i) += dx_nbeta * u1 * v1;
720 result2(i) -= dx_nbeta * u1 * v2;
721 }
722 else
723 {
724 result1(i) += dx_nbeta * u2 * v1;
725 result2(i) -= dx_nbeta * u2 * v2;
726 }
727 }
728 }
729 }
730
731 } // namespace Advection
732} // namespace LocalIntegrators
733
734
736
737#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_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
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:601
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:74
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:130
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:382
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:309
Library of integrals over cells and faces.
Definition advection.h:34
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)