deal.II version GIT relicensing-2291-g668cd86249 2024-12-24 11:30:00+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 - 2024 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>
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>
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 DEAL_II_DEPRECATED 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>
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
218 template <int dim>
224 const FEValuesBase<dim> &fe1,
225 const FEValuesBase<dim> &fe2,
226 double penalty,
227 double factor1 = 1.,
228 double factor2 = -1.)
229 {
230 const unsigned int n_dofs = fe1.dofs_per_cell;
231 AssertDimension(M11.n(), n_dofs);
232 AssertDimension(M11.m(), n_dofs);
233 AssertDimension(M12.n(), n_dofs);
234 AssertDimension(M12.m(), n_dofs);
235 AssertDimension(M21.n(), n_dofs);
236 AssertDimension(M21.m(), n_dofs);
237 AssertDimension(M22.n(), n_dofs);
238 AssertDimension(M22.m(), n_dofs);
239
240 const double fi = factor1;
241 const double fe = (factor2 < 0) ? factor1 : factor2;
242 const double f = .5 * (fi + fe);
243
244 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
245 {
246 const double dx = fe1.JxW(k);
247 const Tensor<1, dim> n = fe1.normal_vector(k);
248 for (unsigned int i = 0; i < n_dofs; ++i)
249 for (unsigned int j = 0; j < n_dofs; ++j)
250 {
251 double uni = 0.;
252 double une = 0.;
253 double vni = 0.;
254 double vne = 0.;
255 const double divui =
256 fe1[FEValuesExtractors::Vector(0)].divergence(j, k);
257 const double divue =
258 fe2[FEValuesExtractors::Vector(0)].divergence(j, k);
259 const double divvi =
260 fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
261 const double divve =
262 fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
263
264 for (unsigned int d = 0; d < dim; ++d)
265 {
266 uni += fe1.shape_value_component(j, k, d) * n[d];
267 une += fe2.shape_value_component(j, k, d) * n[d];
268 vni += fe1.shape_value_component(i, k, d) * n[d];
269 vne += fe2.shape_value_component(i, k, d) * n[d];
270 }
271 M11(i, j) +=
272 dx * (-.5 * fi * divvi * uni - .5 * fi * divui * vni +
273 f * penalty * uni * vni);
274 M12(i, j) +=
275 dx * (.5 * fi * divvi * une - .5 * fe * divue * vni -
276 f * penalty * vni * une);
277 M21(i, j) +=
278 dx * (-.5 * fe * divve * uni + .5 * fi * divui * vne -
279 f * penalty * uni * vne);
280 M22(i, j) +=
281 dx * (.5 * fe * divve * une + .5 * fe * divue * vne +
282 f * penalty * une * vne);
283 }
284 }
285 }
286
298 template <int dim>
301 Vector<double> &result2,
302 const FEValuesBase<dim> &fe1,
303 const FEValuesBase<dim> &fe2,
304 const ArrayView<const std::vector<double>> &input1,
305 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput1,
306 const ArrayView<const std::vector<double>> &input2,
307 const ArrayView<const std::vector<Tensor<1, dim>>> &Dinput2,
308 double pen,
309 double int_factor = 1.,
310 double ext_factor = -1.)
311 {
312 const unsigned int n1 = fe1.dofs_per_cell;
313
314 AssertDimension(fe1.get_fe().n_components(), dim);
319
320 const double fi = int_factor;
321 const double fe = (ext_factor < 0) ? int_factor : ext_factor;
322 const double penalty = .5 * pen * (fi + fe);
323
324
325 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
326 {
327 const double dx = fe1.JxW(k);
328 const Tensor<1, dim> n = fe1.normal_vector(k);
329 double uni = 0.;
330 double une = 0.;
331 double divui = 0.;
332 double divue = 0.;
333 for (unsigned int d = 0; d < dim; ++d)
334 {
335 uni += input1[d][k] * n[d];
336 une += input2[d][k] * n[d];
337 divui += Dinput1[d][k][d];
338 divue += Dinput2[d][k][d];
339 }
340
341 for (unsigned int i = 0; i < n1; ++i)
342 {
343 double vni = 0.;
344 double vne = 0.;
345 const double divvi =
346 fe1[FEValuesExtractors::Vector(0)].divergence(i, k);
347 const double divve =
348 fe2[FEValuesExtractors::Vector(0)].divergence(i, k);
349 for (unsigned int d = 0; d < dim; ++d)
350 {
351 vni += fe1.shape_value_component(i, k, d) * n[d];
352 vne += fe2.shape_value_component(i, k, d) * n[d];
353 }
354
355 result1(i) += dx * (-.5 * fi * divvi * uni -
356 .5 * fi * divui * vni + penalty * uni * vni);
357 result1(i) += dx * (.5 * fi * divvi * une -
358 .5 * fe * divue * vni - penalty * vni * une);
359 result2(i) += dx * (-.5 * fe * divve * uni +
360 .5 * fi * divui * vne - penalty * uni * vne);
361 result2(i) += dx * (.5 * fe * divve * une +
362 .5 * fe * divue * vne + penalty * une * vne);
363 }
364 }
365 }
366 } // namespace GradDiv
367} // namespace LocalIntegrators
368
370
371
372#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_DEPRECATED
Definition config.h:205
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define AssertVectorVectorDimension(VEC, DIM1, DIM2)
#define AssertDimension(dim1, dim2)
std::vector< index_type > data
Definition mpi.cc:735
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:300
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:220
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