This example is to show the quadratic edge element approximation of the electric field of the time harmonic Maxwell equation.
\begin{align} \nabla \times (\mu^{-1}\nabla \times u) - \omega^2 \varepsilon \, u &= J \quad \text{ in } \quad \Omega, \\ n \times u &= n \times g_D \quad \text{ on } \quad \Gamma_D,\\ n \times (\mu^{-1}\nabla \times u) &= n \times g_N \quad \text{ on } \quad \Gamma_N. \end{align}based on the weak formulation
$$(\mu^{-1}\nabla \times u, \nabla \times v) - (\omega^2\varepsilon u,v) = (J,v) - \langle n \times g_N,v \rangle_{\Gamma_N}.$$Reference
Subroutines:
- Maxwell2
- cubeMaxwell2
- femMaxwell3
- Maxwell2femrate
The method is implemented in Maxwell2
subroutine and tested in cubeMaxwell2
. Together with other elements (ND0,ND1,ND2), femMaxwell3
provides a concise interface to solve Maxwell equation. The ND1 element is tested in Maxwell2femrate
. This doc is based on Maxwell2femrate
.
Locally we construct locBasesIdx
to record the local index used in the bases. Globally we use ascend ordering for each element and thus the orientation of the edge is consistent. No need of elem2edgeSign
. Read Simplicial complex in three dimensions for more discussion of indexing, ordering and orientation.
In addition to the edge structure, we need face and the corresponding pointers.
[node,elem] = cubemesh([0,1,0,1,0,1],1);
[elem2edge,edge] = dof3edge(elem);
[elem2face,face] = dof3face(elem);
Furthermore we need pointer from face to edge. face
is given by auxtructure3
and is sorted according to global indices.
face2edge
is used to compute uI
. So it is consistent with the local index system in edgeinterpolate2
, i.e., if face is (i,j,k)
with i<j<k
, then the
three edges are [i j], [i k], [j k]
.
face2edge = zeros(size(face,1),3,'int32');
face2edge(elem2face(:,1),:) = elem2edge(:,[4 5 6]);
face2edge(elem2face(:,2),:) = elem2edge(:,[2 3 6]);
face2edge(elem2face(:,3),:) = elem2edge(:,[1 3 5]);
face2edge(elem2face(:,4),:) = elem2edge(:,[1 2 4]);
locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];
locBasesIdx = [1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % phi
1 2 0; 1 3 0; 1 4 0; 2 3 0; 2 4 0; 3 4 0; ... % psi
3 2 4; 3 1 4; 2 1 4; 2 1 3; ...
4 2 3; 4 1 3; 4 1 2; 3 1 2]; % chi
Locally we construct locBasesIdx
to record the local index used in the bases. For example, for basis $\chi_i =\lambda_{i_1}\phi _{i_2i_3}$ for i=4
, we can get i1,i2,i3
by:
i = 4+12;
i1 = locBasesIdx(i,1); i2 = locBasesIdx(i,2); i3 = locBasesIdx(i,3);
Suppose [i,j]
is the kth edge and i<j
. The basis is given by
Inside one tetrahedron, the 6 bases functions along with their curl
corresponding to 6 local edges [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]
are
The additional 6 bases for the second family are:
$$ \psi_k = \lambda_i\nabla \lambda_j + \lambda_j \nabla \lambda_i,\qquad \nabla \times \psi_k = 0.$$$$ \psi_1 = \lambda_1\nabla\lambda_2 + \lambda_2\nabla\lambda_1,$$$$ \psi_2 = \lambda_1\nabla\lambda_3 + \lambda_3\nabla\lambda_1,$$$$ \psi_3 = \lambda_1\nabla\lambda_4 + \lambda_4\nabla\lambda_1,$$$$ \psi_4 = \lambda_2\nabla\lambda_3 + \lambda_3\nabla\lambda_2,$$$$ \psi_5 = \lambda_2\nabla\lambda_4 + \lambda_4\nabla\lambda_2,$$$$ \psi_6 = \lambda_3\nabla\lambda_4 + \lambda_4\nabla\lambda_3.$$Suppose i,j,k
are the vertices of the l
-th face and i<j<k
. The two
basis associated to this face are
Inside one tetrahedron, the 8 bases functions assocaited to the four
local faces [2 3 4; 1 3 4; 1 2 4; 1 2 3]
are:
Reference: See page 12, Table 9.2. Arnold, Douglas N. and Falk, Richard S. and Winther, Ragnar. Geometric decompositions and local bases for spaces of finite element differential forms. Comput. Methods Appl. Mech. Engrg. 198():1660--1672, 2009.
Locally, we order the local bases in the following way:
$$\{\chi_1^1,~\,\chi_1^2,~\,\chi_2^1,~\,\chi_2^2,~\,\chi_3^1,~\,\chi_3^2, ~\,\chi_4^1,~\,\chi_4^2.\}$$and rewrite the local bases as:
13- 20: $$\{\chi_1,~\,\chi_2,~\,\chi_3,~\,\chi_4,~\,\chi_5,~\,\chi_6,~\, \chi_7,~\,\chi_8.\}$$
Suppose [i,j]
is the kth edge and i<j
. The corresponding degree of freedom is
It is dual to the basis $\{\phi_k\}$ in the sense that
$$l_{\ell}(\phi _k) = \delta_{k,\ell}.$$The 6 degree of freedoms for $\psi_k$ are:
$$l_k^1 (v) = 3\int_{e_k} v\cdot t(\lambda _i - \lambda_j) \, {\rm d}s \approx \frac{1}{2}[v(i) - v(j)]\cdot e_{k}.$$The 8 degree of freedoms for $\chi_i$ is given in edgeinterpolate2
.
%% Setting
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
mesh = struct('node',node,'elem',elem);
option.L0 = 0;
option.maxIt = 4;
option.elemType = 'ND2';
option.printlevel = 1;
option.plotflag = 1;
imatlab_export_fig('print-png') % Static png figures.
%% Dirichlet boundary condition.
fprintf('Dirichlet boundary conditions. \n');
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Dirichlet');
femMaxwell3(mesh,pde,option);
fprintf('Neumann boundary condition. \n');
option.plotflag = 0;
pde = Maxwelldata2;
mesh.bdFlag = setboundary3(node,elem,'Neumann');
femMaxwell3(mesh,pde,option);
Both the H(curl)-norm and the L2-norm is 2nd order.
MGCG using HX preconditioner converges uniformly in all cases.