Recover mapvtk funcionality from old observation module

Author: adrianarce-elemwave

src_output/mapVTKOutput.F90

@@ -72,6 +72,7 @@ subroutine store_relevant_coordinates(this, problemInfo)
       this%nPoints = counter
       call alloc_and_init(this%coords, 3, this%nPoints, -99)
       call alloc_and_init(this%materialTag, this%nPoints, -1)
+      call alloc_and_init(this%currentType, this%nPoints, -1)
 
       counter = 0
       do k = this%mainCoords%Z, this%auxCoords%Z
@@ -81,13 +82,13 @@ subroutine store_relevant_coordinates(this, problemInfo)
             if (isWithinBounds(Hfield, i, j, k, problemInfo)) then
                if (isMaterialExceptPML(Hfield, i, j, k, problemInfo)) then
                   counter = counter + 1
-                  call writeFaceTagInfo(this, counter, i, j, k, problemInfo%materialTag%getFaceTag(Hfield, i, j, k))
+                  call writeFaceTagInfo(this, counter, i, j, k, Hfield, problemInfo%materialTag%getFaceTag(Hfield, i, j, k))
                end if
                if (problemInfo%materialTag%getFaceTag(Hfield, i, j, k) < 0 &
                    .and. (btest(iabs(problemInfo%materialTag%getFaceTag(Hfield, i, j, k)), Hfield - 1)) &
                    .and. .not. isPML(Hfield, i, j, k, problemInfo)) then
                   counter = counter + 1
-                  call writeFaceTagInfo(this, counter, i, j, k, problemInfo%materialTag%getFaceTag(Hfield, i, j, k))
+                  call writeFaceTagInfo(this, counter, i, j, k, Hfield, problemInfo%materialTag%getFaceTag(Hfield, i, j, k))
                end if
             end if
          end do
@@ -103,12 +104,13 @@ logical function isMaterialExceptPML(field, i, j, k, problemInfo)
       isMaterialExceptPML = isMaterialExceptPML .and. (.not. isPML(field, i, j, k, problemInfo))
    end function isMaterialExceptPML
 
-   subroutine writeFaceTagInfo(this, counter, i, j, k, tag)
+   subroutine writeFaceTagInfo(this, counter, i, j, k, field, tag)
       type(mapvtk_output_t), intent(inout) :: this
-      integer, intent(in) :: i, j, k, counter, tag
+      integer, intent(in) :: i, j, k, counter, tag, field
       this%coords(1, counter) = i
       this%coords(2, counter) = j
       this%coords(3, counter) = k
+      this%currentType(counter) = currentType(field)
       this%materialTag(counter) = tag
    end subroutine
 
@@ -135,29 +137,44 @@ subroutine create_geometry_simulation_vtu(this, control)
       vtuPath = join_path(this%path, get_last_component(this%path))//vtuFileExtension
 
       call createUnstructuredDataForVTU(this%nPoints, this%coords, this%currentType, nodes, edges, quads, numNodes, numEdges, numQuads)
-      call ugrid%add_points(nodes)
+      call ugrid%add_points(real(nodes, 4))
 
       allocate (conn(2*numEdges + 4*numQuads))
       conn(1:2*numEdges) = reshape(edges, [2*numEdges])
       conn(2*numEdges + 1:2*numEdges + 4*numQuads) = reshape(quads, [4*numQuads])
 
       allocate (offsets(numEdges + numQuads))
-      do i = 1, numEdges
-         if (i == 1) then
-            offsets(i) = 2
+      do i = 1, numEdges + numQuads
+         if (i <= numEdges) then
+            if (i == 1) then
+               offsets(i) = 2
+            else
+               offsets(i) = offsets(i-1) + 2
+            end if
          else
-            offsets(i) = offsets(i - 1) + 2
+            if (i == 1) then
+               offsets(i) = 4
+            else
+               offsets(i) = offsets(i-1) + 4
+            end if
          end if
       end do
-      do i = 1, numQuads
-         offsets(numEdges + i) = offsets(numEdges + i - 1) + 4
-      end do
 
       allocate(types(numEdges+numQuads))
       types(1:numEdges) = 3
       types(numEdges+1:numEdges+numQuads) = 9
 
       call ugrid%add_cell_connectivity(conn, offsets, types)
+      if (size(offsets) /= numQuads) then
+         print *, "Problema con offsets"
+      end if
+      if (size(types) /= numQuads) then
+         print *, "Problema con types"
+      end if
+      if (offsets(numQuads) /= size(conn)) then
+         print *, "Tenemos un problema con conn y offset"
+      end if
+
       call ugrid%write_file(vtuPath)
 
       !---------------------------------------------

src_output/output.F90

@@ -188,7 +188,7 @@ subroutine init_outputs(sgg, media, sinpml_fullsize, materialTags, bounds, contr
 
                allocate (outputs(outputCount)%mapvtkOutput)
                call init_solver_output(outputs(outputCount)%mapvtkOutput, lowerBound, upperBound, outputRequestType, outputTypeExtension, control%mpidir, problemInfo)
-               call create_geometry_simulation_vtk(outputs(outputCount)%mapvtkOutput, problemInfo, control)
+               call create_geometry_simulation_vtu(outputs(outputCount)%mapvtkOutput, control)
                !! call adjust_computation_range --- Required due to issues in mpi region edges
 
             case (iEx, iEy, iEz, iHx, iHy, iHz)

src_output/outputUtils.F90

@@ -31,6 +31,7 @@ module mod_outputUtils
    public :: computeJ2
    public :: fieldo
    public :: create_data_file
+   public :: currentType
    !===========================
 
    !===========================
@@ -575,6 +576,19 @@ integer function u(field1, field2)
       end if
    end function
 
+   integer function currentType(field)
+        integer(kind=4) :: field
+        select case (field)
+        case (iEx); currentType = iJx
+        case (iEy); currentType = iJy
+        case (iEz); currentType = iJz
+        case (iHx); currentType = iBloqueJx
+        case (iHy); currentType = iBloqueJy
+        case (iHz); currentType = iBloqueJz
+        case default; call StopOnError(0, 0, 'field is not a E or H field')
+        end select
+      end function
+
    function get_field(field, i, j, k, fields_reference) result(res)
       implicit none
       real(kind=rkind) :: res

src_output/volumicProbeUtils.F90

@@ -91,47 +91,43 @@ subroutine createUnstructuredDataForVTU(counter, coords, currentType, Nodes, Edg
       real(kind=RKIND), allocatable, dimension(:, :), intent(out) :: Nodes
       integer(kind=4), allocatable, dimension(:, :), intent(out) ::  Edges, Quads
 
-      if (counter /= 0) then
-         Allocate (Nodes(3, counter))
-      else
-         return
-      end if
+      if (counter == 0) return
 
       call countElements(counter, currentType, numEdges, numQuads)
 
       allocate (Edges(2, numEdges))
       allocate (Quads(4, numQuads))
-      allocate (Nodes(3, counter*(numEdges + numQuads)))
+      allocate (Nodes(3, 2*numEdges + 4*numQuads))
 
       call registerElements(counter, coords, currentType, Nodes, Edges, Quads)
-
+      return
    end subroutine
 
    subroutine registerNode(nodes, nodeIx, x, y, z)
       real(kind=RKIND), dimension(:, :), intent(inout) :: nodes
       integer(kind=SINGLE), intent(in) :: nodeIx, x, y, z
-
-      nodes(1, nodeIx) = x*1.0_RKIND
-      nodes(2, nodeIx) = y*1.0_RKIND
-      nodes(3, nodeIx) = z*1.0_RKIND
+      !We need to avoid using idx 0
+      nodes(1, nodeIx + 1) = x*1.0_RKIND
+      nodes(2, nodeIx + 1) = y*1.0_RKIND
+      nodes(3, nodeIx + 1) = z*1.0_RKIND
    end subroutine
 
    subroutine registerEdge(edges, edgeIdx, startNodeIdx, endNodeIdx)
       integer(kind=SINGLE), dimension(:, :), intent(inout) :: edges
       integer(kind=SINGLE), intent(in) :: edgeIdx, startNodeIdx, endNodeIdx
 
-      edges(1, edgeIdx) = startNodeIdx
-      edges(2, edgeIdx) = endNodeIdx
+      edges(1, edgeIdx + 1) = startNodeIdx
+      edges(2, edgeIdx + 1) = endNodeIdx
    end subroutine
 
    subroutine registerQuad(quads, quadIdx, firstNodeIdx, secondNodeIdx, thirdNodeIdx, fourthNodeIdx)
       integer(kind=SINGLE), dimension(:, :), intent(inout) :: quads
       integer(kind=SINGLE), intent(in) :: quadIdx, firstNodeIdx, secondNodeIdx, thirdNodeIdx, fourthNodeIdx
 
-      quads(1, quadIdx) = firstNodeIdx
-      quads(2, quadIdx) = secondNodeIdx
-      quads(2, quadIdx) = thirdNodeIdx
-      quads(2, quadIdx) = fourthNodeIdx
+      quads(1, quadIdx + 1) = firstNodeIdx
+      quads(2, quadIdx + 1) = secondNodeIdx
+      quads(3, quadIdx + 1) = thirdNodeIdx
+      quads(4, quadIdx + 1) = fourthNodeIdx
    end subroutine
 
    subroutine countElements(counter, currentType, numEdges, numQuads)
@@ -158,9 +154,9 @@ subroutine registerElements(counter, coords, currentType, Nodes, Edges, Quads)
       integer :: nodeIdx, quadIdx, edgeIdx
       integer :: i
 
-      nodeIdx = 0
-      quadIdx = 0
-      edgeIdx = 0
+      nodeIdx = -1
+      quadIdx = -1
+      edgeIdx = -1
 
       do i = 1, counter
          select case (currentType(i))