Added center point procedure to edge and eps in inside to avoid some rounding errors

src/modules/init/moduleInput.f90

@@ -1255,7 +1255,7 @@ MODULE moduleInput
         findCellColl => findCellCollMesh
         ! Link edges with cells in meshColl
         DO e=1, mesh%numEdges
-          nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
+          nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%center())
           IF (nVolColl > 0) THEN
             mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj
 

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -28,6 +28,7 @@ MODULE moduleMesh1DCart
       PROCEDURE, PASS::   getNodes     => getNodes1DCart
       PROCEDURE, PASS::   intersection => intersection1DCart
       PROCEDURE, PASS::   randPos      => randPosEdge
+      procedure, pass::   center       => centerEdgePoint
       procedure, nopass:: fPsi         => fPsiPoint
 
   END TYPE meshEdge1DCart
@@ -164,6 +165,8 @@ MODULE moduleMesh1DCart
 
     !Calculate a 'random' position in edge
     FUNCTION randPosEdge(self) RESULT(r)
+      implicit none
+
       CLASS(meshEdge1DCart), INTENT(in):: self
       REAL(8):: r(1:3)
 
@@ -171,6 +174,16 @@ MODULE moduleMesh1DCart
 
     END FUNCTION randPosEdge
 
+    function centerEdgePoint(self) RESULT(r)
+      implicit none
+    
+      class(meshEdge1DCart), intent(in):: self
+      real(8):: r(1:3)
+
+      r = (/ self%x, 0.D0, 0.D0 /)
+
+    end function centerEdgePoint
+
     !CELL FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init element

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -28,6 +28,7 @@ MODULE moduleMesh1DRad
       PROCEDURE, PASS::   getNodes     => getNodes1DRad
       PROCEDURE, PASS::   intersection => intersection1DRad
       PROCEDURE, PASS::   randPos      => randPos1DRad
+      procedure, pass::   center       => centerEdgePoint
       procedure, nopass:: fPsi         => fPsiPoint
 
   END TYPE meshEdge1DRad
@@ -171,6 +172,14 @@ MODULE moduleMesh1DRad
 
     END FUNCTION randPos1DRad
 
+    function centerEdgePoint(self) RESULT(r)
+      class(meshEdge1DRad), intent(in):: self
+      real(8):: r(1:3)
+
+      r = (/ self%r, 0.D0, 0.D0 /)
+
+    end function centerEdgePoint
+
     !CELL FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init element

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -28,6 +28,7 @@ MODULE moduleMesh2DCart
       PROCEDURE, PASS::   getNodes     => getNodes2DCart
       PROCEDURE, PASS::   intersection => intersection2DCartEdge
       PROCEDURE, PASS::   randPos      => randPosEdge
+      procedure, pass::   center       => centerEdgeSegm
       procedure, nopass:: fPsi         => fPsiSegm
 
   END TYPE meshEdge2DCart
@@ -216,19 +217,37 @@ MODULE moduleMesh2DCart
       IMPLICIT NONE
 
       CLASS(meshEdge2DCart), INTENT(in):: self
-      REAL(8):: rnd
+      real(8):: Xi(1:3)
-      REAL(8):: r(1:3)
+      real(8):: r(1:3)
-      REAL(8):: p1(1:2), p2(1:2)
+      real(8):: fPsi(1:2)
 
-      rnd = random()
+      Xi    = 0.d0
+      Xi(1) = random()
 
-      p1 = (/self%x(1), self%y(1) /)
+      fPsi = self%fPsi(Xi, 2)
-      p2 = (/self%x(2), self%y(2) /)
+
-      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
+      r = (/dot_product(fPsi, self%x), &
-      r(3)   =  0.D0
+            dot_product(fPsi, self%y), &
+            0.d0/)
 
     END FUNCTION randPosEdge
 
+    function centerEdgeSegm(self) result(r)
+      use moduleMeshCommon, only: cenSeg
+      implicit none
+
+      class(meshEdge2DCart), intent(in):: self
+      real(8):: r(1:3)
+      real(8):: fPsi(1:2)
+
+      fPsi = self%fPsi(cenSeg, 2)
+
+      r = (/dot_product(fPsi, self%x), &
+            dot_product(fPsi, self%y), &
+            0.d0/)
+
+    end function centerEdgeSegm
+
     !VOLUME FUNCTIONS
     !QUAD FUNCTIONS
     !Init element
@@ -440,9 +459,13 @@ MODULE moduleMesh2DCart
 
       REAL(8), INTENT(in):: Xi(1:3)
       LOGICAL:: ins
+      real(8):: eps
 
-      ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. &
+      eps = 1.0d-10
-            (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0)
+
+
+      ins = (Xi(1) >= -1.D0-eps .AND. Xi(1) <= 1.D0+eps) .AND. &
+            (Xi(2) >= -1.D0-eps .AND. Xi(2) <= 1.D0+eps)
 
     END FUNCTION insideQuad
 

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -28,6 +28,7 @@ MODULE moduleMesh2DCyl
       PROCEDURE, PASS::   getNodes     => getNodes2DCyl
       PROCEDURE, PASS::   intersection => intersection2DCylEdge
       PROCEDURE, PASS::   randPos      => randPosEdge
+      procedure, pass::   center       => centerEdgeSegm
       procedure, nopass:: fPsi         => fPsiSegm
 
   END TYPE meshEdge2DCyl
@@ -225,19 +226,37 @@ MODULE moduleMesh2DCyl
       IMPLICIT NONE
 
       CLASS(meshEdge2DCyl), INTENT(in):: self
-      REAL(8):: rnd
+      real(8):: Xi(1:3)
-      REAL(8):: r(1:3)
+      real(8):: r(1:3)
-      REAL(8):: p1(1:2), p2(1:2)
+      real(8):: fPsi(1:2)
 
-      rnd = random()
+      Xi    = 0.d0
+      Xi(1) = random()
 
-      p1 = (/self%z(1), self%r(1) /)
+      fPsi = self%fPsi(Xi, 2)
-      p2 = (/self%z(2), self%r(2) /)
+
-      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
+      r = (/dot_product(fPsi, self%z), &
-      r(3)   =  0.D0
+            dot_product(fPsi, self%r), &
+            0.d0/)
 
     END FUNCTION randPosEdge
 
+    function centerEdgeSegm(self) result(r)
+      use moduleMeshCommon, only: cenSeg
+      implicit none
+
+      class(meshEdge2DCyl), intent(in):: self
+      real(8):: r(1:3)
+      real(8):: fPsi(1:2)
+
+      fPsi = self%fPsi(cenSeg, 2)
+
+      r = (/dot_product(fPsi, self%z), &
+            dot_product(fPsi, self%r), &
+            0.d0/)
+
+    end function centerEdgeSegm
+
     !CELL FUNCTIONS
     !QUAD FUNCTIONS
     !Init element
@@ -456,9 +475,12 @@ MODULE moduleMesh2DCyl
 
       REAL(8), INTENT(in):: Xi(1:3)
       LOGICAL:: ins
+      real(8):: eps
 
-      ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. &
+      eps = 1.0d-10
-            (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0)
+
+      ins = (Xi(1) >= -1.D0-eps .AND. Xi(1) <= 1.D0+eps) .AND. &
+            (Xi(2) >= -1.D0-eps .AND. Xi(2) <= 1.D0+eps)
 
     END FUNCTION insideQuad
 

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -29,6 +29,7 @@ MODULE moduleMesh3DCart
       PROCEDURE, PASS::   getNodes     => getNodes3DCartTria
       PROCEDURE, PASS::   intersection => intersection3DCartTria
       PROCEDURE, PASS::   randPos      => randPosEdgeTria
+      procedure, pass::   center       => centerEdgeTria
       !PARTICULAR PROCEDURES
       PROCEDURE, NOPASS:: fPsi => fPsiTria
 
@@ -214,6 +215,22 @@ MODULE moduleMesh3DCart
 
     END FUNCTION randPosEdgeTria
 
+    function centerEdgeTria(self) result(r)
+      use moduleMeshCommon, only: cenTria
+      implicit none
+
+      class(meshEdge3DCartTria), intent(in):: self
+      real(8):: r(1:3)
+      real(8):: fPsi(1:3)
+
+      fPsi = self%fPsi(cenTria, 3)
+
+      r = (/dot_product(fPsi, self%x), &
+            dot_product(fPsi, self%y), &
+            dot_product(fPsi, self%z)/)
+
+    end function centerEdgeTria
+
     !VOLUME FUNCTIONS
     !TETRA FUNCTIONS
     !Init element

src/modules/mesh/moduleMesh.f90

@@ -119,6 +119,7 @@ MODULE moduleMesh
       PROCEDURE(getNodesEdge_interface),     DEFERRED, PASS:: getNodes
       PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
       PROCEDURE(randPosEdge_interface),      DEFERRED, PASS:: randPos
+      procedure(centerEdge_interface),       deferred, pass:: center
       PROCEDURE(fPsi_interface),             DEFERRED, NOPASS:: fPsi
       !Gather value and spatial derivative on the nodes at position Xi
       PROCEDURE, PASS, PRIVATE:: gatherF_edge_scalar
@@ -177,6 +178,14 @@ MODULE moduleMesh
 
     END FUNCTION randPosEdge_interface
 
+    ! Returns the center point of an edge
+    FUNCTION centerEdge_interface(self) RESULT(r)
+      IMPORT:: meshEdge
+      CLASS(meshEdge), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+    END FUNCTION centerEdge_interface
+
   END INTERFACE
 
   !Containers for edges in the mesh

src/modules/mesh/moduleMesh@elements.f90

@@ -443,6 +443,9 @@ submodule(moduleMesh) elements
 
       DO e = 1, self%numCells
         Xi = self%cells(e)%obj%phy2log(r)
+        if (e < 160) then
+          print *, e, Xi
+        end if
         IF(self%cells(e)%obj%inside(Xi)) THEN
           nVol = self%cells(e)%obj%n
           EXIT

src/modules/mesh/moduleMeshCommon.f90

@@ -5,15 +5,23 @@ module moduleMeshCommon
   real(8), parameter:: corSeg(1:3)  = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
   real(8), parameter:: wSeg(1:3)    = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
 
-  ! Quad
-  real(8), parameter:: corQuad(1:3) = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
-  real(8), parameter:: wQuad(1:3)   = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
-
   ! tria
   real(8), parameter:: Xi1Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       3.D0/5.D0,   1.D0/5.D0  /)
   real(8), parameter:: Xi2Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       1.D0/5.D0,   3.D0/5.D0  /)
   real(8), parameter:: wTria(1:4)   = (/      -27.D0/96.D0, 25.D0/96.D0,     25.D0/96.D0, 25.D0/96.D0 /)
 
+  ! Quad
+  real(8), parameter:: corQuad(1:3) = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
+  real(8), parameter:: wQuad(1:3)   = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
+
+  ! Center point in natural coordinates
+  ! Segment
+  real(8), parameter:: cenSeg(1:3)  = (/ 0.5d0,       0.0d0,       0.0d0 /)
+  ! Tria
+  real(8), parameter:: cenTria(1:3) = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 0.0d0 /)
+  ! Quad
+  real(8), parameter:: cenQuad(1:3) = (/ 0.0d0,       0.0d0,       0.0d0 /)
+
   contains
     ! ELEMENT FUNCTIONS
     ! Point