Skip to content
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
109 lines (77 sloc) 4.28 KB
// This code simulates the electromagnetic wave radiation of a half-wave dipole antenna excited at 1 GHz.
// The simulation is performed in 2D (i.e. the fields are constant in the z direction). The goal of this code
// is to provide a lightweight example that can be used with minimal modifications to simulate actual 3D antennas.
// The radiation pattern of this dipole antenna can be obtained by calling 'interpolate' on the time-averaged power
// expression for the points of a circle near the domain boundary and centered at the origin. The interpolated
// values can then be saved to .csv using 'writevector' and 'polarplot' could be used in Matlab for visualization.
#include "sparselizardbase.h"
using namespace mathop;
void sparselizard(void)
wallclock clk;
// The domain regions as defined in 'dipoleantenna.geo':
int air = 1, conductor = 2, feed = 3, boundary = 4;
mesh mymesh("dipoleantenna.msh");
// Define the whole non-conducting region for convenience:
int wholedomain = regionunion({air, feed});
// Edge shape functions 'hcurl' for the electric field E.
// Because of the propagating EM waves the electric field has an
// in-phase and quadrature component (thus harmonics 2 and 3 are used):
// E = E2 * sin(2*pi*1GHz*t) + E3 * cos(2*pi*1GHz*t).
field E("hcurl", {2,3});
field Es = E.harmonic(2), Ec = E.harmonic(3);
// Use interpolation order 2 on the whole domain:
E.setorder(wholedomain, 2);
// Define the speed of light [m/s] and the working frequency [Hz].
// The total half-wave dipole antenna length is 15 cm. This corresponds to a 1GHz frequency:
double c = 3e8, freq = 1e9;
// The dipole is a perfect conductor. We thus force E to 0 on it:
// We force an electric field of 1 V/m in the y direction on the feed port
// for the sin component of E only (the cos component is forced to zero):
Es.setconstraint(feed, array3x1(0,1,0));
formulation maxwell;
// This is the weak formulation in time for electromagnetic waves:
maxwell += integral(wholedomain, -curl(dof(E))*curl(tf(E)) - 1/(c*c)*dtdt(dof(E))*tf(E));
// The OUTWARD pointing normal is required for the wave radiation condition.
// Depending on the 'boundary' lines orientation the normal can be flipped.
// To confirm the normal direction it is written to disk below:
expression n = -normal(boundary);
n.write(boundary, "normal.pos");
// Silver-Mueller radiation condition to force outgoing electromagnetic waves:
maxwell += integral(boundary, -1/c * crossproduct(crossproduct(n, dt(dof(E))), n) * tf(E));
// Generate the algebraic matrices A and right handside b of Ax = b:
// Get the solution vector x of Ax = b:
vec solE = solve(maxwell.A(), maxwell.b());
// Transfer the data in the solution vector to field E:
E.setdata(wholedomain, solE);
clk.print("Resolution time (before post-processing):");
// Save the electric field E and magnetic field H with an order 2 interpolation.
// Since from Maxwell -mu0*dt(H) = curl(E) we have that H = 1/(mu0*w^2) * curl(dt(E)).
double mu0 = 4*getpi()*1e-7;
expression H = 1/(mu0*pow(2*getpi()*freq, 2)) * curl(dt(E));
H.write(wholedomain, "H.pos", 2);
E.write(wholedomain, "E.pos", 2);
// Write the electric field at 50 timesteps of a period for a time visualization:
E.write(wholedomain, "E.pos", 2, 50);
// Write the Poynting vector E x H.
expression S = crossproduct(E, H);
// This involves a product that makes new harmonics appear at 2x1GHz (harmonics 4 and 5) plus a
// constant component (harmonic 1). An FFT is thus required (performed on 6 timesteps below):
S.write(wholedomain, 6, "S.pos", 2);
// Also save in time at 50 timesteps of a period:
S.write(wholedomain, "S.pos", 2, 50);
// Code validation line. Can be removed.
std::cout << (norm(crossproduct(Es,Ec)).max(wholedomain, 5)[0] < 0.00218416 && norm(crossproduct(Es,Ec)).max(wholedomain, 5)[0] > 0.00218413);
int main(void)
return 0;
You can’t perform that action at this time.