Reference documentation for deal.II version 8.1.0
maxwell.h
1 // ---------------------------------------------------------------------
2 // @f$Id: maxwell.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__integrators_maxwell_h
18 #define __deal2__integrators_maxwell_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
28 
30 
31 namespace LocalIntegrators
32 {
65  namespace Maxwell
66  {
93  template <int dim>
96  const Tensor<2,dim> &h0,
97  const Tensor<2,dim> &h1,
98  const Tensor<2,dim> &h2)
99  {
100  Tensor<1,dim> result;
101  switch (dim)
102  {
103  case 2:
104  result[0] = h1[0][1]-h0[1][1];
105  result[1] = h0[0][1]-h1[0][0];
106  break;
107  case 3:
108  result[0] = h1[0][1]+h2[0][2]-h0[1][1]-h0[2][2];
109  result[1] = h2[1][2]+h0[1][0]-h1[2][2]-h1[0][0];
110  result[2] = h0[2][0]+h1[2][1]-h2[0][0]-h2[1][1];
111  break;
112  default:
113  Assert(false, ExcNotImplemented());
114  }
115  return result;
116  }
117 
131  template <int dim>
134  const Tensor<1,dim> &g0,
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] = normal[2]*(g2[1]-g0[2])+normal[1]*(g1[0]-g0[1]);
149  result[1] = normal[0]*(g0[2]-g1[0])+normal[2]*(g2[1]-g1[2]);
150  result[2] = normal[1]*(g1[0]-g2[1])+normal[0]*(g0[2]-g2[0]);
151  break;
152  default:
153  Assert(false, ExcNotImplemented());
154  }
155  return result;
156  }
157 
169  template <int dim>
172  const FEValuesBase<dim> &fe,
173  const double factor = 1.)
174  {
175  const unsigned int n_dofs = fe.dofs_per_cell;
176 
177  AssertDimension(fe.get_fe().n_components(), dim);
178  AssertDimension(M.m(), n_dofs);
179  AssertDimension(M.n(), n_dofs);
180 
181  // Depending on the dimension,
182  // the cross product is either
183  // a scalar (2d) or a vector
184  // (3d). Accordingly, in the
185  // latter case we have to sum
186  // up three bilinear forms, but
187  // in 2d, we don't. Thus, we
188  // need to adapt the loop over
189  // all dimensions
190  const unsigned int d_max = (dim==2) ? 1 : dim;
191 
192  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
193  {
194  const double dx = factor * fe.JxW(k);
195  for (unsigned int i=0; i<n_dofs; ++i)
196  for (unsigned int j=0; j<n_dofs; ++j)
197  for (unsigned int d=0; d<d_max; ++d)
198  {
199  const unsigned int d1 = (d+1)%dim;
200  const unsigned int d2 = (d+2)%dim;
201 
202  const double cv = fe.shape_grad_component(i,k,d1)[d2] - fe.shape_grad_component(i,k,d2)[d1];
203  const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
204 
205  M(i,j) += dx * cu * cv;
206  }
207  }
208  }
209 
223  template <int dim>
224  void curl_matrix (
226  const FEValuesBase<dim> &fe,
227  const FEValuesBase<dim> &fetest,
228  double factor = 1.)
229  {
230  unsigned int t_comp = (dim==3) ? dim : 1;
231  const unsigned int n_dofs = fe.dofs_per_cell;
232  const unsigned int t_dofs = fetest.dofs_per_cell;
233  AssertDimension(fe.get_fe().n_components(), dim);
234  AssertDimension(fetest.get_fe().n_components(), t_comp);
235  AssertDimension(M.m(), t_dofs);
236  AssertDimension(M.n(), n_dofs);
237 
238  const unsigned int d_max = (dim==2) ? 1 : dim;
239 
240  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
241  {
242  const double dx = fe.JxW(k) * factor;
243  for (unsigned int i=0; i<t_dofs; ++i)
244  for (unsigned int j=0; j<n_dofs; ++j)
245  for (unsigned int d=0; d<d_max; ++d)
246  {
247  const unsigned int d1 = (d+1)%dim;
248  const unsigned int d2 = (d+2)%dim;
249 
250  const double vv = fetest.shape_value_component(i,k,d);
251  const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
252  M(i,j) += dx * cu * vv;
253  }
254  }
255  }
256 
275  template <int dim>
278  const FEValuesBase<dim> &fe,
279  double penalty,
280  double factor = 1.)
281  {
282  const unsigned int n_dofs = fe.dofs_per_cell;
283 
284  AssertDimension(fe.get_fe().n_components(), dim);
285  AssertDimension(M.m(), n_dofs);
286  AssertDimension(M.n(), n_dofs);
287 
288  // Depending on the
289  // dimension, the cross
290  // product is either a scalar
291  // (2d) or a vector
292  // (3d). Accordingly, in the
293  // latter case we have to sum
294  // up three bilinear forms,
295  // but in 2d, we don't. Thus,
296  // we need to adapt the loop
297  // over all dimensions
298  const unsigned int d_max = (dim==2) ? 1 : dim;
299 
300  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
301  {
302  const double dx = factor * fe.JxW(k);
303  const Point<dim> &n = fe.normal_vector(k);
304  for (unsigned int i=0; i<n_dofs; ++i)
305  for (unsigned int j=0; j<n_dofs; ++j)
306  for (unsigned int d=0; d<d_max; ++d)
307  {
308  const unsigned int d1 = (d+1)%dim;
309  const unsigned int d2 = (d+2)%dim;
310 
311  const double cv = fe.shape_grad_component(i,k,d1)[d2] - fe.shape_grad_component(i,k,d2)[d1];
312  const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
313  const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
314  const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
315 
316  M(i,j) += dx*(2.*penalty*u*v - cv*u - cu*v);
317  }
318  }
319  }
331  template <int dim>
334  const FEValuesBase<dim> &fe,
335  double factor = 1.)
336  {
337  const unsigned int n_dofs = fe.dofs_per_cell;
338 
339  AssertDimension(fe.get_fe().n_components(), dim);
340  AssertDimension(M.m(), n_dofs);
341  AssertDimension(M.n(), n_dofs);
342 
343  // Depending on the
344  // dimension, the cross
345  // product is either a scalar
346  // (2d) or a vector
347  // (3d). Accordingly, in the
348  // latter case we have to sum
349  // up three bilinear forms,
350  // but in 2d, we don't. Thus,
351  // we need to adapt the loop
352  // over all dimensions
353  const unsigned int d_max = (dim==2) ? 1 : dim;
354 
355  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
356  {
357  const double dx = factor * fe.JxW(k);
358  const Point<dim> &n = fe.normal_vector(k);
359  for (unsigned int i=0; i<n_dofs; ++i)
360  for (unsigned int j=0; j<n_dofs; ++j)
361  for (unsigned int d=0; d<d_max; ++d)
362  {
363  const unsigned int d1 = (d+1)%dim;
364  const unsigned int d2 = (d+2)%dim;
365 
366  const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
367  const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
368 
369  M(i,j) += dx*u*v;
370  }
371  }
372  }
373 
390  template <int dim>
391  inline void ip_curl_matrix (
392  FullMatrix<double> &M11,
393  FullMatrix<double> &M12,
394  FullMatrix<double> &M21,
395  FullMatrix<double> &M22,
396  const FEValuesBase<dim> &fe1,
397  const FEValuesBase<dim> &fe2,
398  const double pen,
399  const double factor1 = 1.,
400  const double factor2 = -1.)
401  {
402  const unsigned int n_dofs = fe1.dofs_per_cell;
403 
404  AssertDimension(fe1.get_fe().n_components(), dim);
405  AssertDimension(fe2.get_fe().n_components(), dim);
406  AssertDimension(M11.m(), n_dofs);
407  AssertDimension(M11.n(), n_dofs);
408  AssertDimension(M12.m(), n_dofs);
409  AssertDimension(M12.n(), n_dofs);
410  AssertDimension(M21.m(), n_dofs);
411  AssertDimension(M21.n(), n_dofs);
412  AssertDimension(M22.m(), n_dofs);
413  AssertDimension(M22.n(), n_dofs);
414 
415  const double nu1 = factor1;
416  const double nu2 = (factor2 < 0) ? factor1 : factor2;
417  const double penalty = .5 * pen * (nu1 + nu2);
418 
419  // Depending on the
420  // dimension, the cross
421  // product is either a scalar
422  // (2d) or a vector
423  // (3d). Accordingly, in the
424  // latter case we have to sum
425  // up three bilinear forms,
426  // but in 2d, we don't. Thus,
427  // we need to adapt the loop
428  // over all dimensions
429  const unsigned int d_max = (dim==2) ? 1 : dim;
430 
431  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
432  {
433  const double dx = fe1.JxW(k);
434  const Point<dim> &n = fe1.normal_vector(k);
435  for (unsigned int i=0; i<n_dofs; ++i)
436  for (unsigned int j=0; j<n_dofs; ++j)
437  for (unsigned int d=0; d<d_max; ++d)
438  {
439  const unsigned int d1 = (d+1)%dim;
440  const unsigned int d2 = (d+2)%dim;
441  // curl u, curl v
442  const double cv1 = nu1*fe1.shape_grad_component(i,k,d1)[d2] - fe1.shape_grad_component(i,k,d2)[d1];
443  const double cv2 = nu2*fe2.shape_grad_component(i,k,d1)[d2] - fe2.shape_grad_component(i,k,d2)[d1];
444  const double cu1 = nu1*fe1.shape_grad_component(j,k,d1)[d2] - fe1.shape_grad_component(j,k,d2)[d1];
445  const double cu2 = nu2*fe2.shape_grad_component(j,k,d1)[d2] - fe2.shape_grad_component(j,k,d2)[d1];
446 
447  // u x n, v x n
448  const double u1= fe1.shape_value_component(j,k,d1)*n(d2) - fe1.shape_value_component(j,k,d2)*n(d1);
449  const double u2=-fe2.shape_value_component(j,k,d1)*n(d2) + fe2.shape_value_component(j,k,d2)*n(d1);
450  const double v1= fe1.shape_value_component(i,k,d1)*n(d2) - fe1.shape_value_component(i,k,d2)*n(d1);
451  const double v2=-fe2.shape_value_component(i,k,d1)*n(d2) + fe2.shape_value_component(i,k,d2)*n(d1);
452 
453  M11(i,j) += .5*dx*(2.*penalty*u1*v1 - cv1*u1 - cu1*v1);
454  M12(i,j) += .5*dx*(2.*penalty*v1*u2 - cv1*u2 - cu2*v1);
455  M21(i,j) += .5*dx*(2.*penalty*u1*v2 - cv2*u1 - cu1*v2);
456  M22(i,j) += .5*dx*(2.*penalty*u2*v2 - cv2*u2 - cu2*v2);
457  }
458  }
459  }
460 
461 
462  }
463 }
464 
465 
466 DEAL_II_NAMESPACE_CLOSE
467 
468 #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:858
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
const Point< spacedim > & normal_vector(const unsigned int i) const
Tensor< 1, dim > curl_curl(const Tensor< 2, dim > &h0, const Tensor< 2, dim > &h1, const Tensor< 2, dim > &h2)
Definition: maxwell.h:95
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:133
void tangential_trace_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double factor=1.)
Definition: maxwell.h:332
Library of integrals over cells and faces.
Definition: laplace.h:31
size_type n() const
void curl_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: maxwell.h:170
#define Assert(cond, exc)
Definition: exceptions.h:299
void curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: maxwell.h:224
size_type m() const
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
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:391
::ExceptionBase & ExcNotImplemented()
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int quadrature_point) const
void nitsche_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: maxwell.h:276