Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20: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
maxwell.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 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_maxwell_h
16#define dealii_integrators_maxwell_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{
64 namespace Maxwell
65 {
91 template <int dim>
92 DEAL_II_DEPRECATED_EARLY Tensor<1, dim>
94 const Tensor<2, dim> &h1,
95 const Tensor<2, dim> &h2)
96 {
97 Tensor<1, dim> result;
98 switch (dim)
99 {
100 case 2:
101 result[0] = h1[0][1] - h0[1][1];
102 result[1] = h0[0][1] - h1[0][0];
103 break;
104 case 3:
105 result[0] = h1[0][1] + h2[0][2] - h0[1][1] - h0[2][2];
106 result[1] = h2[1][2] + h0[1][0] - h1[2][2] - h1[0][0];
107 result[2] = h0[2][0] + h1[2][1] - h2[0][0] - h2[1][1];
108 break;
109 default:
111 }
112 return result;
113 }
114
125 template <int dim>
126 DEAL_II_DEPRECATED_EARLY Tensor<1, dim>
128 const Tensor<1, dim> &g1,
129 const Tensor<1, dim> &g2,
130 const Tensor<1, dim> &normal)
131 {
132 Tensor<1, dim> result;
133
134 switch (dim)
135 {
136 case 2:
137 result[0] = normal[1] * (g1[0] - g0[1]);
138 result[1] = -normal[0] * (g1[0] - g0[1]);
139 break;
140 case 3:
141 result[0] =
142 normal[2] * (g2[1] - g0[2]) + normal[1] * (g1[0] - g0[1]);
143 result[1] =
144 normal[0] * (g0[2] - g1[0]) + normal[2] * (g2[1] - g1[2]);
145 result[2] =
146 normal[1] * (g1[0] - g2[1]) + normal[0] * (g0[2] - g2[0]);
147 break;
148 default:
150 }
151 return result;
152 }
153
162 template <int dim>
163 DEAL_II_DEPRECATED_EARLY void
165 const FEValuesBase<dim> &fe,
166 const double factor = 1.)
167 {
168 const unsigned int n_dofs = fe.dofs_per_cell;
169
171 AssertDimension(M.m(), n_dofs);
172 AssertDimension(M.n(), n_dofs);
173
174 // Depending on the dimension,
175 // the cross product is either
176 // a scalar (2d) or a vector
177 // (3d). Accordingly, in the
178 // latter case we have to sum
179 // up three bilinear forms, but
180 // in 2d, we don't. Thus, we
181 // need to adapt the loop over
182 // all dimensions
183 const unsigned int d_max = (dim == 2) ? 1 : dim;
184
185 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
186 {
187 const double dx = factor * fe.JxW(k);
188 for (unsigned int i = 0; i < n_dofs; ++i)
189 for (unsigned int j = 0; j < n_dofs; ++j)
190 for (unsigned int d = 0; d < d_max; ++d)
191 {
192 const unsigned int d1 = (d + 1) % dim;
193 const unsigned int d2 = (d + 2) % dim;
194
195 const double cv = fe.shape_grad_component(i, k, d2)[d1] -
196 fe.shape_grad_component(i, k, d1)[d2];
197 const double cu = fe.shape_grad_component(j, k, d2)[d1] -
198 fe.shape_grad_component(j, k, d1)[d2];
199
200 M(i, j) += dx * cu * cv;
201 }
202 }
203 }
204
215 template <int dim>
216 DEAL_II_DEPRECATED_EARLY void
218 const FEValuesBase<dim> &fe,
219 const FEValuesBase<dim> &fetest,
220 double factor = 1.)
221 {
222 const unsigned int n_dofs = fe.dofs_per_cell;
223 const unsigned int t_dofs = fetest.dofs_per_cell;
225 // There should be the right number of components (3 in 3d, otherwise 1)
226 // for the curl.
227 AssertDimension(fetest.get_fe().n_components(), (dim == 3) ? dim : 1);
228 AssertDimension(M.m(), t_dofs);
229 AssertDimension(M.n(), n_dofs);
230
231 const unsigned int d_max = (dim == 2) ? 1 : dim;
232
233 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
234 {
235 const double dx = fe.JxW(k) * factor;
236 for (unsigned int i = 0; i < t_dofs; ++i)
237 for (unsigned int j = 0; j < n_dofs; ++j)
238 for (unsigned int d = 0; d < d_max; ++d)
239 {
240 const unsigned int d1 = (d + 1) % dim;
241 const unsigned int d2 = (d + 2) % dim;
242
243 const double vv = fetest.shape_value_component(i, k, d);
244 const double cu = fe.shape_grad_component(j, k, d2)[d1] -
245 fe.shape_grad_component(j, k, d1)[d2];
246 M(i, j) += dx * cu * vv;
247 }
248 }
249 }
250
264 template <int dim>
265 void
267 const FEValuesBase<dim> &fe,
268 const unsigned int face_no,
269 double penalty,
270 double factor = 1.)
271 {
272 const unsigned int n_dofs = fe.dofs_per_cell;
273
275 AssertDimension(M.m(), n_dofs);
276 AssertDimension(M.n(), n_dofs);
277
278 // Depending on the
279 // dimension, the cross
280 // product is either a scalar
281 // (2d) or a vector
282 // (3d). Accordingly, in the
283 // latter case we have to sum
284 // up three bilinear forms,
285 // but in 2d, we don't. Thus,
286 // we need to adapt the loop
287 // over all dimensions
288 const unsigned int d_max = (dim == 2) ? 1 : dim;
289
290 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
291 {
292 const double dx = factor * fe.JxW(k);
293 const Tensor<1, dim> n = fe.normal_vector(k);
294 for (unsigned int i = 0; i < n_dofs; ++i)
295 for (unsigned int j = 0; j < n_dofs; ++j)
296 if (fe.get_fe().has_support_on_face(i, face_no) &&
297 fe.get_fe().has_support_on_face(j, face_no))
298 {
299 for (unsigned int d = 0; d < d_max; ++d)
300 {
301 const unsigned int d1 = (d + 1) % dim;
302 const unsigned int d2 = (d + 2) % dim;
303
304 const double cv = fe.shape_grad_component(i, k, d2)[d1] -
305 fe.shape_grad_component(i, k, d1)[d2];
306 const double cu = fe.shape_grad_component(j, k, d2)[d1] -
307 fe.shape_grad_component(j, k, d1)[d2];
308 const double v =
309 fe.shape_value_component(i, k, d1) * n[d2] -
310 fe.shape_value_component(i, k, d2) * n[d1];
311 const double u =
312 fe.shape_value_component(j, k, d1) * n[d2] -
313 fe.shape_value_component(j, k, d2) * n[d1];
314
315 M(i, j) += dx * (2. * penalty * u * v - cv * u - cu * v);
316 }
317 }
318 }
319 }
327 template <int dim>
328 DEAL_II_DEPRECATED_EARLY void
330 const FEValuesBase<dim> &fe,
331 double factor = 1.)
332 {
333 const unsigned int n_dofs = fe.dofs_per_cell;
334
336 AssertDimension(M.m(), n_dofs);
337 AssertDimension(M.n(), n_dofs);
338
339 // Depending on the
340 // dimension, the cross
341 // product is either a scalar
342 // (2d) or a vector
343 // (3d). Accordingly, in the
344 // latter case we have to sum
345 // up three bilinear forms,
346 // but in 2d, we don't. Thus,
347 // we need to adapt the loop
348 // over all dimensions
349 const unsigned int d_max = (dim == 2) ? 1 : dim;
350
351 for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
352 {
353 const double dx = factor * fe.JxW(k);
354 const Tensor<1, dim> n = fe.normal_vector(k);
355 for (unsigned int i = 0; i < n_dofs; ++i)
356 for (unsigned int j = 0; j < n_dofs; ++j)
357 for (unsigned int d = 0; d < d_max; ++d)
358 {
359 const unsigned int d1 = (d + 1) % dim;
360 const unsigned int d2 = (d + 2) % dim;
361
362 const double v = fe.shape_value_component(i, k, d1) * n(d2) -
363 fe.shape_value_component(i, k, d2) * n(d1);
364 const double u = fe.shape_value_component(j, k, d1) * n(d2) -
365 fe.shape_value_component(j, k, d2) * n(d1);
366
367 M(i, j) += dx * u * v;
368 }
369 }
370 }
371
384 template <int dim>
385 DEAL_II_DEPRECATED_EARLY inline void
390 const FEValuesBase<dim> &fe1,
391 const FEValuesBase<dim> &fe2,
392 const double pen,
393 const double factor1 = 1.,
394 const double factor2 = -1.)
395 {
396 const unsigned int n_dofs = fe1.n_dofs_per_cell();
397
398 AssertDimension(fe1.get_fe().n_components(), dim);
399 AssertDimension(fe2.get_fe().n_components(), dim);
400 AssertDimension(M11.m(), n_dofs);
401 AssertDimension(M11.n(), n_dofs);
402 AssertDimension(M12.m(), n_dofs);
403 AssertDimension(M12.n(), n_dofs);
404 AssertDimension(M21.m(), n_dofs);
405 AssertDimension(M21.n(), n_dofs);
406 AssertDimension(M22.m(), n_dofs);
407 AssertDimension(M22.n(), n_dofs);
408
409 const double nu1 = factor1;
410 const double nu2 = (factor2 < 0) ? factor1 : factor2;
411 const double penalty = .5 * pen * (nu1 + nu2);
412
413 // Depending on the
414 // dimension, the cross
415 // product is either a scalar
416 // (2d) or a vector
417 // (3d). Accordingly, in the
418 // latter case we have to sum
419 // up three bilinear forms,
420 // but in 2d, we don't. Thus,
421 // we need to adapt the loop
422 // over all dimensions
423 const unsigned int d_max = (dim == 2) ? 1 : dim;
424
425 for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
426 {
427 const double dx = fe1.JxW(k);
428 const Tensor<1, dim> n = fe1.normal_vector(k);
429 for (unsigned int i = 0; i < n_dofs; ++i)
430 for (unsigned int j = 0; j < n_dofs; ++j)
431 for (unsigned int d = 0; d < d_max; ++d)
432 {
433 const unsigned int d1 = (d + 1) % dim;
434 const unsigned int d2 = (d + 2) % dim;
435 // curl u, curl v
436 const double cv1 =
437 nu1 * fe1.shape_grad_component(i, k, d2)[d1] -
438 fe1.shape_grad_component(i, k, d1)[d2];
439 const double cv2 =
440 nu2 * fe2.shape_grad_component(i, k, d2)[d1] -
441 fe2.shape_grad_component(i, k, d1)[d2];
442 const double cu1 =
443 nu1 * fe1.shape_grad_component(j, k, d2)[d1] -
444 fe1.shape_grad_component(j, k, d1)[d2];
445 const double cu2 =
446 nu2 * fe2.shape_grad_component(j, k, d2)[d1] -
447 fe2.shape_grad_component(j, k, d1)[d2];
448
449 // u x n, v x n
450 const double u1 =
451 fe1.shape_value_component(j, k, d1) * n(d2) -
452 fe1.shape_value_component(j, k, d2) * n(d1);
453 const double u2 =
454 -fe2.shape_value_component(j, k, d1) * n(d2) +
455 fe2.shape_value_component(j, k, d2) * n(d1);
456 const double v1 =
457 fe1.shape_value_component(i, k, d1) * n(d2) -
458 fe1.shape_value_component(i, k, d2) * n(d1);
459 const double v2 =
460 -fe2.shape_value_component(i, k, d1) * n(d2) +
461 fe2.shape_value_component(i, k, d2) * n(d1);
462
463 M11(i, j) +=
464 .5 * dx * (2. * penalty * u1 * v1 - cv1 * u1 - cu1 * v1);
465 M12(i, j) +=
466 .5 * dx * (2. * penalty * v1 * u2 - cv1 * u2 - cu2 * v1);
467 M21(i, j) +=
468 .5 * dx * (2. * penalty * u1 * v2 - cv2 * u1 - cu1 * v2);
469 M22(i, j) +=
470 .5 * dx * (2. * penalty * u2 * v2 - cv2 * u2 - cu2 * v2);
471 }
472 }
473 }
474
475
476 } // namespace Maxwell
477} // namespace LocalIntegrators
478
479
481
482#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
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int q_point) const
unsigned int n_components() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) 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
const unsigned int v1
#define AssertDimension(dim1, dim2)
#define DEAL_II_NOT_IMPLEMENTED()
void nitsche_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const unsigned int face_no, double penalty, double factor=1.)
Definition maxwell.h:266
void curl_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition maxwell.h:164
void tangential_trace_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition maxwell.h:329
void curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition maxwell.h:217
Tensor< 1, dim > tangential_curl(const Tensor< 1, dim > &g0, const Tensor< 1, dim > &g1, const Tensor< 1, dim > &g2, const Tensor< 1, dim > &normal)
Definition maxwell.h:127
void ip_curl_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const double pen, const double factor1=1., const double factor2=-1.)
Definition maxwell.h:386
Tensor< 1, dim > curl_curl(const Tensor< 2, dim > &h0, const Tensor< 2, dim > &h1, const Tensor< 2, dim > &h2)
Definition maxwell.h:93
Library of integrals over cells and faces.
Definition advection.h:34