Maxwell (EFIE)#

Goal#

Solve a Maxwell electric field integral equation (EFIE) for PEC scattering.

This tutorial corresponds to the example program:

  • examples/maxwell_pec_scatter.cpp

What the code does#

  1. Loads a surface mesh.

  2. Builds RWG basis functions on the mesh.

  3. Assembles the EFIE operator and excitation.

  4. Solves the resulting linear system.

  5. Writes outputs for visualization / post-processing.

How to run#

./build/maxwell_pec_scatter ../sphere_0.1.off 6.283185307179586

Full source#

  1#include "biekit/biekit.h"
  2
  3#include <cmath>
  4#include <complex>
  5#include <iostream>
  6#include <fstream>
  7#include <stdexcept>
  8#include <string>
  9
 10class MaxwellPecScatterProblem {
 11public:
 12    void main(int argc, char** argv) {
 13        std::string mesh_path = "../sphere_0.1.off";
 14
 15        // time convention: e^{-j w t}
 16        const double omega = 2.0 * kPi * 300e6;
 17        const double mu = 4.0 * kPi * 1e-7;
 18        const double eps = 8.854187817e-12;
 19        const std::complex<double> wave_number = omega * std::sqrt(mu * eps);
 20
 21        std::cout << "Loading mesh: " << mesh_path << "\n";
 22        Mesh mesh = load_mesh(mesh_path);
 23
 24        auto grid = build_grid<2, 3>(mesh);
 25
 26        RwgTriFiniteElement<int> fe_tri;
 27        FunctionSpace<SpaceType::Hdiv, 2, 3> trial_space(grid, fe_tri);
 28        FunctionSpace<SpaceType::Hdiv, 2, 3> test_space(grid, fe_tri);
 29
 30        const std::array<double, 3> khat{0.0, 0.0, 1.0};
 31        const std::array<std::complex<double>, 3> E0{
 32            std::complex<double>{1.0, 0.0},
 33            std::complex<double>{0.0, 0.0},
 34            std::complex<double>{0.0, 0.0}
 35        };
 36
 37        std::cout << "Assembling RHS...\n";
 38        const int quad_order = 2;
 39        std::vector<std::complex<double>> rhs = assemble_rhs_incident_E(test_space, wave_number, quad_order, khat, E0);
 40
 41        std::vector<std::complex<double>> sol;
 42        sol.resize(static_cast<std::size_t>(trial_space.dofmap().global_dof_size()));
 43        blas1::fill(std::span<std::complex<double>>(sol), std::complex<double>{0.0, 0.0});
 44
 45        std::cout << "Binding operator...\n";
 46        ElectricFieldOperator<2, 3> A(trial_space, test_space, wave_number, quad_order);
 47
 48        std::cout << "Solving with GMRES...\n";
 49        GmresComplex solver;
 50        const auto info = solver.solve(A, sol, rhs);
 51
 52        std::cout << "Solution norm_inf = " << blas1::norm_inf(sol) << "\n";
 53
 54        std::vector<double> surface_current_mag;
 55        std::vector<double> surface_charge_mag;
 56        MaxwellComputeCellCenterCurrentAndCharge(trial_space, sol, omega, surface_current_mag, surface_charge_mag);
 57
 58        std::cout << "Computing RCS on xoz-plane...\n";
 59        {
 60            std::ofstream ofs("rcs_xoz.csv");
 61            if (!ofs) {
 62                throw std::runtime_error("Failed to open rcs_xoz.csv");
 63            }
 64            ofs.setf(std::ios::scientific);
 65            ofs.precision(16);
 66            ofs << "theta_deg,rcs_db\n";
 67
 68            const int n = 181; // 0..180 deg inclusive
 69            for (int i = 0; i < n; ++i) {
 70                const double theta = (kPi * i) / 180.0;
 71                const double rcs_db = MaxwellFarFieldRcs(
 72                    trial_space,
 73                    sol,
 74                    wave_number,
 75                    std::sqrt(mu / eps),
 76                    theta,
 77                    0.0,
 78                    true);
 79                ofs << (theta * 180.0 / kPi) << "," << rcs_db << "\n";
 80            }
 81
 82            if (!ofs) {
 83                throw std::runtime_error("Failed while writing rcs_xoz.csv");
 84            }
 85        }
 86
 87        VtkWriteOptions vtk_opt;
 88        vtk_opt.title = "maxwell_pec_scatter";
 89        WriteVtkPolyDataCellFields(
 90            mesh,
 91            "maxwell_solution.vtk",
 92            {"surface_current", "surface_charge"},
 93            {surface_current_mag, surface_charge_mag},
 94            vtk_opt);
 95        std::cout << "[success] write vtk: maxwell_solution.vtk\n";
 96        std::cout << "[success] write csv: rcs_xoz.csv\n";
 97        exit_code_ = info.converged ? 0 : 2;
 98    }
 99
100    int exit_code() const { return exit_code_; }
101
102private:
103    int exit_code_ = 0;
104};
105
106int main(int argc, char** argv) {
107    biekit::enable_floating_point_exceptions();
108
109    try {
110        MaxwellPecScatterProblem problem;
111        problem.main(argc, argv);
112        return problem.exit_code();
113    } catch (const std::exception& e) {
114        std::cerr << "Fatal error: " << e.what() << "\n";
115        return 1;
116    }
117}