diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 593af98..077e9d6 100644 --- a/src/modules/init/moduleInput.f90 +++ b/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 diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 1240c03..25a3fd7 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/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 diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index e250bb2..ca09561 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/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 diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 7a84d15..301cc6b 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/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):: r(1:3) - REAL(8):: p1(1:2), p2(1:2) + real(8):: Xi(1:3) + real(8):: r(1:3) + real(8):: fPsi(1:2) - rnd = random() + Xi = 0.d0 + Xi(1) = random() - p1 = (/self%x(1), self%y(1) /) - p2 = (/self%x(2), self%y(2) /) - r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 - r(3) = 0.D0 + fPsi = self%fPsi(Xi, 2) + + r = (/dot_product(fPsi, self%x), & + 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. & - (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0) + eps = 1.0d-10 + + + 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 diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index f14b2c8..0610b87 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/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):: r(1:3) - REAL(8):: p1(1:2), p2(1:2) + real(8):: Xi(1:3) + real(8):: r(1:3) + real(8):: fPsi(1:2) - rnd = random() + Xi = 0.d0 + Xi(1) = random() - p1 = (/self%z(1), self%r(1) /) - p2 = (/self%z(2), self%r(2) /) - r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 - r(3) = 0.D0 + fPsi = self%fPsi(Xi, 2) + + r = (/dot_product(fPsi, self%z), & + 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. & - (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0) + eps = 1.0d-10 + + 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 diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index a89f51c..e7cc060 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/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 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index db2a399..179097f 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/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 diff --git a/src/modules/mesh/moduleMesh@elements.f90 b/src/modules/mesh/moduleMesh@elements.f90 index 93ea8ff..6bd8623 100644 --- a/src/modules/mesh/moduleMesh@elements.f90 +++ b/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 diff --git a/src/modules/mesh/moduleMeshCommon.f90 b/src/modules/mesh/moduleMeshCommon.f90 index d2103b0..ee5f13b 100644 --- a/src/modules/mesh/moduleMeshCommon.f90 +++ b/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