Reference documentation for deal.II version Git 01fb814 2017-08-18 21:00:44 -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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/quadrature.h>
23 #include <deal.II/lac/full_matrix.h>
24 #include <deal.II/fe/mapping.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/meshworker/dof_info.h>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 namespace LocalIntegrators
31 {
63  namespace Maxwell
64  {
91  template <int dim>
94  const Tensor<2,dim> &h0,
95  const Tensor<2,dim> &h1,
96  const Tensor<2,dim> &h2)
97  {
98  Tensor<1,dim> result;
99  switch (dim)
100  {
101  case 2:
102  result[0] = h1[0][1]-h0[1][1];
103  result[1] = h0[0][1]-h1[0][0];
104  break;
105  case 3:
106  result[0] = h1[0][1]+h2[0][2]-h0[1][1]-h0[2][2];
107  result[1] = h2[1][2]+h0[1][0]-h1[2][2]-h1[0][0];
108  result[2] = h0[2][0]+h1[2][1]-h2[0][0]-h2[1][1];
109  break;
110  default:
111  Assert(false, ExcNotImplemented());
112  }
113  return result;
114  }
115 
129  template <int dim>
132  const Tensor<1,dim> &g0,
133  const Tensor<1,dim> &g1,
134  const Tensor<1,dim> &g2,
135  const Tensor<1,dim> &normal)
136  {
137  Tensor<1,dim> result;
138 
139  switch (dim)
140  {
141  case 2:
142  result[0] = normal[1] * (g1[0]-g0[1]);
143  result[1] =-normal[0] * (g1[0]-g0[1]);
144  break;
145  case 3:
146  result[0] = normal[2]*(g2[1]-g0[2])+normal[1]*(g1[0]-g0[1]);
147  result[1] = normal[0]*(g0[2]-g1[0])+normal[2]*(g2[1]-g1[2]);
148  result[2] = normal[1]*(g1[0]-g2[1])+normal[0]*(g0[2]-g2[0]);
149  break;
150  default:
151  Assert(false, ExcNotImplemented());
152  }
153  return result;
154  }
155 
167  template <int dim>
170  const FEValuesBase<dim> &fe,
171  const double factor = 1.)
172  {
173  const unsigned int n_dofs = fe.dofs_per_cell;
174 
175  AssertDimension(fe.get_fe().n_components(), dim);
176  AssertDimension(M.m(), n_dofs);
177  AssertDimension(M.n(), n_dofs);
178 
179  // Depending on the dimension,
180  // the cross product is either
181  // a scalar (2d) or a vector
182  // (3d). Accordingly, in the
183  // latter case we have to sum
184  // up three bilinear forms, but
185  // in 2d, we don't. Thus, we
186  // need to adapt the loop over
187  // all dimensions
188  const unsigned int d_max = (dim==2) ? 1 : dim;
189 
190  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
191  {
192  const double dx = factor * fe.JxW(k);
193  for (unsigned int i=0; i<n_dofs; ++i)
194  for (unsigned int j=0; j<n_dofs; ++j)
195  for (unsigned int d=0; d<d_max; ++d)
196  {
197  const unsigned int d1 = (d+1)%dim;
198  const unsigned int d2 = (d+2)%dim;
199 
200  const double cv = fe.shape_grad_component(i,k,d2)[d1] - fe.shape_grad_component(i,k,d1)[d2];
201  const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];
202 
203  M(i,j) += dx * cu * cv;
204  }
205  }
206  }
207 
221  template <int dim>
222  void curl_matrix (
224  const FEValuesBase<dim> &fe,
225  const FEValuesBase<dim> &fetest,
226  double factor = 1.)
227  {
228  const unsigned int n_dofs = fe.dofs_per_cell;
229  const unsigned int t_dofs = fetest.dofs_per_cell;
230  AssertDimension(fe.get_fe().n_components(), dim);
231  // There should be the right number of components (3 in 3D, otherwise 1)
232  // for the curl.
233  AssertDimension(fetest.get_fe().n_components(), (dim == 3) ? dim : 1);
234  AssertDimension(M.m(), t_dofs);
235  AssertDimension(M.n(), n_dofs);
236 
237  const unsigned int d_max = (dim==2) ? 1 : dim;
238 
239  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
240  {
241  const double dx = fe.JxW(k) * factor;
242  for (unsigned int i=0; i<t_dofs; ++i)
243  for (unsigned int j=0; j<n_dofs; ++j)
244  for (unsigned int d=0; d<d_max; ++d)
245  {
246  const unsigned int d1 = (d+1)%dim;
247  const unsigned int d2 = (d+2)%dim;
248 
249  const double vv = fetest.shape_value_component(i,k,d);
250  const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];
251  M(i,j) += dx * cu * vv;
252  }
253  }
254  }
255 
272  template <int dim>
275  const FEValuesBase<dim> &fe,
276  const unsigned int face_no,
277  double penalty,
278  double factor = 1.)
279  {
280  const unsigned int n_dofs = fe.dofs_per_cell;
281 
282  AssertDimension(fe.get_fe().n_components(), dim);
283  AssertDimension(M.m(), n_dofs);
284  AssertDimension(M.n(), n_dofs);
285 
286  // Depending on the
287  // dimension, the cross
288  // product is either a scalar
289  // (2d) or a vector
290  // (3d). Accordingly, in the
291  // latter case we have to sum
292  // up three bilinear forms,
293  // but in 2d, we don't. Thus,
294  // we need to adapt the loop
295  // over all dimensions
296  const unsigned int d_max = (dim==2) ? 1 : dim;
297 
298  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
299  {
300  const double dx = factor * fe.JxW(k);
301  const Tensor<1,dim> n = fe.normal_vector(k);
302  for (unsigned int i=0; i<n_dofs; ++i)
303  for (unsigned int j=0; j<n_dofs; ++j)
304  if (fe.get_fe().has_support_on_face(i, face_no) && fe.get_fe().has_support_on_face(j, face_no))
305  {
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,d2)[d1] - fe.shape_grad_component(i,k,d1)[d2];
312  const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];
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  }
320  }
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 Tensor<1,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 
389  template <int dim>
390  inline void ip_curl_matrix (
391  FullMatrix<double> &M11,
392  FullMatrix<double> &M12,
393  FullMatrix<double> &M21,
394  FullMatrix<double> &M22,
395  const FEValuesBase<dim> &fe1,
396  const FEValuesBase<dim> &fe2,
397  const double pen,
398  const double factor1 = 1.,
399  const double factor2 = -1.)
400  {
401  const unsigned int n_dofs = fe1.dofs_per_cell;
402 
403  AssertDimension(fe1.get_fe().n_components(), dim);
404  AssertDimension(fe2.get_fe().n_components(), dim);
405  AssertDimension(M11.m(), n_dofs);
406  AssertDimension(M11.n(), n_dofs);
407  AssertDimension(M12.m(), n_dofs);
408  AssertDimension(M12.n(), n_dofs);
409  AssertDimension(M21.m(), n_dofs);
410  AssertDimension(M21.n(), n_dofs);
411  AssertDimension(M22.m(), n_dofs);
412  AssertDimension(M22.n(), n_dofs);
413 
414  const double nu1 = factor1;
415  const double nu2 = (factor2 < 0) ? factor1 : factor2;
416  const double penalty = .5 * pen * (nu1 + nu2);
417 
418  // Depending on the
419  // dimension, the cross
420  // product is either a scalar
421  // (2d) or a vector
422  // (3d). Accordingly, in the
423  // latter case we have to sum
424  // up three bilinear forms,
425  // but in 2d, we don't. Thus,
426  // we need to adapt the loop
427  // over all dimensions
428  const unsigned int d_max = (dim==2) ? 1 : dim;
429 
430  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
431  {
432  const double dx = fe1.JxW(k);
433  const Tensor<1,dim> n = fe1.normal_vector(k);
434  for (unsigned int i=0; i<n_dofs; ++i)
435  for (unsigned int j=0; j<n_dofs; ++j)
436  for (unsigned int d=0; d<d_max; ++d)
437  {
438  const unsigned int d1 = (d+1)%dim;
439  const unsigned int d2 = (d+2)%dim;
440  // curl u, curl v
441  const double cv1 = nu1*fe1.shape_grad_component(i,k,d2)[d1] - fe1.shape_grad_component(i,k,d1)[d2];
442  const double cv2 = nu2*fe2.shape_grad_component(i,k,d2)[d1] - fe2.shape_grad_component(i,k,d1)[d2];
443  const double cu1 = nu1*fe1.shape_grad_component(j,k,d2)[d1] - fe1.shape_grad_component(j,k,d1)[d2];
444  const double cu2 = nu2*fe2.shape_grad_component(j,k,d2)[d1] - fe2.shape_grad_component(j,k,d1)[d2];
445 
446  // u x n, v x n
447  const double u1= fe1.shape_value_component(j,k,d1)*n(d2) - fe1.shape_value_component(j,k,d2)*n(d1);
448  const double u2=-fe2.shape_value_component(j,k,d1)*n(d2) + fe2.shape_value_component(j,k,d2)*n(d1);
449  const double v1= fe1.shape_value_component(i,k,d1)*n(d2) - fe1.shape_value_component(i,k,d2)*n(d1);
450  const double v2=-fe2.shape_value_component(i,k,d1)*n(d2) + fe2.shape_value_component(i,k,d2)*n(d1);
451 
452  M11(i,j) += .5*dx*(2.*penalty*u1*v1 - cv1*u1 - cu1*v1);
453  M12(i,j) += .5*dx*(2.*penalty*v1*u2 - cv1*u2 - cu2*v1);
454  M21(i,j) += .5*dx*(2.*penalty*u1*v2 - cv2*u1 - cu1*v2);
455  M22(i,j) += .5*dx*(2.*penalty*u2*v2 - cv2*u2 - cu2*v2);
456  }
457  }
458  }
459 
460 
461  }
462 }
463 
464 
465 DEAL_II_NAMESPACE_CLOSE
466 
467 #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:1166
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
const unsigned int dofs_per_cell
Definition: fe_values.h:1786
Tensor< 1, dim > curl_curl(const Tensor< 2, dim > &h0, const Tensor< 2, dim > &h1, const Tensor< 2, dim > &h2)
Definition: maxwell.h:93
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:131
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: advection.h:30
size_type n() const
void curl_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: maxwell.h:168
#define Assert(cond, exc)
Definition: exceptions.h:328
void curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, double factor=1.)
Definition: maxwell.h:222
size_type m() const
const unsigned int n_quadrature_points
Definition: fe_values.h:1779
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:390
void nitsche_curl_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const unsigned int face_no, double penalty, double factor=1.)
Definition: maxwell.h:273
static::ExceptionBase & ExcNotImplemented()
const FiniteElement< dim, spacedim > & get_fe() const
double JxW(const unsigned int quadrature_point) const