main.cpp 2 KB
Newer Older
Alexandre Halbach committed
1 2 3 4
// This code simulates the steady-state electromagnetic waves in a cross-shaped 
// 3D waveguide made of a perfect conductor.


5
#include "sparselizardbase.h"
6 7 8 9 10 11


using namespace mathop;

void sparselizard(void)
{	
Alexandre Halbach committed
12
    // The domain regions as defined in 'waveguide3D.geo':
13 14 15 16 17
    int left = 1, skin = 2, wholedomain = 3;

    mesh mymesh("waveguide3D.msh");

    // Edge shape functions 'hcurl' for the electric field E.
Alexandre Halbach committed
18
    // Fields x, y and z are the x, y and z coordinate fields.
19 20
    field E("hcurl"), x("x"), y("y"), z("z");

Alexandre Halbach committed
21
    // Use interpolation order 2 on the whole domain:
22 23
    E.setorder(wholedomain, 2);
    
Alexandre Halbach committed
24
    // The cutoff frequency for a 0.2x0.2 m^2 cross section is freq = 1.06 GHz in theory. 
Alexandre Halbach committed
25
    // With this code and a fine enough mesh you will get the same value.
26 27
    double freq = 1.2e9, c = 3e8, pi = 3.14159, k = 2*pi*freq/c;
    
Alexandre Halbach committed
28
    // The waveguide is a perfect conductor. We thus force all
29 30 31 32 33 34 35
    // 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;
Alexandre Halbach committed
36 37
    
    // This is the weak formulation for electromagnetic waves:
38 39
    maxwell += integral(wholedomain, -curl(dof(E))*curl(tf(E)) + k*k*dof(E)*tf(E));
    
Alexandre Halbach committed
40
    // The operations below take about a minute on my laptop... be patient!
41 42 43
    maxwell.generate();
    vec solE = solve(maxwell.A(), maxwell.b());
    
Alexandre Halbach committed
44
    E.setdata(wholedomain, solE);    
45 46 47
    // Save the electric field E and magnetic field H with an order 2 interpolation:
    curl(E).write(wholedomain, "H.pos", 2);
    E.write(wholedomain, "E.pos", 2);
48 49
    
    // Code validation line. Can be removed.
Alexandre Halbach committed
50
    std::cout << ((E*curl(E)).integrate(wholedomain, 5) < -5.6255e-08 && (E*curl(E)).integrate(wholedomain, 5) > -5.6256e-08);
51 52 53 54
}

int main(void)
{	
Alexandre Halbach committed
55
    SlepcInitialize(0,{},0,0);
56 57 58

    sparselizard();

Alexandre Halbach committed
59
    SlepcFinalize();
60 61 62 63 64 65 66 67 68 69 70 71 72

    return 0;
}