Reference documentation for deal.II version 9.4.0
\(\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
26#include <deal.II/fe/mapping.h>
27
29
31
33
34namespace 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
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
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;
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 unsigned int dofs_per_cell
Definition: fe_values.h:2450
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
double JxW(const unsigned int quadrature_point) const
const FiniteElement< dim, spacedim > & get_fe() const
unsigned int n_components() const
size_type n() const
size_type m() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition: exceptions.h:1690
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
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 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 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
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 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
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)