This example is to show the linear 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:
- Maxwell1
- cubeMaxwell1
- femMaxwell3
- Maxwell1femrate
The method is implemented in Maxwell1
subroutine and tested in cubeMaxwell1
. Together with other elements (ND0,ND1,ND2), femMaxwell3
provides a concise interface to solve Maxwell equation. The ND1 element is tested in Maxwell1femrate
. This doc is based on Maxwell1femrate
.
Use the function
[elem2dof,edge,elem2edgeSign] = dof3edge(elem);
to construct the pointer from element index to edge index. Read
<dof3edgedoc.html Dof on Edges in Three Dimensions> for details.
node = [1,0,0; 0,1,0; 0,0,0; 0,0,1];
elem = [1 2 3 4];
localEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
imatlab_export_fig('print-png') % Static png figures.
set(gcf,'Units','normal');
set(gcf,'Position',[0.25,0.25,0.25,0.25]);
showmesh3(node,elem);
view(-72,9);
findnode3(node);
findedge(node,localEdge,'all','vec');
The six dofs associated to edges in a tetrahedron is sorted in the ordering [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]
. Here [1 2 3 4]
are local indices of vertices.
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.
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]
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 additional 6 degree of freedoms 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}.$$%% Setting
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
mesh = struct('node',node,'elem',elem);
option.L0 = 1;
option.maxIt = 4;
option.elemType = 'ND1';
option.printlevel = 1;
%% Dirichlet boundary condition.
fprintf('Dirichlet boundary conditions. \n');
pde = Maxwelldata2;
bdFlag = setboundary3(node,elem,'Dirichlet');
femMaxwell3(mesh,pde,option);
%% Pure Neumann boundary condition.
fprintf('Neumann boundary condition. \n');
option.plotflag = 0;
pde = Maxwelldata2;
mesh.bdFlag = setboundary3(node,elem,'Neumann');
femMaxwell3(mesh,pde,option);
The H(curl)-norm is still 1st order but the L2-norm is improved to 2nd order.
MGCG using HX preconditioner converges uniformly in all cases.