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