Implementation of curved elements for Hesthaven Evolution Operator

src/components/DGOperatorFactory.cpp

@@ -63,7 +63,7 @@ namespace maxwell {
 		const BilinearForm& op2,
 		FiniteElementSpace& fes)
 	{
-		auto aux = mfem::Mult(op1.SpMat(), op2.SpMat());
+		auto aux = std::unique_ptr<SparseMatrix>(mfem::Mult(op1.SpMat(), op2.SpMat()));
 		auto res = std::make_unique<BilinearForm>(&fes);
 		res->Assemble();
 		res->Finalize();
@@ -397,15 +397,16 @@ namespace maxwell {
 		GlobalIndices globalId(fes_.GetNDofs());
 		for (auto f : { E, H }) {
 			auto MInv = buildGlobalInverseMassMatrixOperator();
-			auto intBdrFluxZeroNormal = buildByMult(*MInv[f], *buildZeroNormalIBFISubOperator(f), fes_);
+			auto dense = buildByMult(*MInv[f], *buildZeroNormalIBFISubOperator(f), fes_)->SpMat().ToDenseMatrix();
 			for (auto d : { X, Y, Z }) {
 				loadBlockInGlobalAtIndices(
-					*intBdrFluxZeroNormal->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][d], globalId.index[f][d]),
 					-1.0
 				);
 			}
+			delete dense;
 		}
 	}
 	void DGOperatorFactory::addGlobalOneNormalIBFIOperators(SparseMatrix* global)
@@ -416,18 +417,20 @@ namespace maxwell {
 			for (auto x{ X }; x <= Z; x++) {
 				auto y = (x + 1) % 3;
 				auto z = (x + 2) % 3;
+				auto dense = buildByMult(*MInv[f], *buildOneNormalIBFISubOperator(altField(f), { x }), fes_)->SpMat().ToDenseMatrix();
 				loadBlockInGlobalAtIndices(
-					*buildByMult(*MInv[f], *buildOneNormalIBFISubOperator(altField(f), { x }), fes_)->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][y], globalId.index[altField(f)][z]),
-					-1.0 + double(f) * 2.0
+					1.0 - double(f) * 2.0
 				);
 				loadBlockInGlobalAtIndices(
-					*buildByMult(*MInv[f], *buildOneNormalIBFISubOperator(altField(f), { x }), fes_)->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][z], globalId.index[altField(f)][y]),
-					1.0 - double(f) * 2.0
+					-1.0 + double(f) * 2.0
 				);
+				delete dense;
 			}
 		}
 	}
@@ -438,11 +441,13 @@ namespace maxwell {
 			auto MInv = buildGlobalInverseMassMatrixOperator();
 			for (auto d{ X }; d <= Z; d++) {
 				for (auto d2{ X }; d2 <= Z; d2++) {
+					auto dense = buildByMult(*MInv[f], *buildTwoNormalIBFISubOperator(f, { d, d2 }), fes_)->SpMat().ToDenseMatrix();
 					loadBlockInGlobalAtIndices(
-						*buildByMult(*MInv[f], *buildTwoNormalIBFISubOperator(f, { d, d2 }), fes_)->SpMat().ToDenseMatrix(),
+						*dense,
 						*global,
 						std::make_pair(globalId.index[f][d], globalId.index[f][d2])
 					);
+					delete dense;
 				}
 			}
 		}
@@ -455,19 +460,20 @@ namespace maxwell {
 			for (auto x{ X }; x <= Z; x++) {
 				auto y = (x + 1) % 3;
 				auto z = (x + 2) % 3;
-				auto dir = buildByMult(*MInv[f], *buildDerivativeSubOperator(x), fes_);
+				auto dense = buildByMult(*MInv[f], *buildDerivativeSubOperator(x), fes_)->SpMat().ToDenseMatrix();
 				loadBlockInGlobalAtIndices(
-					*dir->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][z], globalId.index[altField(f)][y]),
 					1.0 - double(f) * 2.0
 				);
 				loadBlockInGlobalAtIndices(
-					*dir->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][y], globalId.index[altField(f)][z]),
 					-1.0 + double(f) * 2.0
 				);
+				delete dense;
 			}
 		}
 	}
@@ -476,15 +482,16 @@ namespace maxwell {
 		GlobalIndices globalId(fes_.GetNDofs());
 		for (auto f : { E, H }) {
 			auto MInv = buildGlobalInverseMassMatrixOperator();
-			auto fluxZeroNormal = buildByMult(*MInv[f], *buildZeroNormalSubOperator(f), fes_);
+			auto dense = buildByMult(*MInv[f], *buildZeroNormalSubOperator(f), fes_)->SpMat().ToDenseMatrix();
 			for (auto d : { X, Y, Z }) {
 				loadBlockInGlobalAtIndices(
-					*fluxZeroNormal->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][d], globalId.index[f][d]),
 					-1.0
 				);
 			}
+			delete dense;
 		}
 	}
 	void DGOperatorFactory::addGlobalOneNormalOperators(SparseMatrix* global)
@@ -495,21 +502,20 @@ namespace maxwell {
 			for (auto x{ X }; x <= Z; x++) {
 				auto y = (x + 1) % 3;
 				auto z = (x + 2) % 3;
-
-				auto oneNormal = buildByMult(*MInv[f], *buildOneNormalSubOperator(altField(f), { x }), fes_);
-
+				auto dense = buildByMult(*MInv[f], *buildOneNormalSubOperator(altField(f), { x }), fes_)->SpMat().ToDenseMatrix();
 				loadBlockInGlobalAtIndices(
-					*oneNormal->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][y], globalId.index[altField(f)][z]),
 					1.0 - double(f) * 2.0
 				);
 				loadBlockInGlobalAtIndices(
-					*oneNormal->SpMat().ToDenseMatrix(),
+					*dense,
 					*global,
 					std::make_pair(globalId.index[f][z], globalId.index[altField(f)][y]),
 					-1.0 + double(f) * 2.0
 				);
+				delete dense;
 			}
 		}
 	}
@@ -520,11 +526,13 @@ namespace maxwell {
 			auto MInv = buildGlobalInverseMassMatrixOperator();
 			for (auto d{ X }; d <= Z; d++) {
 				for (auto d2{ X }; d2 <= Z; d2++) {
+					auto dense = buildByMult(*MInv[f], *buildTwoNormalSubOperator(f, { d, d2 }), fes_)->SpMat().ToDenseMatrix();
 					loadBlockInGlobalAtIndices(
-						*buildByMult(*MInv[f], *buildTwoNormalSubOperator(f, { d, d2 }), fes_)->SpMat().ToDenseMatrix(),
+						*dense,
 						*global,
 						std::make_pair(globalId.index[f][d], globalId.index[f][d2])
 					);
+					delete dense;
 				}
 			}
 		}
@@ -533,15 +541,15 @@ namespace maxwell {
 	std::unique_ptr<SparseMatrix> DGOperatorFactory::buildTFSFGlobalOperator()
 	{
 
-		std::unique_ptr<SparseMatrix> res;
+		std::unique_ptr<SparseMatrix> res = std::make_unique<SparseMatrix>(6 * fes_.GetNDofs(), 6 * fes_.GetNDofs());
 
 #ifdef SHOW_TIMER_INFORMATION
 		std::cout << "Elapsed time (ms): " + std::to_string(std::chrono::duration_cast<std::chrono::milliseconds>
 			(std::chrono::high_resolution_clock::now() - startTime).count()) << std::endl;
 		std::cout << "Assembling TFSF Inverse Mass One-Normal Operators" << std::endl;
 #endif
 
-		addGlobalOneNormalIBFIOperators(res.get());
+		addGlobalOneNormalOperators(res.get());
 
 		if (pd_.opts.fluxType == FluxType::Upwind) {
 
@@ -551,7 +559,7 @@ namespace maxwell {
 			std::cout << "Assembling TFSF Inverse Mass Zero-Normal Operators" << std::endl;
 #endif	
 
-			addGlobalZeroNormalIBFIOperators(res.get());
+			addGlobalZeroNormalOperators(res.get());
 
 
 #ifdef SHOW_TIMER_INFORMATION	
@@ -560,7 +568,7 @@ namespace maxwell {
 			std::cout << "Assembling TFSF Inverse Mass Two-Normal Operators" << std::endl;
 #endif	
 
-			addGlobalTwoNormalIBFIOperators(res.get());
+			addGlobalTwoNormalOperators(res.get());
 
 		}
 
@@ -571,7 +579,7 @@ namespace maxwell {
 	std::unique_ptr<SparseMatrix> DGOperatorFactory::buildGlobalOperator()
 	{
 
-		std::unique_ptr<SparseMatrix> res;
+		std::unique_ptr<SparseMatrix> res = std::make_unique<SparseMatrix>(6 * fes_.GetNDofs(), 6 * fes_.GetNDofs());
 
 		if (pd_.model.getInteriorBoundaryToMarker().size() != 0) { //IntBdrConds
 

src/components/DGOperatorFactory.h

@@ -1,7 +1,7 @@
 #pragma once
 
 #include "ProblemDefinition.h"
-#include "evolution/GlobalMethods.h"
+#include "evolution/GlobalEvolution.h"
 #include "evolution/HesthavenEvolutionMethods.h"
 #include "evolution/MaxwellEvolutionMethods.h"
 

src/components/RCSManager.cpp

@@ -69,7 +69,7 @@ double func_exp_real_part_2D(const Vector& x, const double freq, const Phi phi)
 }
 double func_exp_imag_part_2D(const Vector& x, const double freq, const Phi phi)
 {
-	auto landa = physicalConstants::speedOfLight_SI/  freq;
+	auto landa = physicalConstants::speedOfLight_SI / freq;
 	auto wavenumber = 2.0 * M_PI / landa;
 	auto rad_term = wavenumber * (x[0] * cos(phi) + x[1] * sin(phi));
 	return sin(rad_term);
@@ -162,15 +162,15 @@ const double getTime(const std::string& timePath)
 	return std::stod(timeString);
 }
 
-std::map<SphericalAngles, Freq2Value> fillPostDataMaps(const std::vector<double>& frequencies, const std::vector<SphericalAngles>& angleVec)
+std::map<SphericalAngles, Freq2Value> initAngles2FreqValues(const std::vector<double>& frequencies, const std::vector<SphericalAngles>& angleVec)
 {
 	std::map<SphericalAngles, Freq2Value> res;
 	for (const auto& angpair : angleVec) {
-		Freq2Value inner;
+		Freq2Value f2v;
 		for (const auto& f : frequencies) {
-			inner.emplace(f, 0.0);
+			f2v.emplace(f, 0.0);
 		}
-		res.emplace(angpair, inner);
+		res.emplace(angpair, f2v);
 	}
 	return res;
 }
@@ -186,7 +186,7 @@ PlaneWaveData buildPlaneWaveData(const json& json)
 		}
 	}
 
-	if (mean == -1e5 || delay == -1e5) {
+	if (std::abs(mean - 1e5) > 1e-6 || std::abs(delay - 1e5) > 1e-6) {
 		throw std::runtime_error("Verify PlaneWaveData inputs for RCS normalization term.");
 	}
 
@@ -263,9 +263,9 @@ std::vector<double> buildNormalizationTerm(const std::string& json_path, const s
 		res[f] = physicalConstants::vacuumPermittivity_SI * std::pow(std::abs(freq_val), 2.0);
 	}
 	auto max = *std::max_element(std::begin(res), std::end(res));
-	//for (auto f{ 0 }; f < res.size(); f++) {
-	//	res[f] /= max;
-	//}
+	for (auto f{ 0 }; f < res.size(); f++) {
+		res[f] /= max;
+	}
 
 	trimLowMagFreqs(freq2complex, frequencies);
 	exportIncidentGaussData(time, gauss_val, json_path);
@@ -384,16 +384,13 @@ std::pair<std::complex<double>, std::complex<double>> RCSManager::performRCSCalc
 	return std::pair<std::complex<double>, std::complex<double>>(phi_value, theta_value);
 }
 
-RCSManager::RCSManager(const std::string& path, const std::string& json_path, double dt, int steps, const std::vector<SphericalAngles>& angle_vec)
+RCSManager::RCSManager(const std::string& path, const std::string& json_path, std::vector<double>& frequencies, const std::vector<SphericalAngles>& angle_vec)
 {
 
 	Mesh mesh{ Mesh::LoadFromFile(path + "/mesh", 1, 0) };
 	getFESFromGF(mesh, path);
 
-	const double f_max = 2.0 / dt;
-	std::vector<double> frequencies;
-	frequencies = logspace(6.0, 7.7, 100);
-	auto RCSdata{ fillPostDataMaps(frequencies, angle_vec) };
+	auto RCSdata{ initAngles2FreqValues(frequencies, angle_vec) };
 
 	auto normalization_term{ buildNormalizationTerm(json_path, path, frequencies) };
 

src/components/RCSManager.h

@@ -64,7 +64,7 @@ struct RCSData {
 class RCSManager {
 public:
 
-	RCSManager(const std::string& data_path, const std::string& json_path, double f_max, int steps, const std::vector<SphericalAngles>& angle);
+	RCSManager(const std::string& data_path, const std::string& json_path, std::vector<double>& frequencies, const std::vector<SphericalAngles>& angle);
 
 private:
 

src/components/Sources.cpp

@@ -47,7 +47,6 @@ double InitialField::eval(
 	const Position& p, const Time& t,
 	const FieldType& f, const Direction& d) const
 {
-	assert(p.Size() == center_.Size());
 
 	if (f != fieldType_) {
 		return 0.0;

src/components/Sources.h

@@ -41,6 +41,7 @@ class InitialField : public Source {
 	FieldType& fieldType() { return fieldType_; }
 	Polarization& polarization() { return polarization_; }
 	Function* magnitude() { return magnitude_.get(); }
+	Position& center() { return center_; }
 
 private:
 	std::unique_ptr<Function> magnitude_;

src/driver/driver.cpp

@@ -98,9 +98,18 @@ std::unique_ptr<InitialField> buildResonantModeInitialField(
 	const std::vector<std::size_t>& modes = { 1 })
 {
 	Sources res;
-	Source::Position center_((int)modes.size());
-	center_ = 0.0;
-	return std::make_unique<InitialField>(SinusoidalMode{ modes }, ft, p, center_);
+	Source::Position center((int)modes.size());
+	center = 0.0;
+	return std::make_unique<InitialField>(SinusoidalMode{ modes }, ft, p, center);
+}
+
+std::unique_ptr<InitialField> buildBesselJ6InitialField(
+	const FieldType& ft = E,
+	const Source::Polarization& p = Source::Polarization({ 0.0, 0.0, 1.0 }))
+{
+	Sources res;
+	Source::Position center = Source::Position({ 0.0, 0.0, 0.0 });
+	return std::make_unique<InitialField>(BesselJ6(), ft, p, center);
 }
 
 std::unique_ptr<Planewave> buildGaussianPlanewave(
@@ -136,6 +145,12 @@ Sources buildSources(const json& case_data)
 					case_data["sources"][s]["magnitude"]["modes"])
 				);
 			}
+			else if (case_data["sources"][s]["magnitude"]["type"] == "besselj6") {
+				res.add(buildBesselJ6InitialField(
+					assignFieldType(case_data["sources"][s]["field_type"]),
+					assemble3DVector(case_data["sources"][s]["polarization"]))
+				);
+			}
 		}
 		else if (case_data["sources"][s]["type"] == "totalField") {
 			res.add(buildGaussianPlanewave(
@@ -193,6 +208,10 @@ SolverOptions buildSolverOptions(const json& case_data)
 		if (case_data["solver_options"].contains("hesthaven_operator")) {
 			res.setHesthavenOperator(case_data["solver_options"]["hesthaven_operator"]);
 		}
+
+		if (case_data["solver_options"].contains("global_operator")) {
+			res.setGlobalOperator(case_data["solver_options"]["global_operator"]);
+		}
 	}
 	return res;
 }

src/evolution/CMakeLists.txt

@@ -8,7 +8,7 @@ add_library(maxwell-evolution STATIC
 	"Fields.cpp" 
 	"HesthavenEvolutionMethods.cpp"
 	"HesthavenEvolution.cpp" 
-	"GlobalMethods.cpp")
+	"GlobalEvolution.cpp")
 
 get_filename_component(PARENT_DIR ../ ABSOLUTE)
 include_directories(${PARENT_DIR})