* \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
* 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
* \int_Z \nabla\!\times\! u \cdot
* \nabla\!\times\! v \,dx
* @f]
+ * in weak form.
*
- * @ingroup Integrators
* @author Guido Kanschat
* @date 2011
*/
}
/**
- * 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
*/
}
}
}
+
+ /**
+ * 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);
+ }
+ }
+ }
+
+
}
}