Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07:30: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\}}\)
Loading...
Searching...
No Matches
grad_div.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_integrators_grad_div_h
16#define dealii_integrators_grad_div_h
17
18
19#include <deal.II/base/config.h>
20
23
25#include <deal.II/fe/mapping.h>
26
28
30
32
33namespace LocalIntegrators
34{
41 namespace GradDiv
42 {
49 template <int dim>
50 void
52 const FEValuesBase<dim> &fe,
53 double factor = 1.)
54 {
55 const unsigned int n_dofs = fe.dofs_per_cell;
56
58 AssertDimension(M.m(), n_dofs);
59 AssertDimension(M.n(), n_dofs);
60
61 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
62 {
63 const double dx = factor * fe.JxW(k);
64 for (unsigned int i = 0; i < n_dofs; ++i)
65 for (unsigned int j = 0; j < n_dofs; ++j)
66 {
67 const double divu =
68 fe[FEValuesExtractors::Vector(0)].divergence(j, k);
69 const double divv =
70 fe[FEValuesExtractors::Vector(0)].divergence(i, k);
71
72 M(i, j) += dx * divu * divv;
73 }
74 }
75 }
76
83 template <int dim, typename number>
84 void
86 const FEValuesBase<dim> &fetest,
87 const ArrayView<const std::vector<Tensor<1, dim>>> &input,
88 const double factor = 1.)
89 {
90 const unsigned int n_dofs = fetest.dofs_per_cell;
91
92 AssertDimension(fetest.get_fe().n_components(), dim);
94
95 for (unsigned int k = 0; k < fetest.n_quadrature_points; ++k)
96 {
97 const double dx = factor * fetest.JxW(k);
98 for (unsigned int i = 0; i < n_dofs; ++i)
99 {
100 const double divv =
101 fetest[FEValuesExtractors::Vector(0)].divergence(i, k);
102 double du = 0.;
103 for (unsigned int d = 0; d < dim; ++d)
104 du += input[d][k][d];
105
106 result(i) += dx * du * divv;
107 }
108 }
109 }
110
119 template <int dim>
120 inline void
122 const FEValuesBase<dim> &fe,
123 double penalty,
124 double factor = 1.)
125 {
126 const unsigned int n_dofs = fe.dofs_per_cell;
127
129 AssertDimension(M.m(), n_dofs);
130 AssertDimension(M.n(), n_dofs);
131
132 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
133 {
134 const double dx = factor * fe.JxW(k);
135 const Tensor<1, dim> n = fe.normal_vector(k);
136 for (unsigned int i = 0; i < n_dofs; ++i)
137 for (unsigned int j = 0; j < n_dofs; ++j)
138 {
139 const double divu =
140 fe[FEValuesExtractors::Vector(0)].divergence(j, k);
141 const double divv =
142 fe[FEValuesExtractors::Vector(0)].divergence(i, k);
143 double un = 0., vn = 0.;
144 for (unsigned int d = 0; d < dim; ++d)
145 {
146 un += fe.shape_value_component(j, k, d) * n[d];
147 vn += fe.shape_value_component(i, k, d) * n[d];
148 }
149
150 M(i, j) += dx * 2. * penalty * un * vn;
151 M(i, j) -= dx * (divu * vn + divv * un);
152 }
153 }
154 }
155
171 template <int dim>
172 void
174 const FEValuesBase<dim> &fe,
175 const ArrayView<const std::vector<double>> &input,
176 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput,
177 const ArrayView<const std::vector<double>> &data,
178 double penalty,
179 double factor = 1.)
180 {
181 const unsigned int n_dofs = fe.dofs_per_cell;
186
187 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
188 {
189 const double dx = factor * fe.JxW(k);
190 const Tensor<1, dim> n = fe.normal_vector(k);
191
192 double umgn = 0.;
193 double divu = 0.;
194 for (unsigned int d = 0; d < dim; ++d)
195 {
196 umgn += (input[d][k] - data[d][k]) * n[d];
197 divu += Dinput[d][k][d];
198 }
199
200 for (unsigned int i = 0; i < n_dofs; ++i)
201 {
202 double vn = 0.;
203 const double divv =
204 fe[FEValuesExtractors::Vector(0)].divergence(i, k);
205 for (unsigned int d = 0; d < dim; ++d)
206 vn += fe.shape_value_component(i, k, d) * n[d];
207
208 result(i) +=
209 dx * (2. * penalty * umgn * vn - divv * umgn - divu * vn);
210 }
211 }
212 }
213
219 template <int dim>
220 void
225 const FEValuesBase<dim> &fe1,
226 const FEValuesBase<dim> &fe2,
227 double penalty,
228 double factor1 = 1.,
229 double factor2 = -1.)
230 {
231 const unsigned int n_dofs = fe1.dofs_per_cell;
232 AssertDimension(M11.n(), n_dofs);
233 AssertDimension(M11.m(), n_dofs);
234 AssertDimension(M12.n(), n_dofs);
235 AssertDimension(M12.m(), n_dofs);
236 AssertDimension(M21.n(), n_dofs);
237 AssertDimension(M21.m(), n_dofs);
238 AssertDimension(M22.n(), n_dofs);
239 AssertDimension(M22.m(), n_dofs);
240
241 const double fi = factor1;
242 const double fe = (factor2 < 0) ? factor1 : factor2;
243 const double f = .5 * (fi + fe);
244
245 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
246 {
247 const double dx = fe1.JxW(k);
248 const Tensor<1, dim> n = fe1.normal_vector(k);
249 for (unsigned int i = 0; i < n_dofs; ++i)
250 for (unsigned int j = 0; j < n_dofs; ++j)
251 {
252 double uni = 0.;
253 double une = 0.;
254 double vni = 0.;
255 double vne = 0.;
256 const double divui =
257 fe1[FEValuesExtractors::Vector(0)].divergence(j, k);
258 const double divue =
259 fe2[FEValuesExtractors::Vector(0)].divergence(j, k);
260 const double divvi =
261 fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
262 const double divve =
263 fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
264
265 for (unsigned int d = 0; d < dim; ++d)
266 {
267 uni += fe1.shape_value_component(j, k, d) * n[d];
268 une += fe2.shape_value_component(j, k, d) * n[d];
269 vni += fe1.shape_value_component(i, k, d) * n[d];
270 vne += fe2.shape_value_component(i, k, d) * n[d];
271 }
272 M11(i, j) +=
273 dx * (-.5 * fi * divvi * uni - .5 * fi * divui * vni +
274 f * penalty * uni * vni);
275 M12(i, j) +=
276 dx * (.5 * fi * divvi * une - .5 * fe * divue * vni -
277 f * penalty * vni * une);
278 M21(i, j) +=
279 dx * (-.5 * fe * divve * uni + .5 * fi * divui * vne -
280 f * penalty * uni * vne);
281 M22(i, j) +=
282 dx * (.5 * fe * divve * une + .5 * fe * divue * vne +
283 f * penalty * une * vne);
284 }
285 }
286 }
287
299 template <int dim>
300 void
302 Vector<double> &result2,
303 const FEValuesBase<dim> &fe1,
304 const FEValuesBase<dim> &fe2,
305 const ArrayView<const std::vector<double>> &input1,
306 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
307 const ArrayView<const std::vector<double>> &input2,
308 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
309 double pen,
310 double int_factor = 1.,
311 double ext_factor = -1.)
312 {
313 const unsigned int n1 = fe1.dofs_per_cell;
314
315 AssertDimension(fe1.get_fe().n_components(), dim);
320
321 const double fi = int_factor;
322 const double fe = (ext_factor < 0) ? int_factor : ext_factor;
323 const double penalty = .5 * pen * (fi + fe);
324
325
326 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
327 {
328 const double dx = fe1.JxW(k);
329 const Tensor<1, dim> n = fe1.normal_vector(k);
330 double uni = 0.;
331 double une = 0.;
332 double divui = 0.;
333 double divue = 0.;
334 for (unsigned int d = 0; d < dim; ++d)
335 {
336 uni += input1[d][k] * n[d];
337 une += input2[d][k] * n[d];
338 divui += Dinput1[d][k][d];
339 divue += Dinput2[d][k][d];
340 }
341
342 for (unsigned int i = 0; i < n1; ++i)
343 {
344 double vni = 0.;
345 double vne = 0.;
346 const double divvi =
347 fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
348 const double divve =
349 fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
350 for (unsigned int d = 0; d < dim; ++d)
351 {
352 vni += fe1.shape_value_component(i, k, d) * n[d];
353 vne += fe2.shape_value_component(i, k, d) * n[d];
354 }
355
356 result1(i) += dx * (-.5 * fi * divvi * uni -
357 .5 * fi * divui * vni + penalty * uni * vni);
358 result1(i) += dx * (.5 * fi * divvi * une -
359 .5 * fe * divue * vni - penalty * vni * une);
360 result2(i) += dx * (-.5 * fe * divve * uni +
361 .5 * fi * divui * vne - penalty * uni * vne);
362 result2(i) += dx * (.5 * fe * divve * une +
363 .5 * fe * divue * vne + penalty * une * vne);
364 }
365 }
366 }
367 } // namespace GradDiv
368} // namespace LocalIntegrators
369
371
372
373#endif
const unsigned int dofs_per_cell
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
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:173
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition grad_div.h:121
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition grad_div.h:51
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:301
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:221
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:85
Library of integrals over cells and faces.
Definition advection.h:34