Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
advection.h
1 // ---------------------------------------------------------------------
2 // @f$Id: advection.h 30942 2013-09-26 08:58:02Z kanschat @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__integrators_advection_h
18 #define __deal2__integrators_advection_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 namespace LocalIntegrators
32 {
48  namespace Advection
49  {
76  template<int dim>
77  void cell_matrix (
79  const FEValuesBase<dim> &fe,
80  const FEValuesBase<dim> &fetest,
81  const VectorSlice<const std::vector<std::vector<double> > > &velocity,
82  const double factor = 1.)
83  {
84  const unsigned int n_dofs = fe.dofs_per_cell;
85  const unsigned int t_dofs = fetest.dofs_per_cell;
86  const unsigned int n_components = fe.get_fe().n_components();
87 
88  AssertDimension(velocity.size(), dim);
89  // If the size of the
90  // velocity vectors is one,
91  // then do not increment
92  // between quadrature points.
93  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
94 
95  if (v_increment == 1)
96  {
98  }
99 
100  AssertDimension(M.n(), n_dofs);
101  AssertDimension(M.m(), t_dofs);
102 
103  for (unsigned k=0; k<fe.n_quadrature_points; ++k)
104  {
105  const double dx = factor * fe.JxW(k);
106  const unsigned int vindex = k * v_increment;
107 
108  for (unsigned j=0; j<n_dofs; ++j)
109  for (unsigned i=0; i<t_dofs; ++i)
110  for (unsigned int c=0; c<n_components; ++c)
111  {
112  double wgradv = velocity[0][vindex]
113  * fe.shape_grad_component(i,k,c)[0];
114  for (unsigned int d=1; d<dim; ++d)
115  wgradv += velocity[d][vindex]
116  * fe.shape_grad_component(i,k,c)[d];
117  M(i,j) -= dx * wgradv * fe.shape_value_component(j,k,c);
118  }
119  }
120  }
121 
129  template <int dim>
130  inline void
132  Vector<double> &result,
133  const FEValuesBase<dim> &fe,
134  const std::vector<Tensor<1,dim> > &input,
135  const VectorSlice<const std::vector<std::vector<double> > > &velocity,
136  double factor = 1.)
137  {
138  const unsigned int nq = fe.n_quadrature_points;
139  const unsigned int n_dofs = fe.dofs_per_cell;
140  Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
141  Assert(result.size() == n_dofs, 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]
156  * fe.shape_value(i,k) * velocity[d][k * v_increment];
157  }
158  }
159 
160 
169  template <int dim>
170  inline void
172  Vector<double> &result,
173  const FEValuesBase<dim> &fe,
174  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
175  const VectorSlice<const std::vector<std::vector<double> > > &velocity,
176  double factor = 1.)
177  {
178  const unsigned int nq = fe.n_quadrature_points;
179  const unsigned int n_dofs = fe.dofs_per_cell;
180  const unsigned int n_comp = fe.get_fe().n_components();
181 
183  Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
184 
185  AssertDimension(velocity.size(), dim);
186  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
187  if (v_increment == 1)
188  {
190  }
191 
192  for (unsigned k=0; k<nq; ++k)
193  {
194  const double dx = factor * fe.JxW(k);
195  for (unsigned i=0; i<n_dofs; ++i)
196  for (unsigned int c=0; c<n_comp; ++c)
197  for (unsigned int d=0; d<dim; ++d)
198  result(i) += dx * input[c][k][d]
199  * fe.shape_value_component(i,k,c) * velocity[d][k * v_increment];
200  }
201  }
202 
232  template <int dim>
235  const FEValuesBase<dim> &fe,
236  const FEValuesBase<dim> &fetest,
237  const VectorSlice<const std::vector<std::vector<double> > > &velocity,
238  double factor = 1.)
239  {
240  const unsigned int n_dofs = fe.dofs_per_cell;
241  const unsigned int t_dofs = fetest.dofs_per_cell;
242  unsigned int n_components = fe.get_fe().n_components();
243  AssertDimension (M.m(), n_dofs);
244  AssertDimension (M.n(), n_dofs);
245 
246  AssertDimension(velocity.size(), dim);
247  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
248  if (v_increment == 1)
249  {
251  }
252 
253  for (unsigned k=0; k<fe.n_quadrature_points; ++k)
254  {
255  const double dx = factor * fe.JxW(k);
256 
257  double nv = 0.;
258  for (unsigned int d=0; d<dim; ++d)
259  nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
260 
261  if (nv > 0)
262  {
263  for (unsigned i=0; i<t_dofs; ++i)
264  for (unsigned j=0; j<n_dofs; ++j)
265  {
266  if (fe.get_fe().is_primitive())
267  M(i,j) += dx * nv * fe.shape_value(i,k) * fe.shape_value(j,k);
268  else
269  for (unsigned int c=0; c<n_components; ++c)
270  M(i,j) += dx * nv * fetest.shape_value_component(i,k,c)
271  * fe.shape_value_component(j,k,c);
272  }
273  }
274  }
275  }
276 
309  template <int dim>
311  FullMatrix<double> &M11,
312  FullMatrix<double> &M12,
313  FullMatrix<double> &M21,
314  FullMatrix<double> &M22,
315  const FEValuesBase<dim> &fe1,
316  const FEValuesBase<dim> &fe2,
317  const FEValuesBase<dim> &fetest1,
318  const FEValuesBase<dim> &fetest2,
319  const VectorSlice<const std::vector<std::vector<double> > > &velocity,
320  const double factor = 1.)
321  {
322  const unsigned int n1 = fe1.dofs_per_cell;
323  // Multiply the quadrature point
324  // index below with this factor to
325  // have simpler data for constant
326  // velocities.
327  AssertDimension(velocity.size(), dim);
328  const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
329  if (v_increment == 1)
330  {
332  }
333 
334  for (unsigned k=0; k<fe1.n_quadrature_points; ++k)
335  {
336  double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
337  for (unsigned int d=1; d<dim; ++d)
338  nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
339  const double dx_nbeta = factor * nbeta * fe1.JxW(k);
340 
341  for (unsigned i=0; i<n1; ++i)
342  for (unsigned j=0; j<n1; ++j)
343  if (fe1.get_fe().is_primitive())
344  {
345  if (nbeta > 0)
346  {
347  M11(i,j) += dx_nbeta
348  * fe1.shape_value(j,k)
349  * fetest1.shape_value(i,k);
350  M21(i,j) -= dx_nbeta
351  * fe1.shape_value(j,k)
352  * fetest2.shape_value(i,k);
353  }
354  else
355  {
356  M22(i,j) -= dx_nbeta
357  * fe2.shape_value(j,k)
358  * fetest2.shape_value(i,k);
359  M12(i,j) += dx_nbeta
360  * fe2.shape_value(j,k)
361  * fetest1.shape_value(i,k);
362  }
363  }
364  else
365  {
366  for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
367  if (nbeta > 0)
368  {
369  M11(i,j) += dx_nbeta
370  * fe1.shape_value_component(j,k,d)
371  * fetest1.shape_value_component(i,k,d);
372  M21(i,j) -= dx_nbeta
373  * fe1.shape_value_component(j,k,d)
374  * fetest2.shape_value_component(i,k,d);
375  }
376  else
377  {
378  M22(i,j) -= dx_nbeta
379  * fe2.shape_value_component(j,k,d)
380  * fetest2.shape_value_component(i,k,d);
381  M12(i,j) += dx_nbeta
382  * fe2.shape_value_component(j,k,d)
383  * fetest1.shape_value_component(i,k,d);
384  }
385  }
386  }
387  }
388  }
389 }
390 
391 
392 DEAL_II_NAMESPACE_CLOSE
393 
394 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:870
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
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:131
const FiniteElement< dim, spacedim > & get_fe() const
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:299
bool is_primitive(const unsigned int i) const
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Point< spacedim > & normal_vector(const unsigned int i) const
unsigned int n_components() const
size_type m() const
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:233
double JxW(const unsigned int quadrature_point) const
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:77
std::size_t size() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const