]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add trace operators
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 19 Aug 2011 13:45:18 +0000 (13:45 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 19 Aug 2011 13:45:18 +0000 (13:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@24117 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/integrators/maxwell.h

index 54a3f16999e0efbb1d2a999347173508cff75cfd..b7a6c334bafc2988ac0ad4bea8dbae4928d28b25 100644 (file)
@@ -37,7 +37,7 @@ namespace LocalIntegrators
  *   \partial_3 u_2 - \partial_2 u_3 \\
  *   \partial_1 u_3 - \partial_3 u_1 \\
  *   \partial_2 u_1 - \partial_1 u_2
- * \end{pmatrix}
+ * \end{pmatrix}.
  * @f]
  *
  * In two space dimensions, the curl is obtained by extending a vector
@@ -46,7 +46,7 @@ namespace LocalIntegrators
  * curl of a vector function and the vector curl of a scalar
  * function. The current implementation exchanges the sign and we have:
  * @f[
- *  \nabla \times \mathbf u = \partial_1 u_2 - \partial 2 u_1
+ *  \nabla \times \mathbf u = \partial_1 u_2 - \partial 2 u_1,
  *  \qquad
  *  \nabla \times p = \begin{pmatrix}
  *    \partial_2 p \\ -\partial_1 p
@@ -156,8 +156,8 @@ namespace LocalIntegrators
  * \int_Z \nabla\!\times\! u \cdot
  * \nabla\!\times\! v \,dx
  * @f]
+ * in weak form.
  *
- * @ingroup Integrators
  * @author Guido Kanschat
  * @date 2011
  */
@@ -203,14 +203,15 @@ namespace LocalIntegrators
     }
     
 /**
- * The curl operator
+ * The matrix for the curl operator
  * @f[
  * \int_Z \nabla\!\times\! u \cdot v \,dx.
  * @f]
  *
- * This is the standard curl operator in 3D and the scalar curl in 2D.
-  *
- * @ingroup Integrators
+ * This is the standard curl operator in 3D and the scalar curl in
+ * 2D. The vector curl operator can be obtained by exchanging test and
+ * trial functions.
+ *
  * @author Guido Kanschat
  * @date 2011
 */
@@ -247,6 +248,212 @@ namespace LocalIntegrators
                }
        }
     }
+
+                                    /**
+                                     * The matrix for weak boundary
+                                     * condition of Nitsche type for
+                                     * the tangential component in
+                                     * Maxwell systems.
+                                     *
+                                     * @f[
+                                     * \int_F \biggl( 2\gamma
+                                     * (u\times n) (v\times n) -
+                                     * (u\times n)(\nu \nabla\times
+                                     * v) - (v\times
+                                     * n)(\nu \nabla\times u)
+                                     * \biggr)
+                                     * @f]
+                                     *
+                                     * @author Guido Kanschat
+                                     * @date 2011
+                                     */
+      template <int dim>
+      void nitsche_curl_matrix (
+       FullMatrix<double>& M,
+       const FEValuesBase<dim>& fe,
+       double penalty,
+       double factor = 1.)
+      {
+       const unsigned int n_dofs = fe.dofs_per_cell;
+       
+       AssertDimension(fe.get_fe().n_components(), dim);
+       AssertDimension(M.m(), n_dofs);
+       AssertDimension(M.n(), n_dofs);
+       
+                                        // Depending on the
+                                        // dimension, the cross
+                                        // product is either a scalar
+                                        // (2d) or a vector
+                                        // (3d). Accordingly, in the
+                                        // latter case we have to sum
+                                        // up three bilinear forms,
+                                        // but in 2d, we don't. Thus,
+                                        // we need to adapt the loop
+                                        // over all dimensions
+       const unsigned int d_max = (dim==2) ? 1 : dim;
+       
+       for (unsigned k=0;k<fe.n_quadrature_points;++k)
+         {
+           const double dx = factor * fe.JxW(k);
+           const Point<dim>& n = fe.normal_vector(k);
+           for (unsigned i=0;i<n_dofs;++i)
+             for (unsigned j=0;j<n_dofs;++j)
+               for (unsigned int d=0;d<d_max;++d)
+                 {
+                   const unsigned int d1 = (d+1)%dim;
+                   const unsigned int d2 = (d+2)%dim;
+                   
+                   const double cv = fe.shape_grad_component(i,k,d1)[d2] - fe.shape_grad_component(i,k,d2)[d1];
+                   const double cu = fe.shape_grad_component(j,k,d1)[d2] - fe.shape_grad_component(j,k,d2)[d1];
+                   const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
+                   const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
+                   
+                   M(i,j) += dx*(2.*penalty*u*v - cv*u - cu*v);
+                 }
+         }
+      }
+                                    /**
+                                     * The product of two tangential
+                                     * traces,
+                                     * @f[
+                                     * \int_F (u\times n)(v\times n)
+                                     * \, ds.
+                                     * @f]
+                                     *
+                                     * @author Guido Kanschat
+                                     * @date 2011
+                                     */
+      template <int dim>
+      void tangential_trace_matrix (
+       FullMatrix<double>& M,
+       const FEValuesBase<dim>& fe,
+       double factor = 1.)
+      {
+       const unsigned int n_dofs = fe.dofs_per_cell;
+       
+       AssertDimension(fe.get_fe().n_components(), dim);
+       AssertDimension(M.m(), n_dofs);
+       AssertDimension(M.n(), n_dofs);
+       
+                                        // Depending on the
+                                        // dimension, the cross
+                                        // product is either a scalar
+                                        // (2d) or a vector
+                                        // (3d). Accordingly, in the
+                                        // latter case we have to sum
+                                        // up three bilinear forms,
+                                        // but in 2d, we don't. Thus,
+                                        // we need to adapt the loop
+                                        // over all dimensions
+       const unsigned int d_max = (dim==2) ? 1 : dim;
+       
+       for (unsigned k=0;k<fe.n_quadrature_points;++k)
+         {
+           const double dx = factor * fe.JxW(k);
+           const Point<dim>& n = fe.normal_vector(k);
+           for (unsigned i=0;i<n_dofs;++i)
+             for (unsigned j=0;j<n_dofs;++j)
+               for (unsigned int d=0;d<d_max;++d)
+                 {
+                   const unsigned int d1 = (d+1)%dim;
+                   const unsigned int d2 = (d+2)%dim;
+                   
+                   const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
+                   const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);
+                   
+                   M(i,j) += dx*u*v;
+                 }
+         }
+      }
+
+                                      /**
+                                       * The interior penalty fluxes
+                                       * for Maxwell systems.
+                                       *
+                                       * @f[
+                                       * \int_F \biggl( \gamma
+                                       * \{u\times n\}\{v\times n\} -
+                                       * \{u\times n\}\{\nu \nabla\times
+                                       * v\}- \{v\times
+                                       * n\}\{\nu \nabla\times u\}
+                                       * \biggr)\;dx
+                                       * @f]
+                                     *
+                                     * @author Guido Kanschat
+                                     * @date 2011
+                                       */
+      template <int dim>
+      inline 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.)
+      {
+       const unsigned int n_dofs = fe1.dofs_per_cell;
+       
+       AssertDimension(fe1.get_fe().n_components(), dim);
+       AssertDimension(fe2.get_fe().n_components(), dim);
+       AssertDimension(M11.m(), n_dofs);
+       AssertDimension(M11.n(), n_dofs);
+       AssertDimension(M12.m(), n_dofs);
+       AssertDimension(M12.n(), n_dofs);
+       AssertDimension(M21.m(), n_dofs);
+       AssertDimension(M21.n(), n_dofs);
+       AssertDimension(M22.m(), n_dofs);
+       AssertDimension(M22.n(), n_dofs);
+       
+       const double nu1 = factor1;
+       const double nu2 = (factor2 < 0) ? factor1 : factor2;
+       const double penalty = .5 * pen * (nu1 + nu2);
+
+                                        // Depending on the
+                                        // dimension, the cross
+                                        // product is either a scalar
+                                        // (2d) or a vector
+                                        // (3d). Accordingly, in the
+                                        // latter case we have to sum
+                                        // up three bilinear forms,
+                                        // but in 2d, we don't. Thus,
+                                        // we need to adapt the loop
+                                        // over all dimensions
+       const unsigned int d_max = (dim==2) ? 1 : dim;
+       
+       for (unsigned k=0;k<fe1.n_quadrature_points;++k)
+         {
+           const double dx = fe1.JxW(k);
+           const Point<dim>& n = fe1.normal_vector(k);
+           for (unsigned i=0;i<n_dofs;++i)
+             for (unsigned j=0;j<n_dofs;++j)
+               for (unsigned d=0;d<d_max;++d)
+                 {
+                   const unsigned int d1 = (d+1)%dim;
+                   const unsigned int d2 = (d+2)%dim;
+                                                    // curl u, curl v
+                   const double cv1 = nu1*fe1.shape_grad_component(i,k,d1)[d2] - fe1.shape_grad_component(i,k,d2)[d1];
+                   const double cv2 = nu2*fe2.shape_grad_component(i,k,d1)[d2] - fe2.shape_grad_component(i,k,d2)[d1];
+                   const double cu1 = nu1*fe1.shape_grad_component(j,k,d1)[d2] - fe1.shape_grad_component(j,k,d2)[d1];
+                   const double cu2 = nu2*fe2.shape_grad_component(j,k,d1)[d2] - fe2.shape_grad_component(j,k,d2)[d1];
+                   
+                                                    // u x n, v x n
+                   const double u1= fe1.shape_value_component(j,k,d1)*n(d2) - fe1.shape_value_component(j,k,d2)*n(d1);
+                   const double u2=-fe2.shape_value_component(j,k,d1)*n(d2) + fe2.shape_value_component(j,k,d2)*n(d1);
+                   const double v1= fe1.shape_value_component(i,k,d1)*n(d2) - fe1.shape_value_component(i,k,d2)*n(d1);
+                   const double v2=-fe2.shape_value_component(i,k,d1)*n(d2) + fe2.shape_value_component(i,k,d2)*n(d1);
+                   
+                   M11(i,j) += .5*dx*(2.*penalty*u1*v1 - cv1*u1 - cu1*v1);
+                   M12(i,j) += .5*dx*(2.*penalty*v1*u2 - cv1*u2 - cu2*v1);
+                   M21(i,j) += .5*dx*(2.*penalty*u1*v2 - cv2*u1 - cu1*v2);
+                   M22(i,j) += .5*dx*(2.*penalty*u2*v2 - cv2*u2 - cu2*v2);
+                 }
+         }
+      }
+      
+
   }
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.