include/deal.II/integrators/maxwell.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: maxwell.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                        // Depending on the dimension,
00177                                        // the cross product is either
00178                                        // a scalar (2d) or a vector
00179                                        // (3d). Accordingly, in the
00180                                        // latter case we have to sum
00181                                        // up three bilinear forms, but
00182                                        // in 2d, we don't. Thus, we
00183                                        // need to adapt the loop over
00184                                        // all dimensions
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                                          // Depending on the
00284                                          // dimension, the cross
00285                                          // product is either a scalar
00286                                          // (2d) or a vector
00287                                          // (3d). Accordingly, in the
00288                                          // latter case we have to sum
00289                                          // up three bilinear forms,
00290                                          // but in 2d, we don't. Thus,
00291                                          // we need to adapt the loop
00292                                          // over all dimensions
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                                          // Depending on the
00339                                          // dimension, the cross
00340                                          // product is either a scalar
00341                                          // (2d) or a vector
00342                                          // (3d). Accordingly, in the
00343                                          // latter case we have to sum
00344                                          // up three bilinear forms,
00345                                          // but in 2d, we don't. Thus,
00346                                          // we need to adapt the loop
00347                                          // over all dimensions
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                                          // Depending on the
00415                                          // dimension, the cross
00416                                          // product is either a scalar
00417                                          // (2d) or a vector
00418                                          // (3d). Accordingly, in the
00419                                          // latter case we have to sum
00420                                          // up three bilinear forms,
00421                                          // but in 2d, we don't. Thus,
00422                                          // we need to adapt the loop
00423                                          // over all dimensions
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                                                      // curl u, curl v
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                                                      // u x n, v x n
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:10 by doxygen 1.7.3