Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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\}}\)
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 
24 
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
29 
31 
33 
35 {
50  namespace Advection
51  {
73  template <int dim>
74  void
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  {
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
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  {
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
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 
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  {
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
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  {
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
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 
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  {
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
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  {
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
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  {
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
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 
462 
463  AssertDimension(velocity.size(), dim);
464  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
465  if (v_increment == 1)
466  {
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
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  {
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
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,
613  Assert(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  {
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
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;
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  {
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 
737 
738 #endif
const double & shape_value(const unsigned int i, const unsigned int q_point) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
const unsigned int dofs_per_cell
double shape_value_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const unsigned int n_quadrature_points
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
unsigned int n_components() const
bool is_primitive() const
size_type n() const
size_type m() const
Definition: vector.h:110
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
const unsigned int v1
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1843
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
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
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
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
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
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
Library of integrals over cells and faces.
Definition: advection.h:35
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)