Merge pull request OpenSEMBA#6 from lmdiazangulo/main

Compilation issue. Probe and Material rework.
EMinsight · May 18, 2022 · 38d9b9f · 38d9b9f
2 parents 3d085fa + d36d5f8
commit 38d9b9f
Show file tree

Hide file tree

Showing 13 changed files with 61 additions and 80 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -14,10 +14,6 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 
 project(maxwell)
 
-
-find_package(mfem CONFIG REQUIRED)
-include_directories(${MFEM_INCLUDE_DIRS})
-
 add_subdirectory(src)
 
 include_directories(${CMAKE_CURRENT_LIST_DIR}/src)	  

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -1,4 +1,3 @@
 cmake_minimum_required(VERSION 3.0)
 
-add_subdirectory(maxwell)
-
+add_subdirectory(maxwell)
diff --git a/src/maxwell/CMakeLists.txt b/src/maxwell/CMakeLists.txt
@@ -1,12 +1,18 @@
 cmake_minimum_required(VERSION 3.0)
 
+find_package(mfem CONFIG REQUIRED)
+include_directories(${MFEM_INCLUDE_DIRS})
+
 add_library(maxwell_solvers STATIC 
-	"Solver.cpp" "Solver.h" 
-	"FiniteElementEvolution.cpp" "FiniteElementEvolution.h" 
-	"BilinearIntegrators.cpp" "BilinearIntegrators.h"    
-	"Types.h" "Options.h"
+	"Solver.h" "Solver.cpp" 
+	"FiniteElementEvolution.h" "FiniteElementEvolution.cpp" 
+	"BilinearIntegrators.h" "BilinearIntegrators.cpp"
+	"Types.h" 
+	"Options.h"
 	"Material.h" "Material.cpp"
 	"Model.h" "Model.cpp" 
 	"Probes.h" "Probes.cpp" 
-	"Sources.h" "Sources.cpp"  )
+	"Sources.h" "Sources.cpp")
+
+target_link_libraries(maxwell_solvers mfem)
 
diff --git a/src/maxwell/FiniteElementEvolution.cpp b/src/maxwell/FiniteElementEvolution.cpp
@@ -8,9 +8,6 @@ FiniteElementEvolutionNoCond::FiniteElementEvolutionNoCond(FiniteElementSpace* f
 	fes_(fes),
 	model_(model)
 {
-	initializeMaterialParameterVectors();
-	getMaterialParameterVectors();
-
 	for (int fInt = FieldType::E; fInt <= FieldType::H; fInt++) {
 		FieldType f = static_cast<FieldType>(fInt);
 		for (int fInt2 = FieldType::E; fInt2 <= FieldType::H; fInt2++) {
@@ -25,44 +22,25 @@ FiniteElementEvolutionNoCond::FiniteElementEvolutionNoCond(FiniteElementSpace* f
 	}
 }
 
-void FiniteElementEvolutionNoCond::initializeMaterialParameterVectors()
-{
-	std::vector<std::pair<attribute, Material>> matMap = model_.getAttToMatVec();
-	int max = 1;
-	for (int i = 0; i < matMap.size(); i++) {
-		if (matMap[i].first > max) {
-			max = matMap[i].first;
-		}
-	}
-	eps_.SetSize(max);
-	mu_.SetSize(max);
-	eps_ = 1.0; mu_ = 1.0;
-}
-
-void FiniteElementEvolutionNoCond::getMaterialParameterVectors()
-{
-	std::vector<std::pair<attribute, Material>> matVec = model_.getAttToMatVec();
-	for(const auto & it : matVec) {
-		eps_[it.first-1] = it.second.getPermittivity();
-		mu_[it.first-1] = it.second.getPermeability();
-	}
-}
-
 Vector
 FiniteElementEvolutionNoCond::buildPieceWiseArgVector(const FieldType& f) const
 {
 	Vector res;
+	res.SetSize(model_.getAttToMat().size());
 
-	switch (f) {
-	case FieldType::E:
-		res.SetSize(eps_.Size());
-		res = eps_;
-		break;
-	case FieldType::H:
-		res.SetSize(mu_.Size());
-		res = mu_;
-		break;
+	std::size_t i = 0;
+	for (auto const& kv : model_.getAttToMat()) {
+		switch (f) {
+		case FieldType::E:
+			res[i] = kv.second.getPermittivity();
+			break;
+		case FieldType::H:
+			res[i] = kv.second.getPermeability();
+			break;
+		}
+		i++;
 	}
+
 	return res;
 }
 

diff --git a/src/maxwell/FiniteElementEvolution.h b/src/maxwell/FiniteElementEvolution.h
@@ -10,8 +10,6 @@
 
 namespace maxwell {
 
-using namespace mfem;
-
 class FiniteElementEvolutionNoCond : public TimeDependentOperator {
 public:
 
@@ -37,7 +35,6 @@ class FiniteElementEvolutionNoCond : public TimeDependentOperator {
 
 	FiniteElementSpace* fes_;
 	Options opts_;
-	Vector eps_, mu_;
 	Model model_;
 
 	std::array<std::array<Operator, 3>, 2> MS_;
@@ -54,9 +51,6 @@ class FiniteElementEvolutionNoCond : public TimeDependentOperator {
 
 	Operator buildByMult(const BilinearForm*, const BilinearForm*) const;
 
-	void getMaterialParameterVectors();
-	void initializeMaterialParameterVectors();
-
 	FluxCoefficient interiorFluxCoefficient() const;
 	FluxCoefficient interiorPenaltyFluxCoefficient() const;
 	FluxCoefficient boundaryFluxCoefficient(const FieldType&) const;

diff --git a/src/maxwell/Model.cpp b/src/maxwell/Model.cpp
@@ -2,7 +2,7 @@
 
 namespace maxwell {
 
-Model::Model(Mesh& mesh, std::vector<std::pair<attribute, Material>>& matVec) :
+Model::Model(Mesh& mesh, const AttributeToMaterial& matVec) :
 mesh_(mesh),
 attToMatVec_(matVec)
 {

diff --git a/src/maxwell/Model.h b/src/maxwell/Model.h
@@ -8,20 +8,20 @@ using namespace mfem;
 
 namespace maxwell {
 
-	using attribute = std::size_t;
+using Attribute = std::size_t;
+using AttributeToMaterial = std::map<Attribute, Material>;
 
 class Model {
 public:
 
-	Model(Mesh& mesh, std::vector<std::pair<attribute, Material>>& attToMatVec);
+	Model(Mesh& mesh, const AttributeToMaterial& attToMatVec);
 	Mesh& getMesh() { return mesh_; };
 	const Mesh& getConstMesh() const { return mesh_; };
-	const std::vector<std::pair<attribute, Material>>& getAttToMatVec() const { return attToMatVec_; }
+	const AttributeToMaterial& getAttToMat() const { return attToMatVec_; }
 
 private:
-
 	Mesh mesh_;
-	std::vector<std::pair<attribute, Material>> attToMatVec_;
+	AttributeToMaterial attToMatVec_;
 
 };
 

diff --git a/src/maxwell/Solver.cpp b/src/maxwell/Solver.cpp
@@ -133,9 +133,9 @@ const std::vector<std::vector<IntegrationPoint>> Solver::buildIntegrationPointsS
 const std::vector<std::vector<std::array<double, 3>>> Solver::saveFieldAtPointsForAllProbes()
 {
 	auto maxDir = model_.getConstMesh().Dimension();
-	std::vector<FieldByVDIM> res;
+	std::vector<FieldFrame> res;
 	for (int i = 0; i < probes_.getProbeVector().size(); i++) {
-		FieldByVDIM aux;
+		FieldFrame aux;
 		aux.resize(elemIds_.at(i).Size());
 		for (int j = 0; j < elemIds_.at(i).Size(); j++) {
 			for (int dir = Direction::X; dir != maxDir; dir++) {

diff --git a/src/maxwell/Solver.h b/src/maxwell/Solver.h
@@ -16,8 +16,11 @@ class Solver {
 public:
 
     using IntegrationPointsSet = std::vector<std::vector<IntegrationPoint>>;
-    using FieldByVDIM = std::vector<std::array<double, 3>>;
-    using TimeFieldPair = std::vector<std::pair<double, FieldByVDIM>>;
+
+    using Time = double;
+    using CVec3 = std::array<double, 3>;
+    using FieldFrame = std::vector<CVec3>;
+    using FieldMovie = std::map<Time, FieldFrame>;
 
     struct Options {
         int order = 2;
@@ -33,7 +36,7 @@ class Solver {
     const Vector& getMaterialProperties(const Material&) const;
 
     mfem::Mesh& getMesh() { return mesh_; }
-    std::vector<TimeFieldPair>& getFieldAtPoint() { return timeField_; }
+    std::vector<FieldMovie>& getFieldForProbe(std::size_t probeNum) { return timeField_; }
 
     void run();
 
@@ -63,8 +66,8 @@ class Solver {
     std::vector<IntegrationPointsSet> integPointSet_;
     std::vector<FieldType> fieldToExtract_;
     double timeRecord_;
-    std::vector<FieldByVDIM> fieldRecord_;
-    std::vector<std::vector<std::pair<double, FieldByVDIM>>> timeField_;
+    std::vector<FieldFrame> fieldRecord_;
+    std::vector<std::vector<std::pair<double, FieldFrame>>> timeField_;
 
     std::unique_ptr<mfem::ParaViewDataCollection> pd_;
 
@@ -74,7 +77,7 @@ class Solver {
 
     std::vector<std::pair<Array<int>, Array<IntegrationPoint>>> Solver::buildElemAndIntegrationPointArrays(DenseMatrix& physPoints);
     const IntegrationPointsSet Solver::buildIntegrationPointsSet(const Array<IntegrationPoint>& ipArray) const;
-    const std::vector<FieldByVDIM> saveFieldAtPointsForAllProbes();
+    const std::vector<FieldFrame> saveFieldAtPointsForAllProbes();
 
     void initializeParaviewData();
     //void initializeGLVISData();

diff --git a/src/maxwell/Types.h b/src/maxwell/Types.h
@@ -4,6 +4,8 @@
 
 namespace maxwell {
 
+using namespace mfem;
+
 using Position = mfem::Vector;
 
 enum FieldType {

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -8,8 +8,12 @@ include_directories(${MFEM_INCLUDE_DIRS})
 
 include_directories(${maxwell_solvers_INCLUDE_DIRS})
 
-add_executable(maxwell_tests "TestAux.cpp" "maxwell/TestSolver1D.cpp" "maxwell/TestBilinearIntegrators.cpp" "maxwell/TestMaterial.cpp")
+add_executable(maxwell_tests 
+"TestAux.cpp" 
+"maxwell/TestSolver1D.cpp" 
+"maxwell/TestBilinearIntegrators.cpp" 
+"maxwell/TestMaterial.cpp")
 
-target_link_libraries(maxwell_tests maxwell_solvers mfem GTest::gtest GTest::gtest_main) 
+target_link_libraries(maxwell_tests maxwell_solvers GTest::gtest GTest::gtest_main) 
 
 
diff --git a/test/TestAux.cpp b/test/TestAux.cpp
@@ -49,7 +49,7 @@ namespace HelperFunctions {
 
 	SparseMatrix* rotateMatrixLexico(BilinearForm& matrix)
 	{
-		const Operator* rotatorOperator = matrix.GetFES()->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+		const Operator* rotatorOperator = matrix.FESpace()->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
 		const SparseMatrix rotatorMatrix = HelperFunctions::operatorToSparseMatrix(rotatorOperator);
 		const SparseMatrix matrixSparse = matrix.SpMat();
 		SparseMatrix* res;
@@ -245,8 +245,8 @@ TEST_F(Auxiliary, checkMassMatrixIsSameForH1andDG)
 		ASSERT_EQ(rotatedMassMatrixH1Sparse->NumRows(), massMatrixDGSparse.NumRows());
 		ASSERT_EQ(rotatedMassMatrixH1Sparse->NumCols(), massMatrixDGSparse.NumCols());
 
-		for (std::size_t i = 0; i < massMatrixDGSparse.NumRows(); i++) {
-			for (std::size_t j = 0; j < massMatrixDGSparse.NumCols(); j++) {
+		for (int i = 0; i < massMatrixDGSparse.NumRows(); i++) {
+			for (int j = 0; j < massMatrixDGSparse.NumCols(); j++) {
 				EXPECT_NEAR(rotatedMassMatrixH1Sparse->Elem(i, j), massMatrixDGSparse.Elem(i, j), 1e-5);
 			}
 		}
@@ -375,8 +375,8 @@ TEST_F(Auxiliary, checkKOperators)
 	ASSERT_EQ(matrixK->NumRows(), matrixF->NumRows());
 	ASSERT_EQ(matrixK->NumCols(), matrixF->NumCols());
 
-	for (std::size_t i = 0; i < matrixK->NumRows(); i++) {
-		for (std::size_t j = 0; j < matrixK->NumCols(); j++) {
+	for (int i = 0; i < matrixK->NumRows(); i++) {
+		for (int j = 0; j < matrixK->NumCols(); j++) {
 			EXPECT_NEAR(matrixK->Elem(i, j), matrixS->Elem(i, j) + matrixF->Elem(i, j), tol);
 		}
 	}
@@ -485,7 +485,6 @@ TEST_F(Auxiliary, printGLVISDataForBasisFunctionNodes)
 
 	Vector nodalVector(order + 1);
 	Vector dofVector(order + 1);
-	IntegrationPoint integPoint;
 	Array<int> vdofs;
 
 	Mesh mesh = HelperFunctions::buildCartesianMeshForOneElement(1, Element::SEGMENT);

diff --git a/test/maxwell/TestSolver1D.cpp b/test/maxwell/TestSolver1D.cpp
@@ -67,11 +67,11 @@ namespace HelperFunctions {
 		return res;
 	}
 
-	std::vector<std::pair<attribute, Material>> buildAttToMatVec(const std::vector<attribute>& attVec, const std::vector<Material>& matVec)
+	AttributeToMaterial buildAttToMatVec(const std::vector<Attribute>& attVec, const std::vector<Material>& matVec)
 	{
-		std::vector<std::pair<attribute, Material>> res;
+		AttributeToMaterial res;
 		for (int i = 0; i < attVec.size(); i++) {
-			res.push_back(std::make_pair(attVec[i], matVec[i]));
+			res.emplace(attVec[i], matVec[i]);
 		}
 		return res;
 	}
@@ -88,12 +88,12 @@ class TestMaxwellSolver1D : public ::testing::Test {
 	Material mat11 = Material(1.0, 1.0); Material mat12 = Material(1.0, 2.0);
 	Material mat21 = Material(2.0, 1.0); Material mat22 = Material(2.0, 2.0);
 
-	std::vector<attribute> attArrSingle = std::vector<attribute>({ 1 });
-	std::vector<attribute> attArrMultiple = std::vector<attribute>({ 1, 2, 3, 4 });
+	std::vector<Attribute> attArrSingle = std::vector<Attribute>({ 1 });
+	std::vector<Attribute> attArrMultiple = std::vector<Attribute>({ 1, 2, 3, 4 });
 	std::vector<Material> matArrSimple = std::vector<Material>({ mat11 });
 	std::vector<Material> matArrMultiple = std::vector<Material>({ mat11,mat12,mat21,mat22 });
 
-	std::vector<std::pair<attribute, Material>> attToMatVec = HelperFunctions::buildAttToMatVec(attArrSingle, matArrSimple);
+	AttributeToMaterial attToMatVec = HelperFunctions::buildAttToMatVec(attArrSingle, matArrSimple);
 
 	Model testModel = Model(mesh1D, attToMatVec);