Remove the old skew code. Add possibility to use the shorter (longer)…

… diagonal when splitting 3D prismatic meshes into tets.
ElmerCSC · Dec 5, 2024 · b56d41b · b56d41b
1 parent 8c6b14c
commit b56d41b
Show file tree

Hide file tree

Showing 2 changed files with 60 additions and 113 deletions.
diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90
@@ -18446,9 +18446,7 @@ FUNCTION MeshExtrude(Mesh_in, Vlist) RESULT(Mesh_out)
       CALL AppendMissingBCs(CurrentModel,j)
     END IF
 
-    IF(.NOT. ListGetLogical(CurrentModel % Simulation,'Old Skew',Found ) ) THEN
-      CALL SetMeshSkew(Mesh_out, CurrentModel % Simulation )
-    END IF
+    CALL SetMeshSkew(Mesh_out, CurrentModel % Simulation )
 
     CALL PrepareMesh( CurrentModel, Mesh_out, isParallel )
 
@@ -19009,11 +19007,7 @@ FUNCTION MeshExtrudeSlices(Mesh_in, Vlist) RESULT(Mesh_out)
       DEALLOCATE(TmpBCs)
     END BLOCK
 
-
-    IF(.NOT. ListGetLogical(CurrentModel % Simulation,'Old Skew',Found ) ) THEN
-      CALL SetMeshSkew(Mesh_out, CurrentModel % Simulation )
-    END IF
-
+    CALL SetMeshSkew(Mesh_out, CurrentModel % Simulation )
 
     ExtrudedMeshName = ListGetString(CurrentModel % Simulation,'Extruded Mesh Name',Found)
     IF(Found) THEN
@@ -23906,6 +23900,41 @@ SUBROUTINE SplitMeshQuads(Mesh, Vlist)
     END DO
     CutCorner = MinCorner
 
+    BLOCK
+      INTEGER :: DoMax
+      INTEGER, POINTER :: Inds(:)
+      REAL(KIND=dp) :: d13,d24
+
+      DoMax = 0
+      IF(LIstGetLogical(Vlist,'Split Mesh Prisms Min',Found )) DoMax = 1
+      IF(LIstGetLogical(Vlist,'Split Mesh Prisms Max',Found )) DoMax = -1
+
+      ! Optionally cut the 3D meshes such that the shorter (longer) diagonal
+      ! is used to cut the quad faces. 
+      IF(DoMax /= 0) THEN
+        DO i=1,Mesh % NumberOfFaces
+          Face => Mesh % Faces(i)
+          IF(Face % TYPE % ElementCode /= 404) CYCLE
+          Inds => Face % NodeIndexes
+
+          ! Compute |r1-r3|^2 
+          d13 = (x(Inds(1))-x(Inds(3)))**2 + &
+              (y(Inds(1))-y(Inds(3)))**2 + (z(Inds(1))-z(Inds(3)))**2 
+
+          ! Compute |r2-r4|^2 
+          d24 = (x(Inds(2))-x(Inds(4)))**2 + &
+              (y(Inds(2))-y(Inds(4)))**2 + (z(Inds(2))-z(Inds(4)))**2 
+
+          IF(DoMax * d13 < DoMax * d24) THEN
+            CutCorner(i) = MIN(Inds(1),Inds(3))
+          ELSE
+            CutCorner(i) = MIN(Inds(2),Inds(4))
+          END IF
+        END DO
+      END IF
+    END BLOCK
+
+
 
     DO WHILE(.TRUE.)
       CutChanges = 0

diff --git a/fem/src/modules/RigidMeshMapper.F90 b/fem/src/modules/RigidMeshMapper.F90
@@ -70,8 +70,7 @@ SUBROUTINE RigidMeshMapper( Model,Solver,dt,Transient )
       RotorBodies(:), IntArray(:) => NULL()
   REAL(KIND=dp) :: x0(4), x1(4), RotMatrix(4,4),TrsMatrix(4,4),SclMatrix(4,4), &
       TrfMatrix(4,4),Identity(4,4), Origin(4),Angles(3),Scaling(3),alpha, Coord(3), &
-      dCoord(3), Norm, dx(3), zmin, zmax, zloc
-  REAL(KIND=dp) :: RotorSkew, StatorSkew
+      dCoord(3), Norm, dx(3) 
   REAL(KIND=dp) :: at0,at1,at2,Coeff,Source,relax(1),MaxDeform,AngleCoeff, RotorRad, RotorAngle
   REAL(KIND=dp), POINTER :: Xorig(:),Yorig(:),Zorig(:),Xnew(:),Ynew(:),Znew(:),&
       RelaxField(:),VeloVal(:), PArray(:,:) => NULL()
@@ -81,7 +80,7 @@ SUBROUTINE RigidMeshMapper( Model,Solver,dt,Transient )
   LOGICAL :: Found,GotMatrix,GotRotate,GotTranslate,GotScale,Visited=.FALSE.,&
       UseOriginalMesh, Cumulative, GotRelaxField=.FALSE., &
       CalculateVelocity,TranslateBeforeRotate, StoreOriginalMesh, &
-      RotorMode, WholeMode, DoIt, GotSkew, GotRotorAngle, GotSkewFun
+      RotorMode, WholeMode, DoIt, GotRotorAngle
   LOGICAL :: AnyMeshMatrix,AnyMeshRotate,AnyMeshTranslate,AnyMeshScale,&
       AnyMeshOrigin, AnyRelax, ConstantMap, GotMap, IsRotor
   LOGICAL, POINTER :: NodeDone(:)
@@ -150,28 +149,27 @@ SUBROUTINE RigidMeshMapper( Model,Solver,dt,Transient )
 
 
   IF( RotorMode .AND. .NOT. Visited ) THEN
-    IF(.NOT. ListGetLogical(CurrentModel % Simulation,'Old Skew',Found ) ) THEN
-      CALL SetMeshSkew(Mesh, CurrentModel % Simulation) 
-
-      ALLOCATE(RotorElement(Mesh % NumberOfBulkElements))
-      RotorElement = .FALSE.
-
-      DO elem = 1,Mesh % NumberOfBulkElements      
-        Element => Mesh % Elements(elem)
-        n = GetElementNOFNodes(Element)
-        CALL GetElementNodes( Nodes, Element )
-        Coord(1) = SUM(Nodes % x(1:n)) / n
-        Coord(2) = SUM(Nodes % y(1:n)) / n
-        Coord(3) = SUM(Nodes % z(1:n)) / n
-        IF(ASSOCIATED(RotorBodies)) THEN
-          IsRotor = ANY( RotorBodies == Element % BodyId ) 
-        ELSE
-          IsRotor = (Coord(1)**2+Coord(2)**2 < RotorRad**2) 
-        END IF
-        RotorElement(elem) = IsRotor
-      END DO
-      NodeDone = .FALSE.
-    END IF
+    ! Skew is moved to the library so that it can be done in different ways. 
+    CALL SetMeshSkew(Mesh, CurrentModel % Simulation) 
+
+    ALLOCATE(RotorElement(Mesh % NumberOfBulkElements))
+    RotorElement = .FALSE.
+
+    DO elem = 1,Mesh % NumberOfBulkElements      
+      Element => Mesh % Elements(elem)
+      n = GetElementNOFNodes(Element)
+      CALL GetElementNodes( Nodes, Element )
+      Coord(1) = SUM(Nodes % x(1:n)) / n
+      Coord(2) = SUM(Nodes % y(1:n)) / n
+      Coord(3) = SUM(Nodes % z(1:n)) / n
+      IF(ASSOCIATED(RotorBodies)) THEN
+        IsRotor = ANY( RotorBodies == Element % BodyId ) 
+      ELSE
+        IsRotor = (Coord(1)**2+Coord(2)**2 < RotorRad**2) 
+      END IF
+      RotorElement(elem) = IsRotor
+    END DO
+    NodeDone = .FALSE.
   END IF
 
   ! If using original mesh as a reference mesh it must be saved,
@@ -198,86 +196,6 @@ SUBROUTINE RigidMeshMapper( Model,Solver,dt,Transient )
     Zorig => Znew
   END IF
 
-  IF( RotorMode .AND. .NOT. Visited ) THEN
-    IF(ListGetLogical(CurrentModel % Simulation,'Old Skew',Found ) ) THEN
-      RotorSkew = AngleCoeff * ListGetCReal(CurrentModel % Simulation,'Rotor Skew',GotSkew )
-      GotSkewFun = ListCheckPresent( CurrentModel % Simulation,'Rotor Skew Function')
-      StatorSkew = AngleCoeff * ListGetCReal(CurrentModel % Simulation,'Stator Skew',Found )
-      GotSkew = GotSkew .OR. GotSkewFun .OR. Found
-      IF( GotSkew ) THEN
-        zmax = ListGetCReal( CurrentModel % Simulation,'Rotor Skew Max Coordinate',Found ) 
-        IF(.NOT. Found) THEN
-          zmax = ListGetCReal( CurrentModel % Simulation,'Extruded Max Coordinate',Found ) 
-          IF(.NOT. Found) zmax = ParallelReduction(MAXVAL(Zorig))        
-        END IF
-        IF(.NOT. Found) THEN
-          CALL Fatal(Caller,'"Rotor Skew" currently requires "Extruded Max Coordinate" to be given!')
-        END IF
-        zmin = ListGetCReal( CurrentModel % Simulation,'Rotor Skew Min Coordinate',Found ) 
-        IF(.NOT. Found) THEN
-          zmin = ListGetCReal( CurrentModel % Simulation,'Extruded Min Coordinate',Found ) 
-          IF(.NOT. Found) zmin = ParallelReduction(MINVAL(Zorig))
-        END IF
-        IF(InfoActive(20)) THEN
-          PRINT *,'RotorSkew:',RotorSkew, StatorSkew, zmin, zmax, GotSkew, GotSkewFun 
-        END IF
-      END IF
-
-      ALLOCATE(RotorElement(Mesh % NumberOfBulkElements))
-      RotorElement = .FALSE.
-
-      DO elem = 1,Mesh % NumberOfBulkElements      
-        Element => Mesh % Elements(elem)
-        n = GetElementNOFNodes(Element)
-
-        CALL GetElementNodes( Nodes, Element )
-        Coord(1) = SUM(Nodes % x(1:n)) / n
-        Coord(2) = SUM(Nodes % y(1:n)) / n
-        Coord(3) = SUM(Nodes % z(1:n)) / n
-        IF(ASSOCIATED(RotorBodies)) THEN
-          IsRotor = ANY( RotorBodies == Element % BodyId ) 
-        ELSE
-          IsRotor = (Coord(1)**2+Coord(2)**2 < RotorRad**2) 
-        END IF
-
-        RotorElement(elem) = IsRotor
-
-        IF(GotSkew) THEN
-          DO Node=1,n
-            NodeIndex(1) = Element % NodeIndexes(Node)
-            NodeI = NodeIndex(1)
-            IF(.NOT. NodeDone(NodeI)) THEN
-              Coord(1) = Xorig(NodeI)
-              Coord(2) = Yorig(NodeI)
-              Coord(3) = Zorig(NodeI)
-
-              ! Skew is not constant, perform it for each node 1st if requested. 
-              zloc = (coord(3)-zmin)/(zmax-zmin)
-
-              ! By construction this must be in [0,1]
-              zloc = MAX(0.0_dp,MIN(1.0_dp,zloc))
-
-              IF( IsRotor ) THEN
-                IF(GotSkewFun) THEN
-                  alpha = AngleCoeff * ListGetFun( CurrentModel % Simulation,'Rotor Skew Function',zloc)                
-                ELSE
-                  alpha = (zloc-0.5_dp) * RotorSkew
-                END IF
-              ELSE
-                alpha = (zloc-0.5_dp) * StatorSkew 
-              END IF
-
-              Xorig(NodeI) = Coord(1)*COS(alpha) - Coord(2)*SIN(alpha)
-              Yorig(Nodei) = Coord(1)*SIN(alpha) + Coord(2)*COS(alpha)        
-              NodeDone(NodeI) = .TRUE.
-            END IF
-          END DO
-        END IF
-      END DO
-      NodeDone = .FALSE.
-    END IF
-  END IF
-
 
   CalculateVelocity = GetLogical( SolverParams,'Calculate Mesh Velocity',Found)
   IF( CalculateVelocity ) THEN