Reference documentation for deal.II version 9.5.1
\(\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// 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
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
32DEAL_II_NAMESPACE_OPEN
33
34namespace LocalIntegrators
35{
42 namespace GradDiv
43 {
50 template <int dim>
51 void
52 cell_matrix(FullMatrix<double> & M,
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
86 cell_residual(Vector<number> & result,
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);
94 AssertVectorVectorDimension(input, dim, fetest.n_quadrature_points);
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
122 nitsche_matrix(FullMatrix<double> & M,
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
174 nitsche_residual(Vector<double> & result,
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)
184 AssertVectorVectorDimension(input, dim, fe.n_quadrature_points);
185 AssertVectorVectorDimension(Dinput, dim, fe.n_quadrature_points);
186 AssertVectorVectorDimension(data, dim, fe.n_quadrature_points);
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
222 ip_matrix(FullMatrix<double> & M11,
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
302 ip_residual(Vector<double> & result1,
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);
317 AssertVectorVectorDimension(input1, dim, fe1.n_quadrature_points);
318 AssertVectorVectorDimension(Dinput1, dim, fe1.n_quadrature_points);
319 AssertVectorVectorDimension(input2, dim, fe2.n_quadrature_points);
320 AssertVectorVectorDimension(Dinput2, dim, fe2.n_quadrature_points);
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
371DEAL_II_NAMESPACE_CLOSE
372
373
374#endif
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
fe_values.h
full_matrix.h
AssertVectorVectorDimension
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
Definition exceptions.h:1810
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
dof_info.h
LocalIntegrators::GradDiv::nitsche_residual
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
LocalIntegrators::GradDiv::nitsche_matrix
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition grad_div.h:122
LocalIntegrators::GradDiv::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition grad_div.h:52
LocalIntegrators::GradDiv::ip_residual
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
LocalIntegrators::GradDiv::ip_matrix
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
LocalIntegrators::GradDiv::cell_residual
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
LocalIntegrators
Library of integrals over cells and faces.
Definition advection.h:35
Generated by doxygen 1.9.8