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