Reference documentation for deal.II version Git 6d63218887 2020-10-30 17:17:53 -0400
\(\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\}}\)
maxwell.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_maxwell_h
17 #define dealii_integrators_maxwell_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
29 
31 
33 
34 namespace LocalIntegrators
35 {
65  namespace Maxwell
66  {
90  template <int dim>
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>
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
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
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
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
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
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 
480 
481 #endif
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
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
const unsigned int dofs_per_cell
Definition: fe_values.h:2097
const FiniteElement< dim, spacedim > & get_fe() const
void tangential_trace_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: maxwell.h:328
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
Library of integrals over cells and faces.
Definition: advection.h:34
void curl_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: maxwell.h:163
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1423
void curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: maxwell.h:216
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int n_quadrature_points
Definition: fe_values.h:2090
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
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
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
double JxW(const unsigned int quadrature_point) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
static ::ExceptionBase & ExcNotImplemented()
Tensor< 1, dim > curl_curl(const Tensor< 2, dim > &h0, const Tensor< 2, dim > &h1, const Tensor< 2, dim > &h2)
Definition: maxwell.h:92
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const