This example is to show the lowest order 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:
- Maxwell
- cubeMaxwell
- femMaxwell3
- Maxwellfemrate
The method is implemented in Maxwell
subroutine and tested in cubeMaxwell
. Together with other elements (ND0,ND1,ND2), femMaxwell3
provides a concise interface to solve Maxwell equation. The ND0 element is tested in Maxwellfemrate
. This doc is based on Maxwellfemrate
.
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];
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
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}.$$%% 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 = 'ND0';
option.printlevel = 1;
%% 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);
The optimal rate of convergence of both the H(curl)-norm and L2-norm (1st order) is observed. The 2nd order convergent rate between two discrete functions $\| \nabla \times (u_I - u_h) \|_{\infty}$ is known as superconvergence.
MGCG using HX preconditioner converges uniformly in all cases.