Py bem

srcs/BEM/boundaryElements.cpp

-#include <gmsh.h>
-#include <iostream>
-#include <map>
-#include <sstream>
-#include <Eigen/Dense>
-#include <cmath>
-#include <iomanip>
-#include <stdio.h>
-#include <chrono>
 #include "functionsBEM.hpp"
-#include <omp.h>
-#include <algorithm>
-#include <list>
 
 // Fills the "boundaryConditions" map associating to each physical group (a string) the type of boundary condition ('true':
 // dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
@@ -44,12 +32,13 @@ void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
 }
 
 // Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the
-// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also fills the
-// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current
-// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of
-// 'dirichlet elements' (useful for the construction/resolution of the system).
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also returns
+// "domainPhysicalGroupDim" and "domainPhysicalGroupTag", the dimension and the tag (respectivelly) of the BEM domain of
+// interest. "domainName" is the name of the current BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement
+// to a serie of node tags. Returns the number of 'dirichlet elements' (useful for the construction/resolution of the system).
 std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
-                              std::map<std::string, std::pair<int, int>> &physicalGroups,
+                              int &domainPhysicalGroupDim,
+                              int &domainPhysicalGroupTag,
                               std::map<int, std::pair<double, double>> &nodalDisplacements,
                               const std::string domainName)
 {
@@ -68,11 +57,14 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
         int physicalGroupDim = physicalGroupsDimTags[i].first;
         int physicalGroupTag = physicalGroupsDimTags[i].second;
         std::string physicalGroupName;
-        gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
-        physicalGroups[physicalGroupName] = {physicalGroupDim, physicalGroupTag};
-
         gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
         
+        if(physicalGroupName.compare(domainName) == 0)
+        {
+            domainPhysicalGroupDim = physicalGroupDim;
+            domainPhysicalGroupTag = physicalGroupTag;
+        }
+
         // If no boundary condition on a given physical curve, don't do anything.
         if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end())
             continue;
@@ -103,7 +95,8 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
                 std::vector<std::size_t> nodeTags;
                 std::vector<double> nodeCoords;
                 std::vector<double> parametricNodeCoords;
-                gmsh::model::mesh::getNodesByElementType(elementTypes[j], nodeTags, nodeCoords, parametricNodeCoords, entityTags[i], false);
+                gmsh::model::mesh::getNodesByElementType(elementTypes[j], nodeTags, nodeCoords, parametricNodeCoords,
+                                                         entityTags[i], false);
 
                 // Constructs a map associating the index of each node tag in the "nodeCoords" vector.
                 std::map<std::size_t, std::size_t> nodeIndices;
@@ -156,8 +149,10 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
                     // the vector.
                     #pragma omp critical
                     {
-                        elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end());
-                        elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end());
+                        elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(),
+                                             localDirichletElementVector.end());
+                        elementVector.insert(elementVector.end(), localNeumannElementVector.begin(),
+                                             localNeumannElementVector.end());
                         nbDirichletElements += localDirichletElementVector.size();
                     }
                 }
@@ -291,7 +286,11 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
         }
     }
     if(!setEpsilon)
-        std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon << " [F/m].\n";
+    {
+        std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName;
+        std::cout << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon;
+        std::cout << " [F/m].\n";
+    }
 
     #pragma omp parallel for schedule(guided)
     for(size_t i = 0; i < elementVector.size(); i++)