diff --git a/src/SOLPS2IMAS.jl b/src/SOLPS2IMAS.jl index bb135a0..c618855 100644 --- a/src/SOLPS2IMAS.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 180d030..8f41840 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using YAML: load_file as YAML_load_file using ArgParse: ArgParse import OMAS as IMASDD +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