00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__integrators_maxwell_h
00013 #define __deal2__integrators_maxwell_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/quadrature.h>
00019 #include <deal.II/lac/full_matrix.h>
00020 #include <deal.II/fe/mapping.h>
00021 #include <deal.II/fe/fe_values.h>
00022 #include <deal.II/meshworker/dof_info.h>
00023
00024 DEAL_II_NAMESPACE_OPEN
00025
00026 namespace LocalIntegrators
00027 {
00060 namespace Maxwell
00061 {
00088 template <int dim>
00089 Tensor<1,dim>
00090 curl_curl (
00091 const Tensor<2,dim>& h0,
00092 const Tensor<2,dim>& h1,
00093 const Tensor<2,dim>& h2)
00094 {
00095 Tensor<1,dim> result;
00096 switch (dim)
00097 {
00098 case 2:
00099 result[0] = h1[0][1]-h0[1][1];
00100 result[1] = h0[0][1]-h1[0][0];
00101 break;
00102 case 3:
00103 result[0] = h1[0][1]+h2[0][2]-h0[1][1]-h0[2][2];
00104 result[1] = h2[1][2]+h0[1][0]-h1[2][2]-h1[0][0];
00105 result[2] = h0[2][0]+h1[2][1]-h2[0][0]-h2[1][1];
00106 break;
00107 default:
00108 Assert(false, ExcNotImplemented());
00109 }
00110 return result;
00111 }
00112
00126 template <int dim>
00127 Tensor<1,dim>
00128 tangential_curl (
00129 const Tensor<1,dim>& g0,
00130 const Tensor<1,dim>& g1,
00131 const Tensor<1,dim>& g2,
00132 const Tensor<1,dim>& normal)
00133 {
00134 Tensor<1,dim> result;
00135
00136 switch (dim)
00137 {
00138 case 2:
00139 result[0] = normal[1] * (g1[0]-g0[1]);
00140 result[1] =-normal[0] * (g1[0]-g0[1]);
00141 break;
00142 case 3:
00143 result[0] = normal[2]*(g2[1]-g0[2])+normal[1]*(g1[0]-g0[1]);
00144 result[1] = normal[0]*(g0[2]-g1[0])+normal[2]*(g2[1]-g1[2]);
00145 result[2] = normal[1]*(g1[0]-g2[1])+normal[0]*(g0[2]-g2[0]);
00146 break;
00147 default:
00148 Assert(false, ExcNotImplemented());
00149 }
00150 return result;
00151 }
00152
00164 template <int dim>
00165 void curl_curl_matrix (
00166 FullMatrix<double>& M,
00167 const FEValuesBase<dim>& fe,
00168 const double factor = 1.)
00169 {
00170 const unsigned int n_dofs = fe.dofs_per_cell;
00171
00172 AssertDimension(fe.get_fe().n_components(), dim);
00173 AssertDimension(M.m(), n_dofs);
00174 AssertDimension(M.n(), n_dofs);
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 const unsigned int d_max = (dim==2) ? 1 : dim;
00186
00187 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00188 {
00189 const double dx = factor * fe.JxW(k);
00190 for (unsigned i=0;i<n_dofs;++i)
00191 for (unsigned j=0;j<n_dofs;++j)
00192 for (unsigned int d=0;d<d_max;++d)
00193 {
00194 const unsigned int d1 = (d+1)%dim;
00195 const unsigned int d2 = (d+2)%dim;
00196
00197 const double cv = fe.shape_grad_component(i,k,d1)[d2] - fe.shape_grad_component(i,k,d2)[d1];
00198 const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
00199
00200 M(i,j) += dx * cu * cv;
00201 }
00202 }
00203 }
00204
00218 template <int dim>
00219 void curl_matrix (
00220 FullMatrix<double>& M,
00221 const FEValuesBase<dim>& fe,
00222 const FEValuesBase<dim>& fetest,
00223 double factor = 1.)
00224 {
00225 unsigned int t_comp = (dim==3) ? dim : 1;
00226 const unsigned int n_dofs = fe.dofs_per_cell;
00227 const unsigned int t_dofs = fetest.dofs_per_cell;
00228 AssertDimension(fe.get_fe().n_components(), dim);
00229 AssertDimension(fetest.get_fe().n_components(), t_comp);
00230 AssertDimension(M.m(), t_dofs);
00231 AssertDimension(M.n(), n_dofs);
00232
00233 const unsigned int d_max = (dim==2) ? 1 : dim;
00234
00235 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00236 {
00237 const double dx = fe.JxW(k) * factor;
00238 for (unsigned i=0;i<t_dofs;++i)
00239 for (unsigned j=0;j<n_dofs;++j)
00240 for (unsigned int d=0;d<d_max;++d)
00241 {
00242 const unsigned int d1 = (d+1)%dim;
00243 const unsigned int d2 = (d+2)%dim;
00244
00245 const double vv = fetest.shape_value_component(i,k,d);
00246 const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
00247 M(i,j) += dx * cu * vv;
00248 }
00249 }
00250 }
00251
00270 template <int dim>
00271 void nitsche_curl_matrix (
00272 FullMatrix<double>& M,
00273 const FEValuesBase<dim>& fe,
00274 double penalty,
00275 double factor = 1.)
00276 {
00277 const unsigned int n_dofs = fe.dofs_per_cell;
00278
00279 AssertDimension(fe.get_fe().n_components(), dim);
00280 AssertDimension(M.m(), n_dofs);
00281 AssertDimension(M.n(), n_dofs);
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293 const unsigned int d_max = (dim==2) ? 1 : dim;
00294
00295 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00296 {
00297 const double dx = factor * fe.JxW(k);
00298 const Point<dim>& n = fe.normal_vector(k);
00299 for (unsigned i=0;i<n_dofs;++i)
00300 for (unsigned j=0;j<n_dofs;++j)
00301 for (unsigned int d=0;d<d_max;++d)
00302 {
00303 const unsigned int d1 = (d+1)%dim;
00304 const unsigned int d2 = (d+2)%dim;
00305
00306 const double cv = fe.shape_grad_component(i,k,d1)[d2] - fe.shape_grad_component(i,k,d2)[d1];
00307 const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
00308 const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
00309 const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
00310
00311 M(i,j) += dx*(2.*penalty*u*v - cv*u - cu*v);
00312 }
00313 }
00314 }
00326 template <int dim>
00327 void tangential_trace_matrix (
00328 FullMatrix<double>& M,
00329 const FEValuesBase<dim>& fe,
00330 double factor = 1.)
00331 {
00332 const unsigned int n_dofs = fe.dofs_per_cell;
00333
00334 AssertDimension(fe.get_fe().n_components(), dim);
00335 AssertDimension(M.m(), n_dofs);
00336 AssertDimension(M.n(), n_dofs);
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348 const unsigned int d_max = (dim==2) ? 1 : dim;
00349
00350 for (unsigned k=0;k<fe.n_quadrature_points;++k)
00351 {
00352 const double dx = factor * fe.JxW(k);
00353 const Point<dim>& n = fe.normal_vector(k);
00354 for (unsigned i=0;i<n_dofs;++i)
00355 for (unsigned j=0;j<n_dofs;++j)
00356 for (unsigned int d=0;d<d_max;++d)
00357 {
00358 const unsigned int d1 = (d+1)%dim;
00359 const unsigned int d2 = (d+2)%dim;
00360
00361 const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
00362 const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
00363
00364 M(i,j) += dx*u*v;
00365 }
00366 }
00367 }
00368
00385 template <int dim>
00386 inline void ip_curl_matrix (
00387 FullMatrix<double>& M11,
00388 FullMatrix<double>& M12,
00389 FullMatrix<double>& M21,
00390 FullMatrix<double>& M22,
00391 const FEValuesBase<dim>& fe1,
00392 const FEValuesBase<dim>& fe2,
00393 const double pen,
00394 const double factor1 = 1.,
00395 const double factor2 = -1.)
00396 {
00397 const unsigned int n_dofs = fe1.dofs_per_cell;
00398
00399 AssertDimension(fe1.get_fe().n_components(), dim);
00400 AssertDimension(fe2.get_fe().n_components(), dim);
00401 AssertDimension(M11.m(), n_dofs);
00402 AssertDimension(M11.n(), n_dofs);
00403 AssertDimension(M12.m(), n_dofs);
00404 AssertDimension(M12.n(), n_dofs);
00405 AssertDimension(M21.m(), n_dofs);
00406 AssertDimension(M21.n(), n_dofs);
00407 AssertDimension(M22.m(), n_dofs);
00408 AssertDimension(M22.n(), n_dofs);
00409
00410 const double nu1 = factor1;
00411 const double nu2 = (factor2 < 0) ? factor1 : factor2;
00412 const double penalty = .5 * pen * (nu1 + nu2);
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424 const unsigned int d_max = (dim==2) ? 1 : dim;
00425
00426 for (unsigned k=0;k<fe1.n_quadrature_points;++k)
00427 {
00428 const double dx = fe1.JxW(k);
00429 const Point<dim>& n = fe1.normal_vector(k);
00430 for (unsigned i=0;i<n_dofs;++i)
00431 for (unsigned j=0;j<n_dofs;++j)
00432 for (unsigned d=0;d<d_max;++d)
00433 {
00434 const unsigned int d1 = (d+1)%dim;
00435 const unsigned int d2 = (d+2)%dim;
00436
00437 const double cv1 = nu1*fe1.shape_grad_component(i,k,d1)[d2] - fe1.shape_grad_component(i,k,d2)[d1];
00438 const double cv2 = nu2*fe2.shape_grad_component(i,k,d1)[d2] - fe2.shape_grad_component(i,k,d2)[d1];
00439 const double cu1 = nu1*fe1.shape_grad_component(j,k,d1)[d2] - fe1.shape_grad_component(j,k,d2)[d1];
00440 const double cu2 = nu2*fe2.shape_grad_component(j,k,d1)[d2] - fe2.shape_grad_component(j,k,d2)[d1];
00441
00442
00443 const double u1= fe1.shape_value_component(j,k,d1)*n(d2) - fe1.shape_value_component(j,k,d2)*n(d1);
00444 const double u2=-fe2.shape_value_component(j,k,d1)*n(d2) + fe2.shape_value_component(j,k,d2)*n(d1);
00445 const double v1= fe1.shape_value_component(i,k,d1)*n(d2) - fe1.shape_value_component(i,k,d2)*n(d1);
00446 const double v2=-fe2.shape_value_component(i,k,d1)*n(d2) + fe2.shape_value_component(i,k,d2)*n(d1);
00447
00448 M11(i,j) += .5*dx*(2.*penalty*u1*v1 - cv1*u1 - cu1*v1);
00449 M12(i,j) += .5*dx*(2.*penalty*v1*u2 - cv1*u2 - cu2*v1);
00450 M21(i,j) += .5*dx*(2.*penalty*u1*v2 - cv2*u1 - cu1*v2);
00451 M22(i,j) += .5*dx*(2.*penalty*u2*v2 - cv2*u2 - cu2*v2);
00452 }
00453 }
00454 }
00455
00456
00457 }
00458 }
00459
00460
00461 DEAL_II_NAMESPACE_CLOSE
00462
00463 #endif
00464