Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
grad_div.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_grad_div_h
17 #define dealii_integrators_grad_div_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 {
44  namespace GradDiv
45  {
55  template <int dim>
56  void
58  const FEValuesBase<dim> &fe,
59  double factor = 1.)
60  {
61  const unsigned int n_dofs = fe.dofs_per_cell;
62 
63  AssertDimension(fe.get_fe().n_components(), dim);
64  AssertDimension(M.m(), n_dofs);
65  AssertDimension(M.n(), n_dofs);
66 
67  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
68  {
69  const double dx = factor * fe.JxW(k);
70  for (unsigned int i = 0; i < n_dofs; ++i)
71  for (unsigned int j = 0; j < n_dofs; ++j)
72  {
73  const double divu =
74  fe[FEValuesExtractors::Vector(0)].divergence(j, k);
75  const double divv =
76  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
77 
78  M(i, j) += dx * divu * divv;
79  }
80  }
81  }
82 
92  template <int dim, typename number>
93  void
95  const FEValuesBase<dim> & fetest,
96  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
97  const double factor = 1.)
98  {
99  const unsigned int n_dofs = fetest.dofs_per_cell;
100 
101  AssertDimension(fetest.get_fe().n_components(), dim);
103 
104  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
105  {
106  const double dx = factor * fetest.JxW(k);
107  for (unsigned int i = 0; i < n_dofs; ++i)
108  {
109  const double divv =
110  fetest[FEValuesExtractors::Vector(0)].divergence(i, k);
111  double du = 0.;
112  for (unsigned int d = 0; d < dim; ++d)
113  du += input[d][k][d];
114 
115  result(i) += dx * du * divv;
116  }
117  }
118  }
119 
128  template <int dim>
129  inline void
131  const FEValuesBase<dim> &fe,
132  double penalty,
133  double factor = 1.)
134  {
135  const unsigned int n_dofs = fe.dofs_per_cell;
136 
137  AssertDimension(fe.get_fe().n_components(), dim);
138  AssertDimension(M.m(), n_dofs);
139  AssertDimension(M.n(), n_dofs);
140 
141  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
142  {
143  const double dx = factor * fe.JxW(k);
144  const Tensor<1, dim> n = fe.normal_vector(k);
145  for (unsigned int i = 0; i < n_dofs; ++i)
146  for (unsigned int j = 0; j < n_dofs; ++j)
147  {
148  const double divu =
149  fe[FEValuesExtractors::Vector(0)].divergence(j, k);
150  const double divv =
151  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
152  double un = 0., vn = 0.;
153  for (unsigned int d = 0; d < dim; ++d)
154  {
155  un += fe.shape_value_component(j, k, d) * n[d];
156  vn += fe.shape_value_component(i, k, d) * n[d];
157  }
158 
159  M(i, j) += dx * 2. * penalty * un * vn;
160  M(i, j) -= dx * (divu * vn + divv * un);
161  }
162  }
163  }
164 
183  template <int dim>
184  void
186  const FEValuesBase<dim> & fe,
187  const ArrayView<const std::vector<double>> & input,
188  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
189  const ArrayView<const std::vector<double>> & data,
190  double penalty,
191  double factor = 1.)
192  {
193  const unsigned int n_dofs = fe.dofs_per_cell;
194  AssertDimension(fe.get_fe().n_components(), dim)
198 
199  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
200  {
201  const double dx = factor * fe.JxW(k);
202  const Tensor<1, dim> n = fe.normal_vector(k);
203 
204  double umgn = 0.;
205  double divu = 0.;
206  for (unsigned int d = 0; d < dim; ++d)
207  {
208  umgn += (input[d][k] - data[d][k]) * n[d];
209  divu += Dinput[d][k][d];
210  }
211 
212  for (unsigned int i = 0; i < n_dofs; ++i)
213  {
214  double vn = 0.;
215  const double divv =
216  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
217  for (unsigned int d = 0; d < dim; ++d)
218  vn += fe.shape_value_component(i, k, d) * n[d];
219 
220  result(i) +=
221  dx * (2. * penalty * umgn * vn - divv * umgn - divu * vn);
222  }
223  }
224  }
225 
234  template <int dim>
235  void
237  FullMatrix<double> & M12,
238  FullMatrix<double> & M21,
239  FullMatrix<double> & M22,
240  const FEValuesBase<dim> &fe1,
241  const FEValuesBase<dim> &fe2,
242  double penalty,
243  double factor1 = 1.,
244  double factor2 = -1.)
245  {
246  const unsigned int n_dofs = fe1.dofs_per_cell;
247  AssertDimension(M11.n(), n_dofs);
248  AssertDimension(M11.m(), n_dofs);
249  AssertDimension(M12.n(), n_dofs);
250  AssertDimension(M12.m(), n_dofs);
251  AssertDimension(M21.n(), n_dofs);
252  AssertDimension(M21.m(), n_dofs);
253  AssertDimension(M22.n(), n_dofs);
254  AssertDimension(M22.m(), n_dofs);
255 
256  const double fi = factor1;
257  const double fe = (factor2 < 0) ? factor1 : factor2;
258  const double f = .5 * (fi + fe);
259 
260  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
261  {
262  const double dx = fe1.JxW(k);
263  const Tensor<1, dim> n = fe1.normal_vector(k);
264  for (unsigned int i = 0; i < n_dofs; ++i)
265  for (unsigned int j = 0; j < n_dofs; ++j)
266  {
267  double uni = 0.;
268  double une = 0.;
269  double vni = 0.;
270  double vne = 0.;
271  const double divui =
272  fe1[FEValuesExtractors::Vector(0)].divergence(j, k);
273  const double divue =
274  fe2[FEValuesExtractors::Vector(0)].divergence(j, k);
275  const double divvi =
276  fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
277  const double divve =
278  fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
279 
280  for (unsigned int d = 0; d < dim; ++d)
281  {
282  uni += fe1.shape_value_component(j, k, d) * n[d];
283  une += fe2.shape_value_component(j, k, d) * n[d];
284  vni += fe1.shape_value_component(i, k, d) * n[d];
285  vne += fe2.shape_value_component(i, k, d) * n[d];
286  }
287  M11(i, j) +=
288  dx * (-.5 * fi * divvi * uni - .5 * fi * divui * vni +
289  f * penalty * uni * vni);
290  M12(i, j) +=
291  dx * (.5 * fi * divvi * une - .5 * fe * divue * vni -
292  f * penalty * vni * une);
293  M21(i, j) +=
294  dx * (-.5 * fe * divve * uni + .5 * fi * divui * vne -
295  f * penalty * uni * vne);
296  M22(i, j) +=
297  dx * (.5 * fe * divve * une + .5 * fe * divue * vne +
298  f * penalty * une * vne);
299  }
300  }
301  }
302 
317  template <int dim>
318  void
320  Vector<double> & result2,
321  const FEValuesBase<dim> & fe1,
322  const FEValuesBase<dim> & fe2,
323  const ArrayView<const std::vector<double>> & input1,
324  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
325  const ArrayView<const std::vector<double>> & input2,
326  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
327  double pen,
328  double int_factor = 1.,
329  double ext_factor = -1.)
330  {
331  const unsigned int n1 = fe1.dofs_per_cell;
332 
333  AssertDimension(fe1.get_fe().n_components(), dim);
338 
339  const double fi = int_factor;
340  const double fe = (ext_factor < 0) ? int_factor : ext_factor;
341  const double penalty = .5 * pen * (fi + fe);
342 
343 
344  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
345  {
346  const double dx = fe1.JxW(k);
347  const Tensor<1, dim> n = fe1.normal_vector(k);
348  double uni = 0.;
349  double une = 0.;
350  double divui = 0.;
351  double divue = 0.;
352  for (unsigned int d = 0; d < dim; ++d)
353  {
354  uni += input1[d][k] * n[d];
355  une += input2[d][k] * n[d];
356  divui += Dinput1[d][k][d];
357  divue += Dinput2[d][k][d];
358  }
359 
360  for (unsigned int i = 0; i < n1; ++i)
361  {
362  double vni = 0.;
363  double vne = 0.;
364  const double divvi =
365  fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
366  const double divve =
367  fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
368  for (unsigned int d = 0; d < dim; ++d)
369  {
370  vni += fe1.shape_value_component(i, k, d) * n[d];
371  vne += fe2.shape_value_component(i, k, d) * n[d];
372  }
373 
374  result1(i) += dx * (-.5 * fi * divvi * uni -
375  .5 * fi * divui * vni + penalty * uni * vni);
376  result1(i) += dx * (.5 * fi * divvi * une -
377  .5 * fe * divue * vni - penalty * vni * une);
378  result2(i) += dx * (-.5 * fe * divve * uni +
379  .5 * fi * divui * vne - penalty * uni * vne);
380  result2(i) += dx * (.5 * fe * divve * une +
381  .5 * fe * divue * vne + penalty * une * vne);
382  }
383  }
384  }
385  } // namespace GradDiv
386 } // namespace LocalIntegrators
387 
388 DEAL_II_NAMESPACE_CLOSE
389 
390 
391 #endif
size_type m() const
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: grad_div.h:130
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const unsigned int dofs_per_cell
Definition: fe_values.h:2106
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: grad_div.h:57
void ip_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, double penalty, double factor1=1., double factor2=-1.)
Definition: grad_div.h:236
const FiniteElement< dim, spacedim > & get_fe() const
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1595
LinearAlgebra::distributed::Vector< Number > Vector
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void nitsche_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const ArrayView< const std::vector< double >> &input, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput, const ArrayView< const std::vector< double >> &data, double penalty, double factor=1.)
Definition: grad_div.h:185
Library of integrals over cells and faces.
Definition: advection.h:34
size_type n() const
void cell_residual(Vector< number > &result, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< Tensor< 1, dim >>> &input, const double factor=1.)
Definition: grad_div.h:94
const unsigned int n_quadrature_points
Definition: fe_values.h:2099
double JxW(const unsigned int quadrature_point) const
void ip_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const ArrayView< const std::vector< double >> &input1, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput1, const ArrayView< const std::vector< double >> &input2, const ArrayView< const std::vector< Tensor< 1, dim >>> &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
Definition: grad_div.h:319
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const