Reference documentation for deal.II version GIT d6cf33b98c 2023-09-22 19: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\}}\)
grad_div.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_grad_div_h
17 #define dealii_integrators_grad_div_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 
34 namespace LocalIntegrators
35 {
42  namespace GradDiv
43  {
50  template <int dim>
51  void
53  const FEValuesBase<dim> &fe,
54  double factor = 1.)
55  {
56  const unsigned int n_dofs = fe.dofs_per_cell;
57 
58  AssertDimension(fe.get_fe().n_components(), dim);
59  AssertDimension(M.m(), n_dofs);
60  AssertDimension(M.n(), n_dofs);
61 
62  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
63  {
64  const double dx = factor * fe.JxW(k);
65  for (unsigned int i = 0; i < n_dofs; ++i)
66  for (unsigned int j = 0; j < n_dofs; ++j)
67  {
68  const double divu =
69  fe[FEValuesExtractors::Vector(0)].divergence(j, k);
70  const double divv =
71  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
72 
73  M(i, j) += dx * divu * divv;
74  }
75  }
76  }
77 
84  template <int dim, typename number>
85  void
87  const FEValuesBase<dim> &fetest,
88  const ArrayView<const std::vector<Tensor<1, dim>>> &input,
89  const double factor = 1.)
90  {
91  const unsigned int n_dofs = fetest.dofs_per_cell;
92 
93  AssertDimension(fetest.get_fe().n_components(), dim);
95 
96  for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
97  {
98  const double dx = factor * fetest.JxW(k);
99  for (unsigned int i = 0; i < n_dofs; ++i)
100  {
101  const double divv =
102  fetest[FEValuesExtractors::Vector(0)].divergence(i, k);
103  double du = 0.;
104  for (unsigned int d = 0; d < dim; ++d)
105  du += input[d][k][d];
106 
107  result(i) += dx * du * divv;
108  }
109  }
110  }
111 
120  template <int dim>
121  inline void
123  const FEValuesBase<dim> &fe,
124  double penalty,
125  double factor = 1.)
126  {
127  const unsigned int n_dofs = fe.dofs_per_cell;
128 
129  AssertDimension(fe.get_fe().n_components(), dim);
130  AssertDimension(M.m(), n_dofs);
131  AssertDimension(M.n(), n_dofs);
132 
133  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
134  {
135  const double dx = factor * fe.JxW(k);
136  const Tensor<1, dim> n = fe.normal_vector(k);
137  for (unsigned int i = 0; i < n_dofs; ++i)
138  for (unsigned int j = 0; j < n_dofs; ++j)
139  {
140  const double divu =
141  fe[FEValuesExtractors::Vector(0)].divergence(j, k);
142  const double divv =
143  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
144  double un = 0., vn = 0.;
145  for (unsigned int d = 0; d < dim; ++d)
146  {
147  un += fe.shape_value_component(j, k, d) * n[d];
148  vn += fe.shape_value_component(i, k, d) * n[d];
149  }
150 
151  M(i, j) += dx * 2. * penalty * un * vn;
152  M(i, j) -= dx * (divu * vn + divv * un);
153  }
154  }
155  }
156 
172  template <int dim>
173  void
175  const FEValuesBase<dim> &fe,
176  const ArrayView<const std::vector<double>> &input,
177  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
178  const ArrayView<const std::vector<double>> &data,
179  double penalty,
180  double factor = 1.)
181  {
182  const unsigned int n_dofs = fe.dofs_per_cell;
183  AssertDimension(fe.get_fe().n_components(), dim)
187 
188  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
189  {
190  const double dx = factor * fe.JxW(k);
191  const Tensor<1, dim> n = fe.normal_vector(k);
192 
193  double umgn = 0.;
194  double divu = 0.;
195  for (unsigned int d = 0; d < dim; ++d)
196  {
197  umgn += (input[d][k] - data[d][k]) * n[d];
198  divu += Dinput[d][k][d];
199  }
200 
201  for (unsigned int i = 0; i < n_dofs; ++i)
202  {
203  double vn = 0.;
204  const double divv =
205  fe[FEValuesExtractors::Vector(0)].divergence(i, k);
206  for (unsigned int d = 0; d < dim; ++d)
207  vn += fe.shape_value_component(i, k, d) * n[d];
208 
209  result(i) +=
210  dx * (2. * penalty * umgn * vn - divv * umgn - divu * vn);
211  }
212  }
213  }
214 
220  template <int dim>
221  void
223  FullMatrix<double> &M12,
224  FullMatrix<double> &M21,
225  FullMatrix<double> &M22,
226  const FEValuesBase<dim> &fe1,
227  const FEValuesBase<dim> &fe2,
228  double penalty,
229  double factor1 = 1.,
230  double factor2 = -1.)
231  {
232  const unsigned int n_dofs = fe1.dofs_per_cell;
233  AssertDimension(M11.n(), n_dofs);
234  AssertDimension(M11.m(), n_dofs);
235  AssertDimension(M12.n(), n_dofs);
236  AssertDimension(M12.m(), n_dofs);
237  AssertDimension(M21.n(), n_dofs);
238  AssertDimension(M21.m(), n_dofs);
239  AssertDimension(M22.n(), n_dofs);
240  AssertDimension(M22.m(), n_dofs);
241 
242  const double fi = factor1;
243  const double fe = (factor2 < 0) ? factor1 : factor2;
244  const double f = .5 * (fi + fe);
245 
246  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
247  {
248  const double dx = fe1.JxW(k);
249  const Tensor<1, dim> n = fe1.normal_vector(k);
250  for (unsigned int i = 0; i < n_dofs; ++i)
251  for (unsigned int j = 0; j < n_dofs; ++j)
252  {
253  double uni = 0.;
254  double une = 0.;
255  double vni = 0.;
256  double vne = 0.;
257  const double divui =
258  fe1[FEValuesExtractors::Vector(0)].divergence(j, k);
259  const double divue =
260  fe2[FEValuesExtractors::Vector(0)].divergence(j, k);
261  const double divvi =
262  fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
263  const double divve =
264  fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
265 
266  for (unsigned int d = 0; d < dim; ++d)
267  {
268  uni += fe1.shape_value_component(j, k, d) * n[d];
269  une += fe2.shape_value_component(j, k, d) * n[d];
270  vni += fe1.shape_value_component(i, k, d) * n[d];
271  vne += fe2.shape_value_component(i, k, d) * n[d];
272  }
273  M11(i, j) +=
274  dx * (-.5 * fi * divvi * uni - .5 * fi * divui * vni +
275  f * penalty * uni * vni);
276  M12(i, j) +=
277  dx * (.5 * fi * divvi * une - .5 * fe * divue * vni -
278  f * penalty * vni * une);
279  M21(i, j) +=
280  dx * (-.5 * fe * divve * uni + .5 * fi * divui * vne -
281  f * penalty * uni * vne);
282  M22(i, j) +=
283  dx * (.5 * fe * divve * une + .5 * fe * divue * vne +
284  f * penalty * une * vne);
285  }
286  }
287  }
288 
300  template <int dim>
301  void
303  Vector<double> &result2,
304  const FEValuesBase<dim> &fe1,
305  const FEValuesBase<dim> &fe2,
306  const ArrayView<const std::vector<double>> &input1,
307  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
308  const ArrayView<const std::vector<double>> &input2,
309  const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
310  double pen,
311  double int_factor = 1.,
312  double ext_factor = -1.)
313  {
314  const unsigned int n1 = fe1.dofs_per_cell;
315 
316  AssertDimension(fe1.get_fe().n_components(), dim);
321 
322  const double fi = int_factor;
323  const double fe = (ext_factor < 0) ? int_factor : ext_factor;
324  const double penalty = .5 * pen * (fi + fe);
325 
326 
327  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
328  {
329  const double dx = fe1.JxW(k);
330  const Tensor<1, dim> n = fe1.normal_vector(k);
331  double uni = 0.;
332  double une = 0.;
333  double divui = 0.;
334  double divue = 0.;
335  for (unsigned int d = 0; d < dim; ++d)
336  {
337  uni += input1[d][k] * n[d];
338  une += input2[d][k] * n[d];
339  divui += Dinput1[d][k][d];
340  divue += Dinput2[d][k][d];
341  }
342 
343  for (unsigned int i = 0; i < n1; ++i)
344  {
345  double vni = 0.;
346  double vne = 0.;
347  const double divvi =
348  fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
349  const double divve =
350  fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
351  for (unsigned int d = 0; d < dim; ++d)
352  {
353  vni += fe1.shape_value_component(i, k, d) * n[d];
354  vne += fe2.shape_value_component(i, k, d) * n[d];
355  }
356 
357  result1(i) += dx * (-.5 * fi * divvi * uni -
358  .5 * fi * divui * vni + penalty * uni * vni);
359  result1(i) += dx * (.5 * fi * divvi * une -
360  .5 * fe * divue * vni - penalty * vni * une);
361  result2(i) += dx * (-.5 * fe * divve * uni +
362  .5 * fi * divui * vne - penalty * uni * vne);
363  result2(i) += dx * (.5 * fe * divve * une +
364  .5 * fe * divue * vne + penalty * une * vne);
365  }
366  }
367  }
368  } // namespace GradDiv
369 } // namespace LocalIntegrators
370 
372 
373 
374 #endif
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 FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
unsigned int n_components() const
size_type n() const
size_type m() const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1812
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
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:86
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: grad_div.h:122
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: grad_div.h:52
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:174
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:222
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:302
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)