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