Several small fixes.

ElmerCSC · Aug 27, 2024 · 757916e · 757916e
1 parent e82613a
commit 757916e
Show file tree

Hide file tree

Showing 6 changed files with 48 additions and 40 deletions.
diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90
@@ -10956,7 +10956,6 @@ SUBROUTINE AddProjectorWeakGeneric()
         END DO
       END IF
 
-
       DO ind=1,BMesh1 % NumberOfBulkElements
 
         ! Optionally save the submesh for specified element, for visualization and debugging
@@ -11139,7 +11138,6 @@ SUBROUTINE AddProjectorWeakGeneric()
         ElemCands = 0
         ElemHits = 0
 
-
         DO indM=1,BMesh2 % NumberOfBulkElements
 
           ElementM => BMesh2 % Elements(indM)        
@@ -12823,16 +12821,20 @@ FUNCTION SegmentOriginDistance2(x1,y1,x2,y2) RESULT ( r2 )
       REAL(KIND=dp) :: x1,y1,x2,y2,r2
       REAL(KIND=dp) :: q,xc,yc
 
-      q = ( x1*(x1-x2) + y1*(y1-y2) ) / &
-          SQRT((x1**2+y1**2) * ((x1-x2)**2+(y1-y2)**2))
-      IF( q <= 0.0_dp ) THEN
-        r2 = x1**2 + y1**2
-      ELSE IF( q >= 1.0_dp ) THEN
-        r2 = x2**2 + y2**2
+      IF(x1==0.AND.y1==0) THEN
+        r2 = 0;
       ELSE
-        xc = x1 + q * (x2-x1)
-        yc = y1 + q * (y2-y1)
-        r2 = xc**2 + yc**2
+        q = ( x1*(x1-x2) + y1*(y1-y2) ) / &
+            SQRT((x1**2+y1**2) * ((x1-x2)**2+(y1-y2)**2))
+        IF( q <= 0.0_dp ) THEN
+          r2 = x1**2 + y1**2
+        ELSE IF( q >= 1.0_dp ) THEN
+          r2 = x2**2 + y2**2
+        ELSE
+          xc = x1 + q * (x2-x1)
+          yc = y1 + q * (y2-y1)
+          r2 = xc**2 + yc**2
+        END IF
       END IF
     END FUNCTION SegmentOriginDistance2
 

diff --git a/fem/src/RadiationFactors.F90 b/fem/src/RadiationFactors.F90
@@ -1280,7 +1280,7 @@ SUBROUTINE CalculateGebhartFactors()
        END IF
 
        ! Scale matrix to unit diagonals
-       Diag = SQRT(1._dp/ABS(Diag))
+       Diag = SQRT(1._dp/MAX(ABS(Diag),1.0d-12))
        DO i=1,RadiationSurfaces
          DO j=G % Rows(i),G % Rows(i+1)-1
            G % Values(j) = G % Values(j)*Diag(i)*Diag(G % Cols(j))

diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90
@@ -762,7 +762,12 @@ SUBROUTINE UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
 
        np = mGetElementDOFs(pIndexes,Element)
        np = MIN(np, n)
-       NormalIndexes(1:np) = NT % BoundaryReorder(pIndexes(1:np))
+       DO i=1,np
+         j = pIndexes(i)
+         IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN
+           NormalIndexes(i) = NT % BoundaryReorder(j)
+         END IF
+       END DO
 
        CALL RotateMatrix( LocalStiffMatrix, LocalForce, n, dim, NDOFs, &
           NormalIndexes, NT % BoundaryNormals, NT % BoundaryTangent1, NT % BoundaryTangent2 )
@@ -816,7 +821,7 @@ SUBROUTINE UpdateGlobalEquationsVec( Gmtr, Lmtr, Gvec, Lvec, n, &
 
      ! Local variables
      INTEGER :: dim, i,j,k,np
-     INTEGER :: Ind(n*NDOFs),pIndexes(64)
+     INTEGER :: NormalIndexes(n),pIndexes(64)
      REAL(KIND=dp) :: Vals(n*NDOFs)
 !DIR$ ATTRIBUTES ALIGN:64::Ind, Vals
 
@@ -848,22 +853,20 @@ SUBROUTINE UpdateGlobalEquationsVec( Gmtr, Lmtr, Gvec, Lvec, n, &
      END IF
 
      IF ( Rotate ) THEN
-       Ind = 0
+       NormalIndexes = 0
 
        np = mGetElementDOFs(pIndexes,Element)
        np = MIN(np,n)
 
        DO i=1,np
          j = pIndexes(i)
-         IF(j > 0 .AND. j <= SIZE(NT % BoundaryReorder) ) THEN
-           Ind(i) = NT % BoundaryReorder(j)
-         ELSE
-           Ind(i) = j
+         IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN
+           NormalIndexes(i) = NT % BoundaryReorder(j)
          END IF
        END DO
 
        ! TODO: See that RotateMatrix is vectorized
-       CALL RotateMatrix( Lmtr, Lvec, n, dim, NDOFs, Ind, NT % BoundaryNormals, &
+       CALL RotateMatrix( Lmtr, Lvec, n, dim, NDOFs, NormalIndexes, NT % BoundaryNormals, &
                     NT % BoundaryTangent1, NT % BoundaryTangent2 )
 
        !IF ( Rotate .AND. NormalTangentialNOFNodes > 0 .AND. ndofs>=dim) THEN
@@ -989,7 +992,13 @@ SUBROUTINE UpdateGlobalForce(ForceVector, LocalForce, n, &
 
        np = mGetElementDOFs(pIndexes,Element)
        np = MIN(np,n)
-       NormalIndexes(1:np) = NT % BoundaryReorder(pIndexes(1:np))
+
+       DO i=1,np
+         j = pIndexes(i)
+         IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN
+           NormalIndexes(i) = NT % BoundaryReorder(j)
+         END IF
+       END DO
 
        CALL RotateMatrix( LocalStiffMatrix, LocalForce, n, dim, NDOFs, &
           NormalIndexes, NT % BoundaryNormals, NT % BoundaryTangent1, NT % BoundaryTangent2 )
@@ -5017,8 +5026,7 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
             Conditional = ActiveCond(BC)
 
             Element => Mesh % Elements(t)
-            IF ( Element % BoundaryInfo % Constraint /= &
-                Model % BCs(BC) % Tag ) CYCLE
+            IF ( Element % BoundaryInfo % Constraint /= Model % BCs(BC) % Tag ) CYCLE
 
             Model % CurrentElement => Element
             IF ( ActivePart(BC) ) THEN
@@ -5045,8 +5053,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
       END IF
     END IF
 
-
-
     BLOCK
       INTEGER,ALLOCATABLE :: ChildBCs(:)
       INTEGER,POINTER :: BCInds(:)
@@ -5215,7 +5221,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
       END IF
     END DO
 
-
     ! Check the boundaries and body forces for possible single nodes BCs that are used to fixed 
     ! the domain for undetermined equations. The loop is slower than optimal in the case that there is 
     ! a large amount of different boundaries that have a node to set. 
@@ -5321,7 +5326,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
       END DO
     END IF
 
-
 !------------------------------------------------------------------------------
 !   Take care of the matrix entries of passive elements
 !------------------------------------------------------------------------------
@@ -5376,7 +5380,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
       DEALLOCATE(PassPerm,NodeIndexes)
     END IF
 
-
     ! Check the boundaries and body forces for possible single nodes BCs that must have a constant
     ! value on that boundary / body force.
     !--------------------------------------------------------------------------------------------
@@ -5578,7 +5581,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
           //I2S(SIZE( A % Cols )),Level=8)
     END IF
 
-
     ! Check the boundaries for a curve bc.
     !--------------------------------------------------------------------------------------------
     BLOCK
@@ -5781,7 +5783,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
 
     END BLOCK
 
-
     ! Check the boundaries and body forces for possible single nodes BCs that must have a constant
     ! value on that boundary / body force.
     !--------------------------------------------------------------------------------------------
@@ -5887,7 +5888,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, &
           //I2S(DirCount),Level=12)
     END IF
 
-
 !------------------------------------------------------------------------------
 
   CONTAINS

diff --git a/fem/src/modules/MagnetoDynamics2D.F90 b/fem/src/modules/MagnetoDynamics2D.F90
@@ -3308,7 +3308,7 @@ SUBROUTINE BulkAssembly()
           BAtIp(8) = AIMAG(imag_value)
           imag_value = CMPLX(BatIp(1), BatIp(3), KIND=dp)
           imag_value2 = CMPLX(BatIp(2), BatIp(4), KIND=dp)
-          BMagnAtIP = SQRT(ABS(imag_value**2._dp) + ABS(imag_value2**2._dp))
+          BMagnAtIP = SQRT(ABS(imag_value*imag_value) + ABS(imag_value2*imag_value2))
         END IF
 
         IF (LorentzForceCompute) THEN
@@ -3320,8 +3320,14 @@ SUBROUTINE BulkAssembly()
           By = CMPLX(BatIp(2), BatIp(4), KIND=dp)
           Jz = CMPLX(BatIp(7), BatIp(8), KIND=dp)
 
-          LorentzForceDensX = ModelDepth * Weight * By / Jz * ABS(Jz)**2._dp
-          LorentzForceDensY = -ModelDepth * Weight * Bx / Jz * ABS(Jz)**2._dp
+          IF(Jz/=0) THEN
+            LorentzForceDensX = ModelDepth * Weight * By / Jz * ABS(Jz)*ABS(Jz)
+            LorentzForceDensY = -ModelDepth * Weight * Bx / Jz * ABS(Jz)*ABS(Jz)
+          ELSE
+            LorentzForceDensX = 0
+            LorentzForceDensY = 0
+          END IF
+
           BodyLorentzForcesRe(1, BodyId) = BodyLorentzForcesRe(1, BodyId) + &
             REAL(LorentzForceDensX)
           BodyLorentzForcesRe(2, BodyId) = BodyLorentzForcesRe(2, BodyId) + &

diff --git a/fem/tests/mgdyn_anisotropic_cond/currents.f90 b/fem/tests/mgdyn_anisotropic_cond/currents.f90
@@ -57,7 +57,7 @@ FUNCTION currdens1( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidR(r, r0, y, turns, I) * sin(theta)
 
@@ -87,7 +87,7 @@ FUNCTION currdens2( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidY(r, r0, y, turns, I)
 
@@ -117,7 +117,7 @@ FUNCTION currdens3( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidR(r, r0, y, turns, I) * cos(theta)
 

diff --git a/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90 b/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90
@@ -57,7 +57,7 @@ FUNCTION currdens1( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidR(r, r0, y, turns, I) * sin(theta)
 
@@ -87,7 +87,7 @@ FUNCTION currdens2( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidY(r, r0, y, turns, I)
 
@@ -117,7 +117,7 @@ FUNCTION currdens3( model, n, args) RESULT(curr)
 
   r = sqrt(x**2 + z**2)
 
-  theta = atan(x/z)
+  theta = atan2(x,z)
 
   curr = currentInToroidR(r, r0, y, turns, I) * cos(theta)