diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
index 45919f256dcb71344900469bc0dd09e5b3b704c6..7977f8348bc33dcca55563884bed8e6405ae62e8 100644
--- a/srcs/BEM/postProcessing.cpp
+++ b/srcs/BEM/postProcessing.cpp
@@ -317,47 +317,42 @@ void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
                          const int domainPhysicalGroupTag)
 {
     // #pragma omp parallel
+    std::vector<std::size_t> localElementTags;  // For parallelization.
+
+    // Loops over each boundary element and retrieves the element in contact with its boundary nodes.
+    // #pragma omp for schedule(guided)
+    for(std::size_t i = 0; i < elementVector.size(); i++)
     {
-        std::vector<std::size_t> localElementTags;  // For parallelization.
+        std::vector<size_t> elementTags;
+        gmsh::model::mesh::getElementsByCoordinates(elementVector[i].x1, elementVector[i].y1, 0, elementTags, 2);
+
 
-        // Loops over each boundary element and retrieves the element in contact with its boundary nodes.
-        // #pragma omp for schedule(guided)
-        for(std::size_t i = 0; i < elementVector.size(); i++)
+        // Adds the elements tags to the local vector. 
+        for(std::size_t j = 0; j < elementTags.size(); j++)
         {
+            // Retrieves the physical group to which the element belongs because it is possible that the element belongs to
+            // an adjacent physical surface.
             int elementType;
             std::vector<size_t> nodeTags;
             int entityDim;
             int entityTag;
             std::vector<int> physicalTags;
-
-            // Adds the elements tags to the local vector. 
-            for(std::size_t j = 0; j < elementTags.size(); j++)
+            // #pragma omp critical
             {
-                // Retrieves the physical group to which the element belongs because it is possible that the element belongs to
-                // an adjacent physical surface.
-                int elementType;
-                std::vector<size_t> nodeTags;
-                int entityDim;
-                int entityTag;
-                std::vector<int> physicalTags;
-                // #pragma omp critical
-                {
-                    gmsh::model::mesh::getElement(elementTags[j], elementType, nodeTags, entityDim, entityTag);
-                    gmsh::model::getPhysicalGroupsForEntity(entityDim, entityTag, physicalTags);
-                }
-                
-                // Add the element tag (to a local 'thread safe' vector) if the physical tag is the right one.
-                for(std::size_t k = 0; k < physicalTags.size(); k++)
-                    if(physicalTags[k] == domainPhysicalGroupTag)
-                        localElementTags.push_back(elementTags[j]);
+                gmsh::model::mesh::getElement(elementTags[j], elementType, nodeTags, entityDim, entityTag);
+                gmsh::model::getPhysicalGroupsForEntity(entityDim, entityTag, physicalTags);
             }
+            
+            // Add the element tag (to a local 'thread safe' vector) if the physical tag is the right one.
+            for(std::size_t k = 0; k < physicalTags.size(); k++)
+                if(physicalTags[k] == domainPhysicalGroupTag)
+                    localElementTags.push_back(elementTags[j]);
         }
-
-        // // Fills the "domainBoundaryElementTags" vector.
-        // #pragma omp critical
-        domainBoundaryElementTags.insert(domainBoundaryElementTags.end(), localElementTags.begin(), localElementTags.end());
     }
-    
+
+    // // Fills the "domainBoundaryElementTags" vector.
+    // #pragma omp critical
+    domainBoundaryElementTags.insert(domainBoundaryElementTags.end(), localElementTags.begin(), localElementTags.end());
 
     // Allows duplicate tags to be removed. If an element has several nodes in common with the boundary, it has been counted
     // several times before.