Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
advection.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 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 
32 DEAL_II_NAMESPACE_OPEN
33 
34 namespace LocalIntegrators
35 {
52  namespace Advection
53  {
78  template <int dim>
79  void
80  cell_matrix(FullMatrix<double> & M,
81  const FEValuesBase<dim> & fe,
82  const FEValuesBase<dim> & fetest,
83  const ArrayView<const std::vector<double>> &velocity,
84  const double factor = 1.)
85  {
86  const unsigned int n_dofs = fe.dofs_per_cell;
87  const unsigned int t_dofs = fetest.dofs_per_cell;
88  const unsigned int n_components = fe.get_fe().n_components();
89 
90  AssertDimension(velocity.size(), dim);
91  // If the size of the
92  // velocity vectors is one,
93  // then do not increment
94  // between quadrature points.
95  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
96 
97  if (v_increment == 1)
98  {
99  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
100  }
101 
102  AssertDimension(M.n(), n_dofs);
103  AssertDimension(M.m(), t_dofs);
104 
105  for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
106  {
107  const double dx = factor * fe.JxW(k);
108  const unsigned int vindex = k * v_increment;
109 
110  for (unsigned j = 0; j < n_dofs; ++j)
111  for (unsigned i = 0; i < t_dofs; ++i)
112  for (unsigned int c = 0; c < n_components; ++c)
113  {
114  double wgradv =
115  velocity[0][vindex] * fe.shape_grad_component(i, k, c)[0];
116  for (unsigned int d = 1; d < dim; ++d)
117  wgradv +=
118  velocity[d][vindex] * fe.shape_grad_component(i, k, c)[d];
119  M(i, j) -= dx * wgradv * fe.shape_value_component(j, k, c);
120  }
121  }
122  }
123 
124 
125 
134  template <int dim>
135  inline void
136  cell_residual(Vector<double> & result,
137  const FEValuesBase<dim> & fe,
138  const std::vector<Tensor<1, dim>> & input,
139  const ArrayView<const std::vector<double>> &velocity,
140  double factor = 1.)
141  {
142  const unsigned int nq = fe.n_quadrature_points;
143  const unsigned int n_dofs = fe.dofs_per_cell;
144  Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
145  Assert(result.size() == n_dofs,
146  ExcDimensionMismatch(result.size(), n_dofs));
147 
148  AssertDimension(velocity.size(), dim);
149  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
150  if (v_increment == 1)
151  {
152  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
153  }
154 
155  for (unsigned k = 0; k < nq; ++k)
156  {
157  const double dx = factor * fe.JxW(k);
158  for (unsigned i = 0; i < n_dofs; ++i)
159  for (unsigned int d = 0; d < dim; ++d)
160  result(i) += dx * input[k][d] * fe.shape_value(i, k) *
161  velocity[d][k * v_increment];
162  }
163  }
164 
165 
166 
177  template <int dim>
178  inline void
179  cell_residual(Vector<double> & result,
180  const FEValuesBase<dim> & fe,
181  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
182  const ArrayView<const std::vector<double>> & velocity,
183  double factor = 1.)
184  {
185  const unsigned int nq = fe.n_quadrature_points;
186  const unsigned int n_dofs = fe.dofs_per_cell;
187  const unsigned int n_comp = fe.get_fe().n_components();
188 
189  AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
190  Assert(result.size() == n_dofs,
191  ExcDimensionMismatch(result.size(), n_dofs));
192 
193  AssertDimension(velocity.size(), dim);
194  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
195  if (v_increment == 1)
196  {
197  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
198  }
199 
200  for (unsigned k = 0; k < nq; ++k)
201  {
202  const double dx = factor * fe.JxW(k);
203  for (unsigned i = 0; i < n_dofs; ++i)
204  for (unsigned int c = 0; c < n_comp; ++c)
205  for (unsigned int d = 0; d < dim; ++d)
206  result(i) += dx * input[c][k][d] *
207  fe.shape_value_component(i, k, c) *
208  velocity[d][k * v_increment];
209  }
210  }
211 
212 
213 
219  template <int dim>
220  inline void
221  cell_residual(Vector<double> & result,
222  const FEValuesBase<dim> & fe,
223  const std::vector<double> & input,
224  const ArrayView<const std::vector<double>> &velocity,
225  double factor = 1.)
226  {
227  const unsigned int nq = fe.n_quadrature_points;
228  const unsigned int n_dofs = fe.dofs_per_cell;
229  Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
230  Assert(result.size() == n_dofs,
231  ExcDimensionMismatch(result.size(), n_dofs));
232 
233  AssertDimension(velocity.size(), dim);
234  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
235  if (v_increment == 1)
236  {
237  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
238  }
239 
240  for (unsigned k = 0; k < nq; ++k)
241  {
242  const double dx = factor * fe.JxW(k);
243  for (unsigned i = 0; i < n_dofs; ++i)
244  for (unsigned int d = 0; d < dim; ++d)
245  result(i) -= dx * input[k] * fe.shape_grad(i, k)[d] *
246  velocity[d][k * v_increment];
247  }
248  }
249 
250 
251 
259  template <int dim>
260  inline void
261  cell_residual(Vector<double> & result,
262  const FEValuesBase<dim> & fe,
263  const ArrayView<const std::vector<double>> &input,
264  const ArrayView<const std::vector<double>> &velocity,
265  double factor = 1.)
266  {
267  const unsigned int nq = fe.n_quadrature_points;
268  const unsigned int n_dofs = fe.dofs_per_cell;
269  const unsigned int n_comp = fe.get_fe().n_components();
270 
271  AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
272  Assert(result.size() == n_dofs,
273  ExcDimensionMismatch(result.size(), n_dofs));
274 
275  AssertDimension(velocity.size(), dim);
276  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
277  if (v_increment == 1)
278  {
279  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
280  }
281 
282  for (unsigned k = 0; k < nq; ++k)
283  {
284  const double dx = factor * fe.JxW(k);
285  for (unsigned i = 0; i < n_dofs; ++i)
286  for (unsigned int c = 0; c < n_comp; ++c)
287  for (unsigned int d = 0; d < dim; ++d)
288  result(i) -= dx * input[c][k] *
289  fe.shape_grad_component(i, k, c)[d] *
290  velocity[d][k * v_increment];
291  }
292  }
293 
294 
295 
313  template <int dim>
314  void
315  upwind_value_matrix(FullMatrix<double> & M,
316  const FEValuesBase<dim> & fe,
317  const FEValuesBase<dim> & fetest,
318  const ArrayView<const std::vector<double>> &velocity,
319  double factor = 1.)
320  {
321  const unsigned int n_dofs = fe.dofs_per_cell;
322  const unsigned int t_dofs = fetest.dofs_per_cell;
323  unsigned int n_components = fe.get_fe().n_components();
324  AssertDimension(M.m(), n_dofs);
325  AssertDimension(M.n(), n_dofs);
326 
327  AssertDimension(velocity.size(), dim);
328  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
329  if (v_increment == 1)
330  {
331  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
332  }
333 
334  for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
335  {
336  const double dx = factor * fe.JxW(k);
337 
338  double nv = 0.;
339  for (unsigned int d = 0; d < dim; ++d)
340  nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
341 
342  if (nv > 0)
343  {
344  for (unsigned i = 0; i < t_dofs; ++i)
345  for (unsigned j = 0; j < n_dofs; ++j)
346  {
347  if (fe.get_fe().is_primitive())
348  M(i, j) +=
349  dx * nv * fe.shape_value(i, k) * fe.shape_value(j, k);
350  else
351  for (unsigned int c = 0; c < n_components; ++c)
352  M(i, j) += dx * nv *
353  fetest.shape_value_component(i, k, c) *
354  fe.shape_value_component(j, k, c);
355  }
356  }
357  }
358  }
359 
360 
361 
386  template <int dim>
387  inline void
388  upwind_value_residual(Vector<double> & result,
389  const FEValuesBase<dim> & fe,
390  const std::vector<double> & input,
391  const std::vector<double> & data,
392  const ArrayView<const std::vector<double>> &velocity,
393  double factor = 1.)
394  {
395  const unsigned int n_dofs = fe.dofs_per_cell;
396 
397  AssertDimension(input.size(), fe.n_quadrature_points);
398  AssertDimension(data.size(), fe.n_quadrature_points);
399 
400  AssertDimension(velocity.size(), dim);
401  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
402  if (v_increment == 1)
403  {
404  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
405  }
406 
407 
408  for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
409  {
410  const double dx = factor * fe.JxW(k);
411 
412  double nv = 0.;
413  for (unsigned int d = 0; d < dim; ++d)
414  nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
415 
416  // Always use the upwind value
417  const double val = (nv > 0.) ? input[k] : -data[k];
418 
419  for (unsigned i = 0; i < n_dofs; ++i)
420  {
421  const double v = fe.shape_value(i, k);
422  result(i) += dx * nv * val * v;
423  }
424  }
425  }
426 
427 
428 
453  template <int dim>
454  inline void
455  upwind_value_residual(Vector<double> & result,
456  const FEValuesBase<dim> & fe,
457  const ArrayView<const std::vector<double>> &input,
458  const ArrayView<const std::vector<double>> &data,
459  const ArrayView<const std::vector<double>> &velocity,
460  double factor = 1.)
461  {
462  const unsigned int n_dofs = fe.dofs_per_cell;
463  const unsigned int n_comp = fe.get_fe().n_components();
464 
465  AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
466  AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);
467 
468  AssertDimension(velocity.size(), dim);
469  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
470  if (v_increment == 1)
471  {
472  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
473  }
474 
475 
476  for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
477  {
478  const double dx = factor * fe.JxW(k);
479 
480  double nv = 0.;
481  for (unsigned int d = 0; d < dim; ++d)
482  nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
483 
484  std::vector<double> val(n_comp);
485 
486  for (unsigned int d = 0; d < n_comp; ++d)
487  {
488  val[d] = (nv > 0.) ? input[d][k] : -data[d][k];
489  for (unsigned i = 0; i < n_dofs; ++i)
490  {
491  const double v = fe.shape_value_component(i, k, d);
492  result(i) += dx * nv * val[d] * v;
493  }
494  }
495  }
496  }
497 
498 
499 
520  template <int dim>
521  void
522  upwind_value_matrix(FullMatrix<double> & M11,
523  FullMatrix<double> & M12,
524  FullMatrix<double> & M21,
525  FullMatrix<double> & M22,
526  const FEValuesBase<dim> & fe1,
527  const FEValuesBase<dim> & fe2,
528  const FEValuesBase<dim> & fetest1,
529  const FEValuesBase<dim> & fetest2,
530  const ArrayView<const std::vector<double>> &velocity,
531  const double factor = 1.)
532  {
533  const unsigned int n1 = fe1.dofs_per_cell;
534  // Multiply the quadrature point
535  // index below with this factor to
536  // have simpler data for constant
537  // velocities.
538  AssertDimension(velocity.size(), dim);
539  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
540  if (v_increment == 1)
541  {
542  AssertVectorVectorDimension(velocity, dim, fe1.n_quadrature_points);
543  }
544 
545  for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
546  {
547  double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
548  for (unsigned int d = 1; d < dim; ++d)
549  nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
550  const double dx_nbeta = factor * std::abs(nbeta) * fe1.JxW(k);
551  FullMatrix<double> &M1 = nbeta > 0. ? M11 : M22;
552  FullMatrix<double> &M2 = nbeta > 0. ? M21 : M12;
553  const FEValuesBase<dim> &fe = nbeta > 0. ? fe1 : fe2;
554  const FEValuesBase<dim> &fetest = nbeta > 0. ? fetest1 : fetest2;
555  const FEValuesBase<dim> &fetestn = nbeta > 0. ? fetest2 : fetest1;
556  for (unsigned i = 0; i < n1; ++i)
557  for (unsigned j = 0; j < n1; ++j)
558  {
559  if (fe1.get_fe().is_primitive())
560  {
561  M1(i, j) += dx_nbeta * fe.shape_value(j, k) *
562  fetest.shape_value(i, k);
563  M2(i, j) -= dx_nbeta * fe.shape_value(j, k) *
564  fetestn.shape_value(i, k);
565  }
566  else
567  {
568  for (unsigned int d = 0; d < fe1.get_fe().n_components();
569  ++d)
570  {
571  M1(i, j) += dx_nbeta *
572  fe.shape_value_component(j, k, d) *
573  fetest.shape_value_component(i, k, d);
574  M2(i, j) -= dx_nbeta *
575  fe.shape_value_component(j, k, d) *
576  fetestn.shape_value_component(i, k, d);
577  }
578  }
579  }
580  }
581  }
582 
583 
584 
605  template <int dim>
606  void
607  upwind_face_residual(Vector<double> & result1,
608  Vector<double> & result2,
609  const FEValuesBase<dim> & fe1,
610  const FEValuesBase<dim> & fe2,
611  const std::vector<double> & input1,
612  const std::vector<double> & input2,
613  const ArrayView<const std::vector<double>> &velocity,
614  const double factor = 1.)
615  {
616  Assert(fe1.get_fe().n_components() == 1,
617  ExcDimensionMismatch(fe1.get_fe().n_components(), 1));
618  Assert(fe2.get_fe().n_components() == 1,
619  ExcDimensionMismatch(fe2.get_fe().n_components(), 1));
620 
621  const unsigned int n1 = fe1.dofs_per_cell;
622  // Multiply the quadrature point
623  // index below with this factor to
624  // have simpler data for constant
625  // velocities.
626  AssertDimension(velocity.size(), dim);
627  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
628  if (v_increment == 1)
629  {
630  AssertVectorVectorDimension(velocity, dim, fe1.n_quadrature_points);
631  }
632 
633  for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
634  {
635  double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
636  for (unsigned int d = 1; d < dim; ++d)
637  nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
638  const double dx_nbeta = factor * nbeta * fe1.JxW(k);
639 
640  for (unsigned i = 0; i < n1; ++i)
641  {
642  const double v1 = fe1.shape_value(i, k);
643  const double v2 = fe2.shape_value(i, k);
644  const double u1 = input1[k];
645  const double u2 = input2[k];
646  if (nbeta > 0)
647  {
648  result1(i) += dx_nbeta * u1 * v1;
649  result2(i) -= dx_nbeta * u1 * v2;
650  }
651  else
652  {
653  result1(i) += dx_nbeta * u2 * v1;
654  result2(i) -= dx_nbeta * u2 * v2;
655  }
656  }
657  }
658  }
659 
660 
661 
682  template <int dim>
683  void
684  upwind_face_residual(Vector<double> & result1,
685  Vector<double> & result2,
686  const FEValuesBase<dim> & fe1,
687  const FEValuesBase<dim> & fe2,
688  const ArrayView<const std::vector<double>> &input1,
689  const ArrayView<const std::vector<double>> &input2,
690  const ArrayView<const std::vector<double>> &velocity,
691  const double factor = 1.)
692  {
693  const unsigned int n_comp = fe1.get_fe().n_components();
694  const unsigned int n1 = fe1.dofs_per_cell;
695  AssertVectorVectorDimension(input1, n_comp, fe1.n_quadrature_points);
696  AssertVectorVectorDimension(input2, n_comp, fe2.n_quadrature_points);
697 
698  // Multiply the quadrature point
699  // index below with this factor to
700  // have simpler data for constant
701  // velocities.
702  AssertDimension(velocity.size(), dim);
703  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
704  if (v_increment == 1)
705  {
706  AssertVectorVectorDimension(velocity, dim, fe1.n_quadrature_points);
707  }
708 
709  for (unsigned k = 0; k < fe1.n_quadrature_points; ++k)
710  {
711  double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
712  for (unsigned int d = 1; d < dim; ++d)
713  nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
714  const double dx_nbeta = factor * nbeta * fe1.JxW(k);
715 
716  for (unsigned i = 0; i < n1; ++i)
717  for (unsigned int d = 0; d < n_comp; ++d)
718  {
719  const double v1 = fe1.shape_value_component(i, k, d);
720  const double v2 = fe2.shape_value_component(i, k, d);
721  const double u1 = input1[d][k];
722  const double u2 = input2[d][k];
723  if (nbeta > 0)
724  {
725  result1(i) += dx_nbeta * u1 * v1;
726  result2(i) -= dx_nbeta * u1 * v2;
727  }
728  else
729  {
730  result1(i) += dx_nbeta * u2 * v1;
731  result2(i) -= dx_nbeta * u2 * v2;
732  }
733  }
734  }
735  }
736 
737  } // namespace Advection
738 } // namespace LocalIntegrators
739 
740 
741 DEAL_II_NAMESPACE_CLOSE
742 
743 #endif
FullMatrix::m
size_type m() const
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
FEValuesBase::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_values.h:2106
FEValuesBase::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1595
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:80
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:136
FEValuesBase::shape_value_component
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
LocalIntegrators
Library of integrals over cells and faces.
Definition: advection.h:34
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:315
FullMatrix::n
size_type n() const
FEValuesBase::shape_value
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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:607
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:388
FEValuesBase::n_quadrature_points
const unsigned int n_quadrature_points
Definition: fe_values.h:2099
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::JxW
double JxW(const unsigned int quadrature_point) const
Vector::size
size_type size() const
FEValuesBase::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
FEValuesBase::normal_vector
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Generated by   doxygen 1.8.14