2D FEM Solution of the Helmholtz Equation (TMz)^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Description:2D solution of the Helmholtz wave equation for the TMz polarization. Arectangular domain is assumed. The side walls are zero Neumann boundaries.The bottom wall is a Dirichlet boundary. The top wall is a Neumann boundary.This version assumes a homogeneous media. The domain is discretized into Ntriangular elements. First order Nodal elements are assumed. The results areplotted and the error is estimated against the exact solution.^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^f300 106⋅16:= c0 2.998 108⋅:= λc0f:= λ 15.989=VDBC01:=k02 π⋅λ:=Read in the coordinates of the nodes from an external file:rNodenodeList.txt:= triNodeelemList.txt:=numNodes rows rNode():= numTris rows triNode():=numNodes 36=numTris 50=rN rNodeT:=Graph the Mesh:meshMapx index 0←node triNodeij,1−←tmpindexrNodenode 0,←index index 1+←j02..∈fori 0 numTris 1−..∈fortmp:= meshMapy index 0←node triNodeij,1−←tmpindexrNodenode 1,←index index 1+←j02..∈fori 0 numTris 1−..∈fortmp:=0 0.5 100.51meshMapymeshMapxRead in the boundary nodes:The node # (1st column) and the Dirichlet boundary ID (2nd Column) are read in.Note that a pre-processor has sorted all nodes so that nodes on Dirichlet boundaries are at theend of the list. boundaryNodesboundaryNodeList.txt:= numBoundaryNodes rows boundaryNodes():=numBoundaryNodes 12=DirchletboundaryNodes index 0←tmpindex 0,boundaryNodesi0,← boundaryNodesi1,1=iftmpindex 1,boundaryNodesi1,← boundaryNodesi1,1=ifindex index 1+← boundaryNodesi1,1=ifi 0 numBoundaryNodes 1−..∈fortmp:=numDBoundaryNodes rows DirchletboundaryNodes():= numDBoundaryNodes 6=numIntNodes numNodes numDBoundaryNodes−:= numIntNodes 30=nodeBctmpi0←i 0 numNodes 1−..∈fornode DirchletboundaryNodesi0,1−←bc DirchletboundaryNodesi1,←tmpnodebc←i 0 numDBoundaryNodes 1−..∈fortmp:=NeumannNodes index 0←tmpindexboundaryNodesi0,← boundaryNodesi1,2=ifindex index 1+← boundaryNodesi1,2=ifi 0 numBoundaryNodes 1−..∈fortmp:=Compute the unitary vectors for the triangles:ii 0 numTris..:=aSub1Vectors index 0←node0 triNodei0,1−←node1 triNodei1,1−←asub1 rNnode1〈〉rNnode0〈〉−←tmpji,asub1j←j02..∈fori 0 numTris 1−..∈fortmp:= aSub2Vectors index 0←node0 triNodei0,1−←node2 triNodei2,1−←asub2 rNnode2〈〉rNnode0〈〉−←tmpji,asub2j←j02..∈fori 0 numTris 1−..∈fortmp:=rootG2d index 0←a1 aSub1Vectorsi〈〉←a2 aSub2Vectorsi〈〉←rg2dia1 a2×←i 0 numTris 1−..∈forrg2d:=Compute the reciprocal unitary vectors for the triangles:zhat001⎛⎜⎜⎝⎞⎠:=aSup1Vectors index 0←asub2 aSub2Vectorsi〈〉←asup1asub2 zhat×rootG2di←tmpji,asup1j←j02..∈fori 0 numTris 1−..∈fortmp:= aSup2Vectors index 0←asub1 aSub1Vectorsi〈〉←asup2zhat asub1×rootG2di←tmpji,asup2j←j02..∈fori 0 numTris 1−..∈fortmp:=Compute the Global K-matrix:Ktmpij,0←i 0 numNodes 1−..∈forj 0 numNodes 1−..∈fornodejtriNodeij,1−←j02..∈forasup1 aSup1Vectorsi〈〉←asup2 aSup2Vectorsi〈〉←asup0 asup1− asup2−←asupj0,asup0j←asupj1,asup1j←asupj2,asup2j←j02..∈foraq asupq〈〉←as asups〈〉←Keaq as⋅2Teqs,−⎛⎜⎝⎞⎠rootG2di⋅←tmpnodeq()nodes,tmpnodeq()nodes,Ke+← nodesnumIntNodes<ifs02..∈for nodeqnumIntNodes<iftmpnodeq()nodeq,1← otherwiseq02..∈forblankspace 0←i 0 numTris 1−..∈fortmp:=K0123401234567891011121.998 -1.001 0 0 0-1.001 3.997 -1.001 0 00 -1.001 3.997 -1.001 00 0 -1.001 3.997 -1.0010 0 0 -1.001 3.9970 0 0 0 -1.001-0.5-4-5.147·10 0 0 00 -1.001-4-5.147·10 0 00 0 -1.001-4-5.147·10 00 0 0 -1.001-4-5.147·100 0 0 0 -1.0010000000000=Tek0212112121211212121⎛⎜⎜⎜⎜⎜⎜⎝⎞⎟⎟⎟⎟⎠⋅:=Compute the Right-Hand-Side:brhsj0←j 0 numNodes 1−..∈fornodejtriNodeij,1−←bcjnodeBcnodej1−←j02..∈forasup1 aSup1Vectorsi〈〉←asup2 aSup2Vectorsi〈〉←asup0 asup1− asup2−←asupj0,asup0j←asupj1,asup1j←asupj2,asup2j←j02..∈foraq asupq〈〉←as asups〈〉←Keaq as⋅2Teqs,−⎛⎜⎝⎞⎠rootG2di⋅←rhsnodeqrhsnodeqKe VDBCbcs⋅−← nodesnumIntNodes≥ifs02..∈for nodeqnumIntNodes<ifrhsnodeqVDBCbcq← otherwiseq02..∈forblankspace 0←i 0 numTris 1−..∈forrhs:=Compute the Solution:b001234567891011121314150.5011.0011.0011.0011.0010.50000000000=c K1−b⋅:=Compute the Error Relative to the Exact Solution:ymax 1:=ExactyrN1i,←tmpicos k0 y ymax−()⋅[]cos k0 ymax⋅()←i 0 numNodes 1−..∈fortmp:=Exact001234567891011121314151.0294577781951.0294577781951.0294577781951.0294577781951.0294577781951.0294577781951.0525601401491.0525601401491.0525601401491.0525601401491.0525601401491.0525601401491.0691644621151.0691644621151.0691644621151.069164462115= c001234567891011121314151.0294720981391.029470855861.0294532785351.0294305997461.0294130009841.0294117193181.0526214997981.0526058002961.0525595105681.0525041119881.0524577964351.0524420296861.0693277009421.0692763203671.0691801131091.069074146415=Error Exact c−:=meanError1numIntNodes0numIntNodes 1−iErrori∑=⋅:=meanError 1.968