Reference documentation for deal.II version Git 2375602 2018-06-23 12:24:17 -0500
maxwell.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2017 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 
32 DEAL_II_NAMESPACE_OPEN
33 
34 namespace LocalIntegrators
35 {
67  namespace Maxwell
68  {
95  template <int dim>
98  const Tensor<2, dim> &h1,
99  const Tensor<2, dim> &h2)
100  {
101  Tensor<1, dim> result;
102  switch (dim)
103  {
104  case 2:
105  result[0] = h1[0][1] - h0[1][1];
106  result[1] = h0[0][1] - h1[0][0];
107  break;
108  case 3:
109  result[0] = h1[0][1] + h2[0][2] - h0[1][1] - h0[2][2];
110  result[1] = h2[1][2] + h0[1][0] - h1[2][2] - h1[0][0];
111  result[2] = h0[2][0] + h1[2][1] - h2[0][0] - h2[1][1];
112  break;
113  default:
114  Assert(false, ExcNotImplemented());
115  }
116  return result;
117  }
118 
132  template <int dim>
135  const Tensor<1, dim> &g1,
136  const Tensor<1, dim> &g2,
137  const Tensor<1, dim> &normal)
138  {
139  Tensor<1, dim> result;
140 
141  switch (dim)
142  {
143  case 2:
144  result[0] = normal[1] * (g1[0] - g0[1]);
145  result[1] = -normal[0] * (g1[0] - g0[1]);
146  break;
147  case 3:
148  result[0] =
149  normal[2] * (g2[1] - g0[2]) + normal[1] * (g1[0] - g0[1]);
150  result[1] =
151  normal[0] * (g0[2] - g1[0]) + normal[2] * (g2[1] - g1[2]);
152  result[2] =
153  normal[1] * (g1[0] - g2[1]) + normal[0] * (g0[2] - g2[0]);
154  break;
155  default:
156  Assert(false, ExcNotImplemented());
157  }
158  return result;
159  }
160 
172  template <int dim>
173  void
175  const FEValuesBase<dim> &fe,
176  const double factor = 1.)
177  {
178  const unsigned int n_dofs = fe.dofs_per_cell;
179 
180  AssertDimension(fe.get_fe().n_components(), dim);
181  AssertDimension(M.m(), n_dofs);
182  AssertDimension(M.n(), n_dofs);
183 
184  // Depending on the dimension,
185  // the cross product is either
186  // a scalar (2d) or a vector
187  // (3d). Accordingly, in the
188  // latter case we have to sum
189  // up three bilinear forms, but
190  // in 2d, we don't. Thus, we
191  // need to adapt the loop over
192  // all dimensions
193  const unsigned int d_max = (dim == 2) ? 1 : dim;
194 
195  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
196  {
197  const double dx = factor * fe.JxW(k);
198  for (unsigned int i = 0; i < n_dofs; ++i)
199  for (unsigned int j = 0; j < n_dofs; ++j)
200  for (unsigned int d = 0; d < d_max; ++d)
201  {
202  const unsigned int d1 = (d + 1) % dim;
203  const unsigned int d2 = (d + 2) % dim;
204 
205  const double cv = fe.shape_grad_component(i, k, d2)[d1] -
206  fe.shape_grad_component(i, k, d1)[d2];
207  const double cu = fe.shape_grad_component(j, k, d2)[d1] -
208  fe.shape_grad_component(j, k, d1)[d2];
209 
210  M(i, j) += dx * cu * cv;
211  }
212  }
213  }
214 
228  template <int dim>
229  void
231  const FEValuesBase<dim> &fe,
232  const FEValuesBase<dim> &fetest,
233  double factor = 1.)
234  {
235  const unsigned int n_dofs = fe.dofs_per_cell;
236  const unsigned int t_dofs = fetest.dofs_per_cell;
237  AssertDimension(fe.get_fe().n_components(), dim);
238  // There should be the right number of components (3 in 3D, otherwise 1)
239  // for the curl.
240  AssertDimension(fetest.get_fe().n_components(), (dim == 3) ? dim : 1);
241  AssertDimension(M.m(), t_dofs);
242  AssertDimension(M.n(), n_dofs);
243 
244  const unsigned int d_max = (dim == 2) ? 1 : dim;
245 
246  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
247  {
248  const double dx = fe.JxW(k) * factor;
249  for (unsigned int i = 0; i < t_dofs; ++i)
250  for (unsigned int j = 0; j < n_dofs; ++j)
251  for (unsigned int d = 0; d < d_max; ++d)
252  {
253  const unsigned int d1 = (d + 1) % dim;
254  const unsigned int d2 = (d + 2) % dim;
255 
256  const double vv = fetest.shape_value_component(i, k, d);
257  const double cu = fe.shape_grad_component(j, k, d2)[d1] -
258  fe.shape_grad_component(j, k, d1)[d2];
259  M(i, j) += dx * cu * vv;
260  }
261  }
262  }
263 
280  template <int dim>
281  void
283  const FEValuesBase<dim> &fe,
284  const unsigned int face_no,
285  double penalty,
286  double factor = 1.)
287  {
288  const unsigned int n_dofs = fe.dofs_per_cell;
289 
290  AssertDimension(fe.get_fe().n_components(), dim);
291  AssertDimension(M.m(), n_dofs);
292  AssertDimension(M.n(), n_dofs);
293 
294  // Depending on the
295  // dimension, the cross
296  // product is either a scalar
297  // (2d) or a vector
298  // (3d). Accordingly, in the
299  // latter case we have to sum
300  // up three bilinear forms,
301  // but in 2d, we don't. Thus,
302  // we need to adapt the loop
303  // over all dimensions
304  const unsigned int d_max = (dim == 2) ? 1 : dim;
305 
306  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
307  {
308  const double dx = factor * fe.JxW(k);
309  const Tensor<1, dim> n = fe.normal_vector(k);
310  for (unsigned int i = 0; i < n_dofs; ++i)
311  for (unsigned int j = 0; j < n_dofs; ++j)
312  if (fe.get_fe().has_support_on_face(i, face_no) &&
313  fe.get_fe().has_support_on_face(j, face_no))
314  {
315  for (unsigned int d = 0; d < d_max; ++d)
316  {
317  const unsigned int d1 = (d + 1) % dim;
318  const unsigned int d2 = (d + 2) % dim;
319 
320  const double cv = fe.shape_grad_component(i, k, d2)[d1] -
321  fe.shape_grad_component(i, k, d1)[d2];
322  const double cu = fe.shape_grad_component(j, k, d2)[d1] -
323  fe.shape_grad_component(j, k, d1)[d2];
324  const double v =
325  fe.shape_value_component(i, k, d1) * n[d2] -
326  fe.shape_value_component(i, k, d2) * n[d1];
327  const double u =
328  fe.shape_value_component(j, k, d1) * n[d2] -
329  fe.shape_value_component(j, k, d2) * n[d1];
330 
331  M(i, j) += dx * (2. * penalty * u * v - cv * u - cu * v);
332  }
333  }
334  }
335  }
346  template <int dim>
347  void
349  const FEValuesBase<dim> &fe,
350  double factor = 1.)
351  {
352  const unsigned int n_dofs = fe.dofs_per_cell;
353 
354  AssertDimension(fe.get_fe().n_components(), dim);
355  AssertDimension(M.m(), n_dofs);
356  AssertDimension(M.n(), n_dofs);
357 
358  // Depending on the
359  // dimension, the cross
360  // product is either a scalar
361  // (2d) or a vector
362  // (3d). Accordingly, in the
363  // latter case we have to sum
364  // up three bilinear forms,
365  // but in 2d, we don't. Thus,
366  // we need to adapt the loop
367  // over all dimensions
368  const unsigned int d_max = (dim == 2) ? 1 : dim;
369 
370  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
371  {
372  const double dx = factor * fe.JxW(k);
373  const Tensor<1, dim> n = fe.normal_vector(k);
374  for (unsigned int i = 0; i < n_dofs; ++i)
375  for (unsigned int j = 0; j < n_dofs; ++j)
376  for (unsigned int d = 0; d < d_max; ++d)
377  {
378  const unsigned int d1 = (d + 1) % dim;
379  const unsigned int d2 = (d + 2) % dim;
380 
381  const double v = fe.shape_value_component(i, k, d1) * n(d2) -
382  fe.shape_value_component(i, k, d2) * n(d1);
383  const double u = fe.shape_value_component(j, k, d1) * n(d2) -
384  fe.shape_value_component(j, k, d2) * n(d1);
385 
386  M(i, j) += dx * u * v;
387  }
388  }
389  }
390 
406  template <int dim>
407  inline void
409  FullMatrix<double> & M12,
410  FullMatrix<double> & M21,
411  FullMatrix<double> & M22,
412  const FEValuesBase<dim> &fe1,
413  const FEValuesBase<dim> &fe2,
414  const double pen,
415  const double factor1 = 1.,
416  const double factor2 = -1.)
417  {
418  const unsigned int n_dofs = fe1.dofs_per_cell;
419 
420  AssertDimension(fe1.get_fe().n_components(), dim);
421  AssertDimension(fe2.get_fe().n_components(), dim);
422  AssertDimension(M11.m(), n_dofs);
423  AssertDimension(M11.n(), n_dofs);
424  AssertDimension(M12.m(), n_dofs);
425  AssertDimension(M12.n(), n_dofs);
426  AssertDimension(M21.m(), n_dofs);
427  AssertDimension(M21.n(), n_dofs);
428  AssertDimension(M22.m(), n_dofs);
429  AssertDimension(M22.n(), n_dofs);
430 
431  const double nu1 = factor1;
432  const double nu2 = (factor2 < 0) ? factor1 : factor2;
433  const double penalty = .5 * pen * (nu1 + nu2);
434 
435  // Depending on the
436  // dimension, the cross
437  // product is either a scalar
438  // (2d) or a vector
439  // (3d). Accordingly, in the
440  // latter case we have to sum
441  // up three bilinear forms,
442  // but in 2d, we don't. Thus,
443  // we need to adapt the loop
444  // over all dimensions
445  const unsigned int d_max = (dim == 2) ? 1 : dim;
446 
447  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
448  {
449  const double dx = fe1.JxW(k);
450  const Tensor<1, dim> n = fe1.normal_vector(k);
451  for (unsigned int i = 0; i < n_dofs; ++i)
452  for (unsigned int j = 0; j < n_dofs; ++j)
453  for (unsigned int d = 0; d < d_max; ++d)
454  {
455  const unsigned int d1 = (d + 1) % dim;
456  const unsigned int d2 = (d + 2) % dim;
457  // curl u, curl v
458  const double cv1 =
459  nu1 * fe1.shape_grad_component(i, k, d2)[d1] -
460  fe1.shape_grad_component(i, k, d1)[d2];
461  const double cv2 =
462  nu2 * fe2.shape_grad_component(i, k, d2)[d1] -
463  fe2.shape_grad_component(i, k, d1)[d2];
464  const double cu1 =
465  nu1 * fe1.shape_grad_component(j, k, d2)[d1] -
466  fe1.shape_grad_component(j, k, d1)[d2];
467  const double cu2 =
468  nu2 * fe2.shape_grad_component(j, k, d2)[d1] -
469  fe2.shape_grad_component(j, k, d1)[d2];
470 
471  // u x n, v x n
472  const double u1 =
473  fe1.shape_value_component(j, k, d1) * n(d2) -
474  fe1.shape_value_component(j, k, d2) * n(d1);
475  const double u2 =
476  -fe2.shape_value_component(j, k, d1) * n(d2) +
477  fe2.shape_value_component(j, k, d2) * n(d1);
478  const double v1 =
479  fe1.shape_value_component(i, k, d1) * n(d2) -
480  fe1.shape_value_component(i, k, d2) * n(d1);
481  const double v2 =
482  -fe2.shape_value_component(i, k, d1) * n(d2) +
483  fe2.shape_value_component(i, k, d2) * n(d1);
484 
485  M11(i, j) +=
486  .5 * dx * (2. * penalty * u1 * v1 - cv1 * u1 - cu1 * v1);
487  M12(i, j) +=
488  .5 * dx * (2. * penalty * v1 * u2 - cv1 * u2 - cu2 * v1);
489  M21(i, j) +=
490  .5 * dx * (2. * penalty * u1 * v2 - cv2 * u1 - cu1 * v2);
491  M22(i, j) +=
492  .5 * dx * (2. * penalty * u2 * v2 - cv2 * u2 - cu2 * v2);
493  }
494  }
495  }
496 
497 
498  } // namespace Maxwell
499 } // namespace LocalIntegrators
500 
501 
502 DEAL_II_NAMESPACE_CLOSE
503 
504 #endif
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
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:134
const unsigned int dofs_per_cell
Definition: fe_values.h:2042
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
void tangential_trace_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: maxwell.h:348
Library of integrals over cells and faces.
Definition: advection.h:34
size_type n() const
void curl_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: maxwell.h:174
#define Assert(cond, exc)
Definition: exceptions.h:1227
void curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: maxwell.h:230
const FiniteElement< dim, spacedim > & get_fe() const
size_type m() const
const unsigned int n_quadrature_points
Definition: fe_values.h:2035
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:408
void nitsche_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const unsigned int face_no, double penalty, double factor=1.)
Definition: maxwell.h:282
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:97
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
double JxW(const unsigned int quadrature_point) const