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