Save Metainfo in VTU files by default. Minor fix in Particle stuff.

fem/src/ParticleUtils.F90

@@ -2447,7 +2447,7 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
 
     TYPE(ValueList_t), POINTER :: Params, BodyForce 
     TYPE(Variable_t), POINTER :: VeloVar, MaskVar, AdvVar
-    TYPE(Element_t), POINTER :: CurrentElement
+    TYPE(Element_t), POINTER :: Element
     TYPE(Mesh_t), POINTER :: Mesh
     TYPE(Nodes_t) :: Nodes
     INTEGER :: Offset, NewParticles,LastParticle,NoElements
@@ -2603,9 +2603,9 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
 
           j = 0
           DO i=1,Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
-            CurrentElement => Mesh % Elements(i)
-            NodeIndexes =>  CurrentElement % NodeIndexes
-            n = CurrentElement % TYPE % NumberOfNodes
+            Element => Mesh % Elements(i)
+            NodeIndexes =>  Element % NodeIndexes
+            n = Element % TYPE % NumberOfNodes
 
             IF( i == Mesh % NumberOfBulkElements ) THEN
               IF( j > 0 ) EXIT
@@ -2824,15 +2824,15 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         j = i
         IF( GotMask ) j = InvPerm(j)
 
-        CurrentElement => Mesh % Elements(j)
-        NodeIndexes =>  CurrentElement % NodeIndexes
-        n = CurrentElement % TYPE % NumberOfNodes
+        Element => Mesh % Elements(j)
+        NodeIndexes =>  Element % NodeIndexes
+        n = Element % TYPE % NumberOfNodes
 
         ! If weight is used see that we have a weight, and that it is positive
         IF( GotWeight ) THEN
           IF( j > Mesh % NumberOfBulkElements ) CYCLE
 
-          body_id = CurrentElement % BodyId
+          body_id = Element % BodyId
           bf_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values,&
               'Body Force',minv=1)
           BodyForce => CurrentModel % BodyForces(bf_id) % Values
@@ -2851,7 +2851,7 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         Nodes % y(1:n) = Mesh % Nodes % y(NodeIndexes)
         Nodes % z(1:n) = Mesh % Nodes % z(NodeIndexes)
 
-        DetJ = ElementSize( CurrentElement, Nodes ) 
+        DetJ = ElementSize( Element, Nodes ) 
         MaxDetJ = MAX( MaxDetJ, DetJ ) 
         MinDetJ = MIN( MinDetJ, DetJ )
       END DO      
@@ -2885,23 +2885,23 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         j = CEILING( NoElements * EvenRandom() )
         IF( GotMask ) j = InvPerm(j)
 
-        CurrentElement => Mesh % Elements(j)
-        NodeIndexes =>  CurrentElement % NodeIndexes
-        n = CurrentElement % TYPE % NumberOfNodes
+        Element => Mesh % Elements(j)
+        NodeIndexes =>  Element % NodeIndexes
+        n = Element % TYPE % NumberOfNodes
 
         Nodes % x(1:n) = Mesh % Nodes % x(NodeIndexes)
         Nodes % y(1:n) = Mesh % Nodes % y(NodeIndexes)
         Nodes % z(1:n) = Mesh % Nodes % z(NodeIndexes)
 
         IF( CheckForSize ) THEN
-          DetJ = ElementSize( CurrentElement, Nodes ) 
+          DetJ = ElementSize( Element, Nodes ) 
 
           ! The weight could be computed really using the integration point
           ! Here we assumes constant weight within the whole element. 
           IF( GotWeight ) THEN
             IF( j > Mesh % NumberOfBulkElements ) CYCLE
 
-            body_id = CurrentElement % BodyId
+            body_id = Element % BodyId
             bf_id = ListGetInteger( CurrentModel % Bodies(body_id) % Values,&
                 'Body Force',minv=1)
             BodyForce => CurrentModel % BodyForces(bf_id) % Values
@@ -2920,7 +2920,7 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         END IF
 
         ! Create a random particle within the element
-        Coord = RandomPointInElement( CurrentElement, Nodes )
+        Coord = RandomPointInElement( Element, Nodes )
 
         i = i + 1
         k = Offset + i
@@ -2950,9 +2950,9 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
           PRINT *,'j too large',j,i,k,(NoElements-1)*(i-1)/(NewParticles-1)+1
         END IF
 
-        CurrentElement => Mesh % Elements(j)
-        NodeIndexes =>  CurrentElement % NodeIndexes
-        n = CurrentElement % TYPE % NumberOfNodes
+        Element => Mesh % Elements(j)
+        NodeIndexes =>  Element % NodeIndexes
+        n = Element % TYPE % NumberOfNodes
         Coord(1) = SUM( Mesh % Nodes % x(NodeIndexes ) ) / n
         Coord(2) = SUM( Mesh % Nodes % y(NodeIndexes ) ) / n
         IF( dim == 3 ) Coord(3) = SUM( Mesh % Nodes % z(NodeIndexes ) ) / n
@@ -2985,9 +2985,9 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
 
       IF( SaveParticleOrigin ) THEN
         DO i=1,Mesh % NumberOfBulkElements
-          CurrentElement => Mesh % Elements(i)
-          NodeIndexes =>  CurrentElement % NodeIndexes
-          n = CurrentElement % TYPE % NumberOfNodes
+          Element => Mesh % Elements(i)
+          NodeIndexes =>  Element % NodeIndexes
+          n = Element % TYPE % NumberOfNodes
           DO j=1,n
             IF( ASSOCIATED( AdvVar ) ) THEN
               k = AdvVar % Perm(NodeIndexes(j))
@@ -3006,9 +3006,9 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         No = AdvVar % Perm( i )
         IF( No == 0 ) CYCLE
 
-        CurrentElement => Mesh % Elements(i)
-        NodeIndexes =>  CurrentElement % NodeIndexes
-        n = CurrentElement % TYPE % NumberOfNodes
+        Element => Mesh % Elements(i)
+        NodeIndexes =>  Element % NodeIndexes
+        n = Element % TYPE % NumberOfNodes
 
         Center(1) = SUM( Mesh % Nodes % x(NodeIndexes ) ) / n
         Center(2) = SUM( Mesh % Nodes % y(NodeIndexes ) ) / n
@@ -3043,9 +3043,9 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
           No = AdvVar % Perm( i )
           IF( No == 0 ) CYCLE
 
-          CurrentElement => Mesh % Elements(i)
-          NodeIndexes =>  CurrentElement % NodeIndexes
-          n = CurrentElement % TYPE % NumberOfNodes
+          Element => Mesh % Elements(i)
+          NodeIndexes =>  Element % NodeIndexes
+          n = Element % TYPE % NumberOfNodes
 
           Nodes % x(1:n) = Mesh % Nodes % x(NodeIndexes)
           Nodes % y(1:n) = Mesh % Nodes % y(NodeIndexes)
@@ -3060,10 +3060,10 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
           m = AdvVar % Perm(i+1) - AdvVar % Perm(i) 
           IF( m == 0 ) CYCLE
 
-          IP = GaussPoints(CurrentElement, m )
+          IP = GaussPoints(Element, m )
 
           DO j = 1, IP % n 
-            stat = ElementInfo( CurrentElement, Nodes, IP % v(j), IP % u(j), IP % w(j), detJ, Basis ) 
+            stat = ElementInfo( Element, Nodes, IP % v(j), IP % u(j), IP % w(j), detJ, Basis ) 
             No = AdvVar % Perm(i) + j
 
             Coord(1) = SUM( Basis(1:n) * Nodes % x(1:n) )
@@ -3102,13 +3102,13 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
         GotScale = ( ABS( DGScale - 1.0_dp ) > TINY( DgScale ) )
 
         DO i=1,NoElements                       
-          CurrentElement => Mesh % Elements(i)
-          NodeIndexes =>  CurrentElement % NodeIndexes
-          n = CurrentElement % TYPE % NumberOfNodes
+          Element => Mesh % Elements(i)
+          NodeIndexes =>  Element % NodeIndexes
+          n = Element % TYPE % NumberOfNodes
 
-          CurrentElement => Mesh % Elements(i)
-          NodeIndexes =>  CurrentElement % NodeIndexes
-          n = CurrentElement % TYPE % NumberOfNodes
+          Element => Mesh % Elements(i)
+          NodeIndexes =>  Element % NodeIndexes
+          n = Element % TYPE % NumberOfNodes
 
           IF( GotScale ) THEN
             Center(1) = SUM( Mesh % Nodes % x(NodeIndexes ) ) / n
@@ -3123,7 +3123,7 @@ SUBROUTINE InitializeParticles( Particles, InitParticles, AppendParticles, &
           END IF
 
           DO j = 1, n
-            No = AdvVar % Perm( CurrentElement % DgIndexes(j) )
+            No = AdvVar % Perm( Element % DgIndexes(j) )
             IF( No == 0 ) CYCLE
             k = NodeIndexes(j)
             Coord(1) = Mesh % Nodes % x(k)
@@ -4404,7 +4404,8 @@ END SUBROUTINE ParticleWallKernel
         AccurateNow = AccurateAlways
       END IF
       FaceIndex0 = 0
-
+      ElementIndex0 = 0
+
 200   ElementIndex = GetParticleElement( Particles, No )
       Rfin = GetParticleCoord( Particles, No )
       Velo = GetParticleVelo( Particles, No )      
@@ -4517,10 +4518,10 @@ END SUBROUTINE LocateParticles
   !> to a global coordinate. Sloppy tolerances are used since we *should* 
   !> have already located the element.
   !--------------------------------------------------------------------------
-  FUNCTION ParticleElementInfo( CurrentElement, GlobalCoord, &
+  FUNCTION ParticleElementInfo( Element, GlobalCoord, &
       SqrtElementMetric, Basis, dBasisdx ) RESULT ( stat )
 
-    TYPE(Element_t), POINTER :: CurrentElement
+    TYPE(Element_t), POINTER :: Element
     REAL(KIND=dp) :: GlobalCoord(:), SqrtElementMetric, LocalDistance
     REAL(KIND=dp) :: Basis(:)
     REAL(KIND=dp), OPTIONAL :: dBasisdx(:,:)
@@ -4535,11 +4536,20 @@ FUNCTION ParticleElementInfo( CurrentElement, GlobalCoord, &
     INTEGER :: n
 
     SAVE ElementNodes
+
+    IF(.NOT. ASSOCIATED(Element) ) THEN
+      CALL Fatal('ParticleElementInfo','Element not associated!')
+    END IF
+
+    IF(.NOT. ASSOCIATED(Element % TYPE) ) THEN
+      PRINT *,'Element:',Element % ElementIndex
+      CALL Fatal('ParticleElementInfo','Element % Type not associated!')
+    END IF
 
-    n = CurrentElement % TYPE % NumberOfNodes
-    CALL GetElementNodes(ElementNodes,CurrentElement)
+    n = Element % TYPE % NumberOfNodes
+    CALL GetElementNodes(ElementNodes,Element)
 
-    Stat = PointInElement( CurrentElement, ElementNodes, &
+    Stat = PointInElement( Element, ElementNodes, &
         GlobalCoord, LocalCoord, GlobalEps = -1.0_dp, LocalEps = 1.0e3_dp, &
 	LocalDistance = LocalDistance ) 
 
@@ -4554,7 +4564,7 @@ FUNCTION ParticleElementInfo( CurrentElement, GlobalCoord, &
         ELSE
           CALL Warn(Caller,'Distance from element higher than expected!')
         END IF
-        PRINT *,'LocalDistance:',LocalDistance,'Element:',CurrentElement % ElementIndex
+        PRINT *,'LocalDistance:',LocalDistance,'Element:',Element % ElementIndex
         PRINT *,'Nodes X:',ElementNodes % x(1:n) - GlobalCoord(1)
         PRINT *,'Nodes Y:',ElementNodes % y(1:n) - GlobalCoord(2)
         PRINT *,'Nodes Z:',ElementNodes % z(1:n) - GlobalCoord(3)
@@ -4566,7 +4576,7 @@ FUNCTION ParticleElementInfo( CurrentElement, GlobalCoord, &
     v = LocalCoord(2)
     w = LocalCoord(3)
 
-    stat = ElementInfo( CurrentElement, ElementNodes, U, V, W, SqrtElementMetric, &
+    stat = ElementInfo( Element, ElementNodes, U, V, W, SqrtElementMetric, &
         Basis, dBasisdx )
     IF(.NOT. Stat) Misses(2) = Misses(2) + 1
 
@@ -4580,10 +4590,10 @@ END FUNCTION ParticleElementInfo
   !> size of individual solvers it has been hard coded here. 
   !--------------------------------------------------------------------------
 
-  SUBROUTINE GetVectorFieldInMesh(Var, CurrentElement, Basis, Velo, dBasisdx, GradVelo )
+  SUBROUTINE GetVectorFieldInMesh(Var, Element, Basis, Velo, dBasisdx, GradVelo )
 
     TYPE(Variable_t), POINTER :: Var
-    TYPE(Element_t) :: CurrentElement
+    TYPE(Element_t) :: Element
     REAL(KIND=dp) :: Basis(:), Velo(:) 
     REAL(KIND=dp), OPTIONAL :: dBasisdx(:,:), GradVelo(:,:)
 
@@ -4623,11 +4633,11 @@ SUBROUTINE GetVectorFieldInMesh(Var, CurrentElement, Basis, Velo, dBasisdx, Grad
     END IF
 
 
-    n = CurrentElement % TYPE % NumberOfNodes    
+    n = Element % TYPE % NumberOfNodes    
     IF( Var % TYPE == Variable_on_nodes_on_elements ) THEN      
-      LocalPerm(1:n) = Var % Perm( CurrentElement % DGIndexes )
+      LocalPerm(1:n) = Var % Perm( Element % DGIndexes )
     ELSE
-      LocalPerm(1:n) = Var % Perm( CurrentElement % NodeIndexes )
+      LocalPerm(1:n) = Var % Perm( Element % NodeIndexes )
     END IF
     npos = COUNT ( LocalPerm(1:n) > 0 )
 
@@ -4697,10 +4707,10 @@ END SUBROUTINE GetVectorFieldInMesh
   !> The routine returns a potential and its gradient.
   !--------------------------------------------------------------------------
 
-  SUBROUTINE GetScalarFieldInMesh(Var, CurrentElement, Basis, Pot, dBasisdx, GradPot )
+  SUBROUTINE GetScalarFieldInMesh(Var, Element, Basis, Pot, dBasisdx, GradPot )
 
     TYPE(Variable_t), POINTER :: Var
-    TYPE(Element_t) :: CurrentElement
+    TYPE(Element_t) :: Element
     REAL(KIND=dp) :: Basis(:), Pot 
     REAL(KIND=dp), OPTIONAL :: dBasisdx(:,:), GradPot(:)
 
@@ -4728,14 +4738,14 @@ SUBROUTINE GetScalarFieldInMesh(Var, CurrentElement, Basis, Pot, dBasisdx, GradP
 
     IF(.NOT. ASSOCIATED( Var ) ) RETURN
 
-    n = CurrentElement % TYPE % NumberOfNodes
+    n = Element % TYPE % NumberOfNodes
     IF( ASSOCIATED( Var % Perm ) ) THEN
-      LocalPerm(1:n) = Var % Perm( CurrentElement % NodeIndexes )
+      LocalPerm(1:n) = Var % Perm( Element % NodeIndexes )
       IF( .NOT. ALL ( LocalPerm(1:n) > 0 )) RETURN
       LocalField(1:n) = Var % Values( LocalPerm(1:n) )
     ELSE
       ! Some variables do not have permutation, most importantly the node coordinates
-      LocalField(1:n) = Var % Values( CurrentElement % NodeIndexes )
+      LocalField(1:n) = Var % Values( Element % NodeIndexes )
     END IF
 
     Pot = SUM( Basis(1:n) * LocalField(1:n) )

fem/src/SOLVER.KEYWORDS

@@ -1209,7 +1209,7 @@ Solver:Logical:     'Save Axis 3'
 Solver:Logical:     'Save Axis'
 Solver:Logical:     'Save Boundaries Only'
 Solver:Logical:     'Save Bulk Only'
-solver:logical:     'save metainfo'
+solver:logical:     'skip metainfo'
 Solver:Logical:     'Save Elemental Fields'
 Solver:Logical:     'Save Nodal Fields'
 Solver:Logical:     'Save Halo Elements Only'

fem/src/modules/ParticleDynamics.F90

@@ -1250,7 +1250,7 @@ SUBROUTINE ParticleFieldInteraction(Particles,dtime,SetParticles,SetFields )
 
        ElementIndex = GetParticleElement( Particles, No )
 
-       IF( ElementIndex == 0 ) CYCLE
+       IF( ElementIndex < 1 ) CYCLE
 
        BulkElement => Mesh % Elements( ElementIndex )
        IF(.NOT. ASSOCIATED( BulkElement ) ) CYCLE

fem/src/modules/ResultOutputSolve/VtuOutputSolver.F90

@@ -355,8 +355,10 @@ SUBROUTINE VtuOutputSolver( Model,Solver,dt,TransientSimulation )
     BinaryOutput = .NOT. AsciiOutput
   END IF
 
-  SaveMetainfo = GetLogical( Params,'Save Metainfo',GotIt )
-
+  ! Do not save metainfo if we use the results of saving for consistency test.
+  SaveMetainfo = .NOT. ( GetLogical( Params,'Skip Metainfo',GotIt ) .OR. &
+      ListCheckPresent( Params,'Reference Values') )
+
   ParallelBase = GetLogical( Params,'Partition Numbering',GotIt )
 
   IF( BinaryOutput ) THEN