Tests for tri mesh and using resize! 4 nodes edges

Correct chosen_tri_edge_order has been identified from the sample
EIRENE files.

Added some verififcation tests for the EIRENE triangular mesh reading.
Still need to check some mismatch in common faces between EIRENE and
SOLPS in the PFR edge region.

Also now only as many ndoes and edges are being initialized as required
in the space. This gets rid of the extra initialization we have been
doing so far.

src/SOLPS2IMAS.jl

@@ -91,9 +91,9 @@ chosen_edge_order = [(1, (1, 2)),
     (4, (3, 1))]
 
 # Following convention is used to index the edges of a triangular cell
-chosen_tri_edge_order = [(1, (2, 3)),
-    (2, (3, 1)),
-    (3, (1, 2))]
+chosen_tri_edge_order = [(1, (1, 2)),
+    (2, (2, 3)),
+    (3, (1, 3))]
 
 """
     solps2imas(
@@ -206,19 +206,21 @@ function solps2imas(
 
             # Resizing objects to hold cell geometry data
             # Should be fewer than this many points, but this way we won't under-fill
-            nodes = resize!(o1.object, ncell * 4)  # Nodes (1D)
-            edges = resize!(o2.object, ncell * 4)  # Edges (2D)
+            # nodes = resize!(o1.object, ncell * 4)  # Nodes (1D)
+            # edges = resize!(o2.object, ncell * 4)  # Edges (2D)
+            nodes = o1.object                  # Nodes (1D)
+            edges = o2.object                  # Edges (2D)
             cells = resize!(o3.object, ncell)  # Cells (3D)
 
             # Initialize geometry for 1D objects(nodes), nodes for 2D objects(edges)
-            for i ∈ 1:(ncell*4)
-                nodes[i].geometry = [0.0, 0.0]
-                edges[i].nodes = [0, 0]
-                resize!(edges[i].boundary, 2)
-                for bnd ∈ edges[i].boundary
-                    bnd.neighbours = Int64[]
-                end
-            end
+            # for i ∈ 1:(ncell*4)
+            #     nodes[i].geometry = [0.0, 0.0]
+            #     edges[i].nodes = [0, 0]
+            #     resize!(edges[i].boundary, 2)
+            #     for bnd ∈ edges[i].boundary
+            #         bnd.neighbours = Int64[]
+            #     end
+            # end
             # Initialize nodes and boundaries for cells
             for i ∈ 1:(ncell)
                 cells[i].nodes = [0, 0, 0, 0]
@@ -248,6 +250,7 @@ function solps2imas(
                             cry[1, icorner, iy, ix],
                         )[1]
                         if i_existing == 0
+                            resize!(nodes, j)
                             nodes[j].geometry =
                                 [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]]
                             cells[ic].nodes[icorner] = j
@@ -266,7 +269,9 @@ function solps2imas(
                         edge_nodes = [cells[ic].nodes[icorner] for icorner ∈ edge_pair]
                         existing_edge_ind = search_edges(edges, edge_nodes)
                         if existing_edge_ind == 0
+                            resize!(edges, edge_ind)
                             edges[edge_ind].nodes = edge_nodes
+                            resize!(edges[edge_ind].boundary, 2)
                             for (ii, edge_bnd) ∈ enumerate(edges[edge_ind].boundary)
                                 edge_bnd.index = edge_nodes[ii]
                             end
@@ -513,7 +518,6 @@ function solps2imas(
         subset_cells = get_grid_subset_with_index(grid_ggd, 5)
         subset_b25nodes = grid_ggd.grid_subset[cur_no_subsets+1]
         subset_b25faces = grid_ggd.grid_subset[cur_no_subsets+2]
-        subset_b25cells = grid_ggd.grid_subset[cur_no_subsets+3]
         subset_trinodes = grid_ggd.grid_subset[cur_no_subsets+4]
         subset_trifaces = grid_ggd.grid_subset[cur_no_subsets+5]
         subset_tricells = grid_ggd.grid_subset[cur_no_subsets+6]
@@ -588,6 +592,9 @@ function solps2imas(
                     resize!(edges, length(edges) + 1)
                     this_edge_ind = length(edges)
                     edges[this_edge_ind].nodes = tri_edge_nodes
+                    for (ii, edge_bnd) ∈ enumerate(edges[this_edge_ind].boundary)
+                        edge_bnd.index = tri_edge_nodes[ii]
+                    end
                     edges[this_edge_ind].measure =
                         distance_between_nodes(nodes, tri_edge_nodes)
                     add_subset_element!(subset_faces, 1, 2, this_edge_ind)

test/runtests.jl

@@ -3,6 +3,7 @@ using Test
 using YAML: load_file as YAML_load_file
 using ArgParse: ArgParse
 using OMAS: OMAS
+import SOLPS2IMAS: get_grid_subset_with_index
 
 allowed_rtol = 1e-4
 
@@ -22,6 +23,9 @@ function parse_commandline()
         ["--parser"],
         Dict(:help => "Stress test file parsing (other than b2 output files)",
             :action => :store_true),
+        ["--fort"],
+        Dict(:help => "Test triangular mesh generation from fort files",
+            :action => :store_true),
     )
     args = ArgParse.parse_args(s)
     if !any(values(args)) # If no flags are set, run all tests
@@ -201,3 +205,63 @@ if args["solps2imas"]
         @test Set(brute_force_separatix_list) == Set(subset_separatix_element_list)
     end
 end
+
+if args["fort"]
+    @testset "Test triangular mesh generation from fort files" begin
+        fort = (
+            "$(@__DIR__)/../samples/fort.33",
+            "$(@__DIR__)/../samples/fort.34",
+            "$(@__DIR__)/../samples/fort.35")
+        b2gmtry = "$(@__DIR__)/../samples/b2fgmtry"
+        b2output = "$(@__DIR__)/../samples/b2time.nc"
+        b2mn = "$(@__DIR__)/../samples/b2mn.dat"
+        ids = SOLPS2IMAS.solps2imas(b2gmtry, b2output; b2mn=b2mn, fort=fort)
+        grid_ggd = ids.edge_profiles.grid_ggd[1]
+        space = grid_ggd.space[1]
+
+        subset_nodes = get_grid_subset_with_index(grid_ggd, 1)
+        subset_faces = get_grid_subset_with_index(grid_ggd, 2)
+        subset_cells = get_grid_subset_with_index(grid_ggd, 5)
+        subset_b25nodes = get_grid_subset_with_index(grid_ggd, -1)
+        subset_b25faces = get_grid_subset_with_index(grid_ggd, -2)
+        subset_b25cells = get_grid_subset_with_index(grid_ggd, -5)
+        subset_trinodes = get_grid_subset_with_index(grid_ggd, -101)
+        subset_trifaces = get_grid_subset_with_index(grid_ggd, -102)
+        subset_tricells = get_grid_subset_with_index(grid_ggd, -105)
+        subset_comnodes = get_grid_subset_with_index(grid_ggd, -201)
+        subset_comfaces = get_grid_subset_with_index(grid_ggd, -202)
+        subset_comcells = get_grid_subset_with_index(grid_ggd, -205)
+
+        nodes = grid_ggd.space[1].objects_per_dimension[1].object
+        edges = grid_ggd.space[1].objects_per_dimension[2].object
+        cells = grid_ggd.space[1].objects_per_dimension[3].object
+
+        gmtry = SOLPS2IMAS.read_b2_output(b2gmtry)
+        nx = gmtry["dim"]["nx"]
+        ny = gmtry["dim"]["ny"]
+        extra_nodes = 2 * (nx + ny)
+
+        @test length(subset_nodes.element) == length(nodes)
+        @test length(subset_faces.element) == length(edges)
+        @test length(subset_cells.element) == length(cells)
+
+        @test length(subset_b25nodes.element) ==
+              length(subset_comnodes.element) + extra_nodes
+
+        @test length(subset_trinodes.element) + extra_nodes == length(nodes)
+
+        for ele ∈ subset_tricells.element
+            tricell_ind = ele.object[1].index
+            for bnd_ind ∈ 1:3
+                if length(cells[tricell_ind].boundary[bnd_ind].neighbours) > 0
+                    common_nodes = intersect(
+                        cells[tricell_ind].nodes,
+                        cells[cells[tricell_ind].boundary[bnd_ind].neighbours[1]].nodes,
+                    )
+                    vert_no = Tuple(indexin(common_nodes, cells[tricell_ind].nodes))
+                    @test SOLPS2IMAS.chosen_tri_edge_order[bnd_ind][2] == vert_no
+                end
+            end
+        end
+    end
+end