Skip to content
Permalink
Fetching contributors…
Cannot retrieve contributors at this time
64 lines (43 sloc) 1.98 KB
// This code simulates the steady-state electromagnetic waves in a cross-shaped
// 3D waveguide made of a perfect conductor.
#include "sparselizardbase.h"
using namespace mathop;
void sparselizard(void)
{
// The domain regions as defined in 'waveguide3D.geo':
int left = 1, skin = 2, wholedomain = 3;
mesh mymesh("waveguide3D.msh");
// Edge shape functions 'hcurl' for the electric field E.
// Fields x, y and z are the x, y and z coordinate fields.
field E("hcurl"), x("x"), y("y"), z("z");
// Use interpolation order 2 on the whole domain:
E.setorder(wholedomain, 2);
// The cutoff frequency for a 0.2x0.2 m^2 cross section is freq = 1.06 GHz in theory.
// With this code and a fine enough mesh you will get the same value.
double freq = 1.2e9, c = 3e8, pi = 3.14159, k = 2*pi*freq/c;
// The waveguide is a perfect conductor. We thus force all
// tangential components of E to 0 on the waveguide skin.
E.setconstraint(skin);
// We force an electric field in the z direction on region 'left'
// that is 0 on the exterior of 'left' and 1 in the center.
E.setconstraint(left, cos(y/0.2*pi)* cos(z/0.2*pi)* array3x1(0,0,1));
formulation maxwell;
// This is the weak formulation for electromagnetic waves:
maxwell += integral(wholedomain, -curl(dof(E))*curl(tf(E)) + k*k*dof(E)*tf(E));
// Generate and solve the algebraic system Ax = b:
maxwell.generate();
vec solE = solve(maxwell.A(), maxwell.b());
// Transfer the data from the solution vector to field E:
E.setdata(wholedomain, solE);
// Save the electric field E with an order 2 interpolation:
E.write(wholedomain, "E.pos", 2);
// Code validation line. Can be removed.
std::cout << ((abs(E)*curl(E)).integrate(wholedomain, 5) < 6.9398e-07 && (abs(E)*curl(E)).integrate(wholedomain, 5) > 6.93978e-07);
}
int main(void)
{
SlepcInitialize(0,{},0,0);
sparselizard();
SlepcFinalize();
return 0;
}
You can’t perform that action at this time.