Discontinuous Galerkin based Radar Cross Section Post-Processing works

src/components/FarField.cpp

@@ -97,7 +97,7 @@ std::unique_ptr<FunctionCoefficient> buildFC_3D(const Frequency freq, const Sphe
 		break;
 	}
 	FunctionCoefficient res(f);
-	return std::make_unique<FunctionCoefficient>(res);
+	return std::make_unique<FunctionCoefficient>(res); 
 }
 
 std::complex<double> complexInnerProduct(ComplexVector& first, ComplexVector& second)
@@ -107,7 +107,7 @@ std::complex<double> complexInnerProduct(ComplexVector& first, ComplexVector& se
 	}
 	std::complex<double> res(0.0, 0.0);
 	for (int i{ 0 }; i < first.size(); i++) {
-		res += first[i] * std::conj(second[i]); //<x,y> = sum(x_i * conj(y_i))
+		res += first[i] * second[i];
 	}
 	return res;
 }
@@ -175,20 +175,24 @@ const Time getTime(const std::string& timePath)
 
 PlaneWaveData buildPlaneWaveData(const json& json)
 {
-	double mean(-1e5), delay(-1e5);
+	double spread(-1e5);
+	mfem::Vector mean;
+	double projMean(-1e5);
 
 	for (auto s{ 0 }; s < json["sources"].size(); s++) {
 		if (json["sources"][s]["type"] == "planewave") {
-			mean = json["sources"][s]["magnitude"]["spread"];
-			delay = json["sources"][s]["magnitude"]["delay"];
+			spread = json["sources"][s]["magnitude"]["spread"];
+			mean = driver::assemble3DVector(json["sources"][s]["magnitude"]["mean"]);
+		    mfem::Vector propagation = driver::assemble3DVector(json["sources"][s]["propagation"]);
+			projMean = mean * propagation / propagation.Norml2();
 		}
 	}
 
-	if (std::abs(mean - 1e5) < 1e-6 || std::abs(delay - 1e5) < 1e-6) {
+	if (std::abs(spread - 1e5) < 1e-6 || std::abs(projMean - 1e5) < 1e-6) {
 		throw std::runtime_error("Verify PlaneWaveData inputs for RCS normalization term.");
 	}
 
-	return PlaneWaveData(mean / physicalConstants::speedOfLight, delay / physicalConstants::speedOfLight);
+	return PlaneWaveData(spread, projMean);
 }
 
 std::vector<double> buildTimeVector(const std::string& data_path)
@@ -205,11 +209,11 @@ std::vector<double> buildTimeVector(const std::string& data_path)
 	return res;
 }
 
-std::vector<double> evaluateGaussianVector(std::vector<Time>& time, double delay, double mean)
+std::vector<double> evaluateGaussianVector(std::vector<Time>& time, double spread, double mean)
 {
 	std::vector<double> res(time.size());
 	for (int t = 0; t < time.size(); ++t) {
-		res[t] = exp(-0.5 * std::pow((time[t] - delay) / mean, 2.0));
+		res[t] = exp(-0.5 * std::pow((time[t] - std::abs(mean)) / spread, 2.0));
 	}
 	return res;
 }
@@ -252,8 +256,8 @@ FreqFields calculateFreqFields(Mesh& mesh, const std::vector<Frequency>& frequen
 		std::vector<GridFunction> A;
 		for (auto const& dir_entry : std::filesystem::directory_iterator(path)) {
 			if (dir_entry.path().generic_string().substr(dir_entry.path().generic_string().size() - 4) != "mesh" &&
-				dir_entry.path().generic_string().substr(dir_entry.path().generic_string().size() - 3) != "rcs" &&
-				dir_entry.path().generic_string().substr(dir_entry.path().generic_string().size() - 8) != "farfield") {
+				dir_entry.path().generic_string().substr(dir_entry.path().generic_string().size() - 4) != "/rcs" &&
+				dir_entry.path().generic_string().substr(dir_entry.path().generic_string().size() - 9) != "/farfield") {
 				A.push_back(getGridFunction(mesh, dir_entry.path().generic_string() + field));
 			}
 		}

src/components/FarField.h

@@ -39,12 +39,12 @@ struct SphericalAngles {
 };
 
 struct PlaneWaveData {
+	double spread;
 	double mean;
-	double delay;
 
-	PlaneWaveData(double m, double dl) :
-		mean(m),
-		delay(dl) {};
+	PlaneWaveData(double s, double m) :
+		spread(s),
+		mean(m) {};
 };
 
 struct FreqFields {
@@ -78,6 +78,8 @@ double func_exp_imag_part_3D(const Position&, const Frequency, const SphericalAn
 std::unique_ptr<FunctionCoefficient> buildFC_2D(const Frequency, const SphericalAngles&, bool isReal);
 std::unique_ptr<FunctionCoefficient> buildFC_3D(const Frequency, const SphericalAngles&, bool isReal);
 
+double calcPsiAngle3D(const Vector& p, const SphericalAngles& angles);
+
 std::complex<double> complexInnerProduct(ComplexVector& first, ComplexVector& second);
 
 std::unique_ptr<FiniteElementSpace> buildFESFromGF(Mesh&, const std::string& data_path);

src/components/RCSDataExtractor.cpp

@@ -216,8 +216,11 @@ RCSDataExtractor::RCSDataExtractor(const std::string data_folder, const std::str
 	}
 
 	// This takes care of assembling the incoming field values for the incident planewave.
-	auto planewave_data{ buildPlaneWaveData(case_data) };
-	Gaussian gauss{ planewave_data.mean, mfem::Vector({-planewave_data.delay}) };
+	auto planewave_spread = case_data["sources"][0]["magnitude"]["spread"];
+	auto planewave_mean = driver::assemble3DVector(case_data["sources"][0]["magnitude"]["mean"]);
+	auto planewave_propagation = driver::assemble3DVector(case_data["sources"][0]["propagation"]);
+	auto proj_mean = planewave_mean * planewave_propagation / planewave_propagation.Norml2();
+	Gaussian gauss{ planewave_spread, mfem::Vector({proj_mean}) };
 	auto sources = driver::buildSources(case_data);
 	std::vector<double> ExInc(time.size()), EyInc(time.size()), EzInc(time.size());
 	for (const auto& source : sources) {

src/components/RCSManager.cpp

@@ -9,9 +9,13 @@ using json = nlohmann::json;
 static std::vector<double> buildNormalizationTerm(const std::string& json_path, const std::string& path, std::vector<double>& frequencies)
 {
 
-	auto planewave_data{ buildPlaneWaveData(driver::parseJSONfile(json_path)) };
+	auto case_data = driver::parseJSONfile(json_path);
+	double spread = case_data["sources"][0]["magnitude"]["spread"];
+	mfem::Vector mean = driver::assemble3DVector(case_data["sources"][0]["magnitude"]["mean"]);
+	mfem::Vector propagation = driver::assemble3DVector(case_data["sources"][0]["propagation"]);
+	double projMean = mean * propagation / propagation.Norml2();;
 	std::vector<double> time{ buildTimeVector(path) };
-	std::vector<double> gauss_val{ evaluateGaussianVector(time, planewave_data.delay, planewave_data.mean) };
+	std::vector<double> gauss_val{ evaluateGaussianVector(time, spread, projMean) };
 
 	std::map<double, std::complex<double>> freq2complex;
 	std::vector<double> res(frequencies.size(), 0.0);

src/driver/driver.cpp

@@ -114,24 +114,29 @@ std::unique_ptr<InitialField> buildBesselJ6InitialField(
 
 std::unique_ptr<TotalField> buildGaussianPlanewave(
 	double spread,
-	double delay,
+	const Source::Position mean,
 	const Source::Polarization& pol,
 	const Source::Propagation& dir,
 	const FieldType ft = FieldType::E
 )
 {
-	Gaussian gauss{ spread, mfem::Vector({-delay}) };
+	Position projMean(3);
+	projMean = 0.0;
+	for (auto v = 0; v < mean.Size(); v++) {
+		projMean[v] = mean[v];
+	}
+	Gaussian gauss{ spread, mfem::Vector({projMean * dir / dir.Norml2()})};
 	Planewave pw(gauss, pol, dir, ft);
 	return std::make_unique<TotalField>(pw);
 }
 
 std::unique_ptr<TotalField> buildDerivGaussDipole(
 	const double length, 
 	const double gaussianSpread, 
-	const double gaussDelay
+	const double gaussMean
 ) 
 {
-	DerivGaussDipole dip(length, gaussianSpread, gaussDelay);
+	DerivGaussDipole dip(length, gaussianSpread, gaussMean);
 	return std::make_unique<TotalField>(dip);
 }
 
@@ -166,7 +171,7 @@ Sources buildSources(const json& case_data)
 		else if (case_data["sources"][s]["type"] == "planewave") {
 			res.add(buildGaussianPlanewave(
 				case_data["sources"][s]["magnitude"]["spread"],
-				case_data["sources"][s]["magnitude"]["delay"],
+				assemble3DVector(case_data["sources"][s]["magnitude"]["mean"]),
 				assemble3DVector(case_data["sources"][s]["polarization"]),
 				assemble3DVector(case_data["sources"][s]["propagation"]),
 				FieldType::E)
@@ -176,7 +181,7 @@ Sources buildSources(const json& case_data)
 			res.add(buildDerivGaussDipole(
 				case_data["sources"][s]["magnitude"]["length"],
 				case_data["sources"][s]["magnitude"]["spread"],
-				case_data["sources"][s]["magnitude"]["delay"])
+				case_data["sources"][s]["magnitude"]["mean"])
 			);
 		}
 		else {

src/driver/driver.h

@@ -11,6 +11,8 @@ using json = nlohmann::json;
 namespace maxwell::driver {
 	json parseJSONfile(const std::string& case_name);
 
+	mfem::Vector assemble3DVector(const json& input);
+
 	maxwell::Solver buildSolverJson(const std::string& case_name, const bool isTest = true);
 	maxwell::Solver buildSolver(const json& case_data, const std::string& case_path, const bool isTest);
 

src/evolution/Fields.cpp

@@ -4,8 +4,9 @@ namespace maxwell {
 
 using namespace mfem;
 
-Fields::Fields(mfem::FiniteElementSpace& fes)
+Fields::Fields(FiniteElementSpace& fes)
 {
+    global_fes_ = std::make_unique<FiniteElementSpace>(fes.GetMesh(), dynamic_cast<const DG_FECollection*>(fes.FEColl()), 3);
     allDOFs_.SetSize(6 * fes.GetNDofs());
     allDOFs_ = 0.0;
     for (int d = X; d <= Z; d++) {
@@ -14,6 +15,17 @@ Fields::Fields(mfem::FiniteElementSpace& fes)
         e_[d].SetSpace(&fes);
         h_[d].SetSpace(&fes);
     }
+
+    e_global_.SetSpace(global_fes_.get());
+    h_global_.SetSpace(global_fes_.get());
+
+    e_global_.SetVector(e_[X], 0);
+    e_global_.SetVector(e_[Y], fes.GetNDofs());
+    e_global_.SetVector(e_[Z], 2 * fes.GetNDofs());
+    h_global_.SetVector(h_[X], 0);
+    h_global_.SetVector(h_[Y], fes.GetNDofs());
+    h_global_.SetVector(h_[Z], 2 * fes.GetNDofs());
+
 }
 
 GridFunction& Fields::get(const FieldType& f, const Direction& d)
@@ -40,5 +52,28 @@ const GridFunction& Fields::get(const FieldType& f, const Direction& d) const
     }
 }
 
+GridFunction& Fields::get(const FieldType& f)
+{
+    assert(f == E || f == H);
+    if (f == E) {
+        return e_global_;
+    }
+    else {
+        return h_global_;
+    }
+}
+
+void Fields::updateGlobal() 
+{
+    auto offset = global_fes_->GetNDofs();
+    e_global_.SetVector(e_[X], 0);
+    e_global_.SetVector(e_[Y], offset);
+    e_global_.SetVector(e_[Z], 2 * offset);
+    h_global_.SetVector(h_[X], 0);
+    h_global_.SetVector(h_[Y], offset);
+    h_global_.SetVector(h_[Z], 2 * offset);
+}
+
+
 
 }

src/evolution/Fields.h

@@ -12,14 +12,20 @@ class Fields {
 
     mfem::GridFunction& get(const FieldType&, const Direction&);
     const mfem::GridFunction& get(const FieldType&, const Direction&) const;
+
+    mfem::GridFunction& get(const FieldType&);
 
     mfem::Vector& allDOFs() { return allDOFs_; }
     const mfem::Vector& allDOFs() const { return allDOFs_; }
 
     double getNorml2() const { return allDOFs_.Norml2(); }
 
+    void updateGlobal();
+
 private:
     mfem::Vector allDOFs_;
+    std::unique_ptr<mfem::FiniteElementSpace> global_fes_;
     std::array<mfem::GridFunction, 3> e_, h_;
+    mfem::GridFunction e_global_, h_global_;
 };
 }