diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 3fbf5c8..0237de1 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -75,10 +75,10 @@ MODULE moduleRandom REAL(8):: rnd0b INTEGER:: rnd, i - rnd0b = random(0.D0, sumWeight) + rnd0b = random() i = 1 DO - IF (rnd0b <= cumWeight(i)) THEN + IF (rnd0b <= cumWeight(i)/sumWeight) THEN rnd = i EXIT diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index dfa3035..dbe182f 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = r2(2)**2 - r1(2)**2 + self%weight = DABS(self%r(2)**2 - self%r(1)**2) !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -223,21 +223,12 @@ MODULE moduleMesh2DCyl CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) - REAL(8):: dr, dz + REAL(8):: p1(1:2), p2(1:2) rnd = random() - dr = self%r(2) - self%r(1) - dz = self%z(2) - self%z(1) - IF (dr /= 0.D0) THEN - r(2) = dr * DSQRT(rnd) + self%r(1) - r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) - - ELSE - r(2) = self%r(1) - r(1) = dz * rnd + self%z(1) - - END IF - + 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 END FUNCTION randPosEdge @@ -572,7 +563,7 @@ MODULE moduleMesh2DCyl !Compute element volume PURE SUBROUTINE volumeQuad(self) - USE moduleConstParam, ONLY: PI8 + USE moduleConstParam, ONLY: PI, PI8 IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(inout):: self @@ -596,15 +587,31 @@ MODULE moduleMesh2DCyl Xi = (/-5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 + IF (self%r(1) == 0.D0) THEN + self%n1%v = self%n1%v * PI / 2.D0 + + END IF Xi = (/ 5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 + IF (self%r(2) == 0.D0) THEN + self%n2%v = self%n2%v * PI / 2.D0 + + END IF Xi = (/ 5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 + IF (self%r(3) == 0.D0) THEN + self%n3%v = self%n3%v * PI / 2.D0 + + END IF Xi = (/-5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 + IF (self%r(4) == 0.D0) THEN + self%n4%v = self%n4%v * PI / 2.D0 + + END IF END SUBROUTINE volumeQuad