From 34e831c53a80db647a90ea00940014d85696e5d6 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Mon, 6 Nov 2023 17:00:49 -0800 Subject: [PATCH 1/9] Add function for identifying cells that are cut off by the wall --- src/supersize_profile.jl | 47 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index cb75739..8914049 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -8,7 +8,7 @@ using Interpolations: Interpolations using GGDUtils: GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary, project_prop_on_subset!, get_subset_centers -using PolygonOps: PolygonOps +using PolygonOps: PolygonOps, inpolygon using JSON: JSON export extrapolate_core @@ -671,6 +671,40 @@ function modify_mesh_ext_near_x!( end end +function identify_cutoff_cells( + dd::OMAS.dd, + r::Matrix{Float64}, + z::Matrix{Float64}, + pfr_transition::Int64, +) + npol, nlvl = size(r) + okay = zeros(Bool, (npol - 1, nlvl - 1)) + border = zeros(Bool, (npol - 1, nlvl - 1)) + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + poly = [[wall_r[i], wall_z[i]] for i ∈ 1:length(wall_r)] + poly = append!(poly, [[wall_r[1], wall_z[1]]]) + inpoly = [[inpolygon([r[i, j], z[i, j]], poly) for j ∈ 1:nlvl] for i ∈ 1:npol] + wl = Matrix(reduce(hcat, inpoly)') .!= 0 + + for i ∈ 2:npol + if (i - 1) != pfr_transition + for j ∈ 2:nlvl + within = [wl[i, j], wl[i-1, j], wl[i-1, j-1], wl[i, j-1]] + okay[i-1, j-1] = all(within) + border[i-1, j-1] = (minimum(within) == 0) & (maximum(within) == 1) + end + end + end + + fz = open("limiter.dat", "w") + for i ∈ 1:length(wall_r) + println(fz, wall_r[i], " ", wall_z[i]) + end + return okay, border +end + #!format off """ record_regular_mesh!() @@ -979,6 +1013,17 @@ function mesh_extension_sol!( # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 pfr_transition = argmax(abs.(diff(pfr))) + okay, border = identify_cutoff_cells(dd,mesh_r,mesh_z,pfr_transition) + fr = open("okay.dat", "w") + fz = open("border.dat", "w") + for i ∈ 1:(npol-1) + for j ∈ 1:(nlvl-1) + print(fr, okay[i, j], " ") + print(fz, border[i, j], " ") + end + println(fr, "") + println(fz, "") + end record_regular_mesh!(grid_ggd, space, mesh_r, mesh_z, pfr_transition) return mesh_r, mesh_z, pfr_transition end From 41cf35052c944be5f1b737e132f4f522be87fdf1 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 09:22:33 -0800 Subject: [PATCH 2/9] improve ID of border cells - Instead of only looking at vertices, which can be deceiving if the wall has a bend in the right place --- Project.toml | 2 ++ src/supersize_profile.jl | 27 ++++++++++++++++++++++----- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index a3e1890..1b5c828 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,10 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 8914049..2f1624f 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -10,6 +10,8 @@ using GGDUtils: project_prop_on_subset!, get_subset_centers using PolygonOps: PolygonOps, inpolygon using JSON: JSON +using LibGEOS: LibGEOS, readgeom, intersection +using GeoInterface: GeoInterface export extrapolate_core export fill_in_extrapolated_core_profile @@ -685,18 +687,33 @@ function identify_cutoff_cells( wall_z = limiter.unit[1].outline.z poly = [[wall_r[i], wall_z[i]] for i ∈ 1:length(wall_r)] poly = append!(poly, [[wall_r[1], wall_z[1]]]) - inpoly = [[inpolygon([r[i, j], z[i, j]], poly) for j ∈ 1:nlvl] for i ∈ 1:npol] - wl = Matrix(reduce(hcat, inpoly)') .!= 0 + wall_poly = LibGEOS.Polygon([poly]) + # inpoly = [[inpolygon([r[i, j], z[i, j]], poly) for j ∈ 1:nlvl] for i ∈ 1:npol] + # wl = Matrix(reduce(hcat, inpoly)') .!= 0 + tolerance = 1e-7 # m^2 for i ∈ 2:npol if (i - 1) != pfr_transition for j ∈ 2:nlvl - within = [wl[i, j], wl[i-1, j], wl[i-1, j-1], wl[i, j-1]] - okay[i-1, j-1] = all(within) - border[i-1, j-1] = (minimum(within) == 0) & (maximum(within) == 1) + cell_r = [r[i, j], r[i, j-1], r[i-1, j-1], r[i-1, j], r[i, j]] + cell_z = [z[i, j], z[i, j-1], z[i-1, j-1], z[i-1, j], z[i, j]] + cell = LibGEOS.Polygon([[[cell_r[k], cell_z[k]] for k in 1:length(cell_r)]]) + # within = [wl[i, j], wl[i-1, j], wl[i-1, j-1], wl[i, j-1]] + # okay[i-1, j-1] = all(within) + # border[i-1, j-1] = (minimum(within) == 0) & (maximum(within) == 1) + cell_area = LibGEOS.area(cell) + inter = LibGEOS.intersection(wall_poly, cell) + inter_area = LibGEOS.area(inter) + okay[i-1, j-1] = abs(inter_area - cell_area) < tolerance + border[i-1, j-1] = 0 < inter_area < (cell_area - tolerance) + # println("i=",i,",j=",j,"okay=",okay[i-1,j-1],",ca=",cell_area,",ia=",inter_area) + + # inter = GeoInterface.coordinates(inter)[1] + end end end + # println("area of wall polygon:", LibGEOS.area(wall_poly)) fz = open("limiter.dat", "w") for i ∈ 1:length(wall_r) From 5b85833a46ccbdf3e0703cff11d03e6ec67943e8 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 11:42:25 -0800 Subject: [PATCH 3/9] wip: mesh ext trimming --- src/supersize_profile.jl | 65 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 61 insertions(+), 4 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 2f1624f..542d6ec 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -715,10 +715,10 @@ function identify_cutoff_cells( end # println("area of wall polygon:", LibGEOS.area(wall_poly)) - fz = open("limiter.dat", "w") - for i ∈ 1:length(wall_r) - println(fz, wall_r[i], " ", wall_z[i]) - end + # fz = open("limiter.dat", "w") + # for i ∈ 1:length(wall_r) + # println(fz, wall_r[i], " ", wall_z[i]) + # end return okay, border end @@ -876,6 +876,63 @@ function record_regular_mesh!( end end +function trim_ext_mesh_with_wall( + dd::OMAS.dd; + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + poly = [[wall_r[i], wall_z[i]] for i ∈ 1:length(wall_r)] + poly = append!(poly, [[wall_r[1], wall_z[1]]]) + wall_poly = LibGEOS.Polygon([poly]) + tolerance = 1e-7 # m^2 + + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + + o0 = space.objects_per_dimension[1] # Nodes + o1 = space.objects_per_dimension[2] # Edges + o2 = space.objects_per_dimension[3] # Cells (2D projection) + + ext_nodes_sub = get_grid_subset_with_index(grid_ggd, -201) + ext_edges_sub = get_grid_subset_with_index(grid_ggd, -202) + ext_xedges_sub = get_grid_subset_with_index(grid_ggd, -203) + ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) + ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) + + for j ∈ 1:length(ext_cells_sub.element) + i = length(ext_cells_sub.element) - (j-1) + nodes = o2.object[ext_cells_sub.element[i].index].nodes + c = [[o0.object[n].geometry.r, o0.object[n].geometry.z] for n in nodes] + cell = LibGEOS.Polygon([c]) + cell_area = LibGEOS.area(cell) + inter = LibGEOS.intersection(wall_poly, cell) + inter_area = LibGEOS.area(inter) + if abs(inter_area - cell_area) > tolerance + # Problem! This cell needs modification. + if cell_area == 0 + # Easy: just delete this guy AND all his nodes + deleteat!(o2.object, ext_cells_sub.element[i].index) + for node in nodes + deleteat!(o1.object, node) + end + deleteat!(ext_cells_sub.element, i) + end + # inter = GeoInterface.coordinates(inter)[1] + end + + end + # println("area of wall polygon:", LibGEOS.area(wall_poly)) + + fz = open("limiter.dat", "w") + for i ∈ 1:length(wall_r) + println(fz, wall_r[i], " ", wall_z[i]) + end + return okay, border +end + """ convert_filename(filename::String) From 1ee716860c207bb7632befcc5957b110143d4eda Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 12:48:20 -0800 Subject: [PATCH 4/9] more mesh trimming --- src/supersize_profile.jl | 40 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 542d6ec..de3e084 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -902,6 +902,8 @@ function trim_ext_mesh_with_wall( ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) + deleted_nodes = [] + deleted_edges = [] for j ∈ 1:length(ext_cells_sub.element) i = length(ext_cells_sub.element) - (j-1) nodes = o2.object[ext_cells_sub.element[i].index].nodes @@ -914,15 +916,47 @@ function trim_ext_mesh_with_wall( # Problem! This cell needs modification. if cell_area == 0 # Easy: just delete this guy AND all his nodes + edges = o2.object[ext_cells_sub.element[i].index].boundary + edges = [edge.object[k].index for k in length(edge.object)] deleteat!(o2.object, ext_cells_sub.element[i].index) - for node in nodes - deleteat!(o1.object, node) + for edge in sort(edges, rev=true) + deleteat!(o1.object, edge) + deleted_edges = [deleted_edges; edge] + end + for node in sort(nodes, rev=true) + deleteat!(o0.object, node) + deleted_nodes = [deleted_nodes; node] end deleteat!(ext_cells_sub.element, i) end # inter = GeoInterface.coordinates(inter)[1] end - + end + sub_edges = [ext_edges_sub.element[i].object[1].index for i in 1:length(ext_edges_sub.element)] + sub_xedges = [ext_xedges_sub.element[i].object[1].index for i in 1:length(ext_xedges_sub.element)] + sub_yedges = [ext_yedges_sub.element[i].object[1].index for i in 1:length(ext_yedges_sub.element)] + for j in 1:length(sub_edges) + i = length(sub_edges) - (j-1) + if sub_edges[i] in deleted_edges + deleteat!(ext_edges_sub, i) + end + if i <= length(sub_xedges) + if sub_xedges[i] in deleted_edges + deleteat!(ext_xedges_sub, i) + end + end + if i <= length(sub_yedges) + if sub_yedges[i] in deleted_edges + deleteat!(ext_yedges_sub, i) + end + end + end + sub_nodes = [ext_nodes_sub.element[i].object[1].index for i in 1:length(ext_nodes_sub)] + for j in 1:length(sub_nodes) + i = length(sub_nodes) - (j-1) + if sub_nodes[i] in deleted_nodes + deleteat!(ext_nodes_sub, i) + end end # println("area of wall polygon:", LibGEOS.area(wall_poly)) From d2d3193dcc844f760365d2a68f60ebb6e8bfea7e Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 15:45:40 -0800 Subject: [PATCH 5/9] wip cell filtering stuff --- src/supersize_profile.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index de3e084..53733b7 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -928,8 +928,24 @@ function trim_ext_mesh_with_wall( deleted_nodes = [deleted_nodes; node] end deleteat!(ext_cells_sub.element, i) + else + # Hard: cell is partly in bounds + # Find which nodes are out of bounds + # Find which edges are out of bounds or intersecting + # Delete all nodes and edges that aren't entirely in bounds + # Delete the cell + # Determine whether new nodes were already added by a previous cell + # The last npol + 1 cells need to be considered based on the order + # in which cells were added. + # Add new nodes that will be unique + # Determine whether new edges were already added by a previous cell + # Add new edges that will be unique + # Create new cell + # Reference the new boundary and nodes + p = LibGEOS.Point() + inter = GeoInterface.coordinates(inter)[1] end - # inter = GeoInterface.coordinates(inter)[1] + end end sub_edges = [ext_edges_sub.element[i].object[1].index for i in 1:length(ext_edges_sub.element)] From 7c6a10e3bbfd63d8f52247858ed5a65fedc3b305 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 8 Nov 2023 12:37:56 -0800 Subject: [PATCH 6/9] Repair equilibrium (add rho if missing) while loading from geqdsk to dd --- src/SD4SOLPS.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 0f87fe4..ab8a12d 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -8,6 +8,7 @@ using Interpolations: Interpolations export find_files_in_allowed_folders export geqdsk_to_imas! +export preparation include("$(@__DIR__)/supersize_profile.jl") include("$(@__DIR__)/repair_eq.jl") @@ -168,6 +169,9 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim limiter.unit[1].outline.z = g.zlim + + # Repair + add_rho_to_equilibrium!(dd) return end From 2d9cd3ebc0599dd5fe0e6ce6e2f721cbc938145a Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 8 Nov 2023 13:15:47 -0800 Subject: [PATCH 7/9] Bugfix: fix PFR split in extended mesh --- src/supersize_profile.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 53733b7..eff7856 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -831,7 +831,7 @@ function record_regular_mesh!( add_subset_element!(all_nodes_sub, space_idx, node_dim, node_idx) # Edges - if (i > 1) & (i != cut) # i-1 to i in the npol direction + if (i > 1) & (i != (cut+1)) # i-1 to i in the npol direction edge_idx1 = edge_start1 + iii * e1_per_i + jj edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] resize!(edges[edge_idx1].boundary, 2) @@ -855,7 +855,7 @@ function record_regular_mesh!( end # Cells - if (i > 1) & (i != cut) & (j > 1) + if (i > 1) & (i != (cut+1)) & (j > 1) cell_idx = cell_start + iii * c_per_i + jjj cells[cell_idx].nodes = [ node_idx, From 16a0d0ec6811aba38e1a78301c468de8af957dd4 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 8 Nov 2023 15:11:39 -0800 Subject: [PATCH 8/9] some progress on blocked cell removal, but there are issues --- data/.gitignore | 2 - src/supersize_profile.jl | 142 +++++++++++++++++++++++++++------------ 2 files changed, 98 insertions(+), 46 deletions(-) delete mode 100644 data/.gitignore diff --git a/data/.gitignore b/data/.gitignore deleted file mode 100644 index f6bc585..0000000 --- a/data/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.mesh_ext.json - diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index eff7856..8c5fca7 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -876,7 +876,7 @@ function record_regular_mesh!( end end -function trim_ext_mesh_with_wall( +function trim_ext_mesh_with_wall!( dd::OMAS.dd; grid_ggd_idx::Int64=1, space_idx::Int64=1, @@ -902,33 +902,42 @@ function trim_ext_mesh_with_wall( ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) + all_nodes_sub = get_grid_subset_with_index(grid_ggd, 1) + all_edges_sub = get_grid_subset_with_index(grid_ggd, 2) + all_xedges_sub = get_grid_subset_with_index(grid_ggd, 3) + all_yedges_sub = get_grid_subset_with_index(grid_ggd, 4) + all_cells_sub = get_grid_subset_with_index(grid_ggd, 5) + deleted_nodes = [] deleted_edges = [] - for j ∈ 1:length(ext_cells_sub.element) - i = length(ext_cells_sub.element) - (j-1) - nodes = o2.object[ext_cells_sub.element[i].index].nodes - c = [[o0.object[n].geometry.r, o0.object[n].geometry.z] for n in nodes] + deleted_cells = [] + ext_cells = length(ext_cells_sub.element) + for j ∈ 1:ext_cells + i = ext_cells - (j-1) + cell_idx = ext_cells_sub.element[i].object[1].index + nodes = o2.object[cell_idx].nodes + c = [[o0.object[n].geometry[1], o0.object[n].geometry[2]] for n in nodes] + c = [c; [c[1]]] cell = LibGEOS.Polygon([c]) cell_area = LibGEOS.area(cell) inter = LibGEOS.intersection(wall_poly, cell) inter_area = LibGEOS.area(inter) if abs(inter_area - cell_area) > tolerance + # print("cell ", i, " has problem: inter=", inter_area, ",cell=",cell_area) # Problem! This cell needs modification. - if cell_area == 0 + if inter_area == 0 + # println(", inter area is 0") # Easy: just delete this guy AND all his nodes - edges = o2.object[ext_cells_sub.element[i].index].boundary - edges = [edge.object[k].index for k in length(edge.object)] - deleteat!(o2.object, ext_cells_sub.element[i].index) - for edge in sort(edges, rev=true) - deleteat!(o1.object, edge) - deleted_edges = [deleted_edges; edge] - end - for node in sort(nodes, rev=true) - deleteat!(o0.object, node) - deleted_nodes = [deleted_nodes; node] - end + edges = o2.object[ext_cells_sub.element[i].object[1].index].boundary + edges = [edge.index for edge in edges] + deleted_edges = [deleted_edges; edges] + deleted_nodes = [deleted_nodes; nodes] + deleteat!(ext_cells_sub.element, i) + # deleteat!(o2.object, cell_idx) + deleted_cells = [deleted_cells; i] else + # println("cell area is not 0") # Hard: cell is partly in bounds # Find which nodes are out of bounds # Find which edges are out of bounds or intersecting @@ -942,45 +951,89 @@ function trim_ext_mesh_with_wall( # Add new edges that will be unique # Create new cell # Reference the new boundary and nodes - p = LibGEOS.Point() + # p = LibGEOS.Point() inter = GeoInterface.coordinates(inter)[1] end end end - sub_edges = [ext_edges_sub.element[i].object[1].index for i in 1:length(ext_edges_sub.element)] - sub_xedges = [ext_xedges_sub.element[i].object[1].index for i in 1:length(ext_xedges_sub.element)] - sub_yedges = [ext_yedges_sub.element[i].object[1].index for i in 1:length(ext_yedges_sub.element)] - for j in 1:length(sub_edges) - i = length(sub_edges) - (j-1) - if sub_edges[i] in deleted_edges - deleteat!(ext_edges_sub, i) + println("# deleted edges = ", length(deleted_edges)) + sub_edges = [ele.object[1].index for ele in ext_edges_sub.element] + sub_xedges = [ele.object[1].index for ele in ext_xedges_sub.element] + sub_yedges = [ele.object[1].index for ele in ext_yedges_sub.element] + asub_edges = [ele.object[1].index for ele in all_edges_sub.element] + asub_xedges = [ele.object[1].index for ele in all_xedges_sub.element] + asub_yedges = [ele.object[1].index for ele in all_yedges_sub.element] + unique!(sort!(deleted_nodes, rev=true)) + unique!(sort!(deleted_edges, rev=true)) + nae = length(asub_edges) + for j in 1:nae + i = nae - (j-1) + if i <= length(sub_edges) + if sub_edges[i] in deleted_edges + deleteat!(ext_edges_sub.element, i) + end end if i <= length(sub_xedges) if sub_xedges[i] in deleted_edges - deleteat!(ext_xedges_sub, i) + deleteat!(ext_xedges_sub.element, i) end end if i <= length(sub_yedges) if sub_yedges[i] in deleted_edges - deleteat!(ext_yedges_sub, i) + deleteat!(ext_yedges_sub.element, i) + end + end + if asub_edges[i] in deleted_edges + deleteat!(all_edges_sub.element, i) + end + if i <= length(asub_xedges) + if asub_xedges[i] in deleted_edges + deleteat!(all_xedges_sub.element, i) + end + end + if i <= length(asub_yedges) + if asub_yedges[i] in deleted_edges + deleteat!(all_yedges_sub.element, i) + end + end + end + sub_nodes = [ele.object[1].index for ele in ext_nodes_sub.element] + asub_nodes = [ele.object[1].index for ele in all_nodes_sub.element] + an = length(asub_nodes) + for j in 1:an + i = an - (j-1) + if asub_nodes[i] in deleted_nodes + deleteat!(all_nodes_sub.element, i) + end + if i <= length(sub_nodes) + if sub_nodes[i] in deleted_nodes + deleteat!(ext_nodes_sub.element, i) end end end - sub_nodes = [ext_nodes_sub.element[i].object[1].index for i in 1:length(ext_nodes_sub)] - for j in 1:length(sub_nodes) - i = length(sub_nodes) - (j-1) - if sub_nodes[i] in deleted_nodes - deleteat!(ext_nodes_sub, i) + asub_cells = [ele.object[1].index for ele in all_cells_sub.element] + nac = length(asub_cells) + for j in 1:nac + i = nac - (j - 1) + if asub_cells[i] in deleted_cells + deleteat!(all_cells_sub.element, i) end end # println("area of wall polygon:", LibGEOS.area(wall_poly)) + # println("length(o0.object)=",length(o0.object), ", maximum(deleted_nodes)=",maximum(deleted_nodes)) + # for node in deleted_nodes + # deleteat!(o0.object, node) + # end + for edge in deleted_edges + deleteat!(o1.object, edge) + end fz = open("limiter.dat", "w") for i ∈ 1:length(wall_r) println(fz, wall_r[i], " ", wall_z[i]) end - return okay, border + return end """ @@ -1137,18 +1190,19 @@ function mesh_extension_sol!( # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 pfr_transition = argmax(abs.(diff(pfr))) - okay, border = identify_cutoff_cells(dd,mesh_r,mesh_z,pfr_transition) - fr = open("okay.dat", "w") - fz = open("border.dat", "w") - for i ∈ 1:(npol-1) - for j ∈ 1:(nlvl-1) - print(fr, okay[i, j], " ") - print(fz, border[i, j], " ") - end - println(fr, "") - println(fz, "") - end + # okay, border = identify_cutoff_cells(dd,mesh_r,mesh_z,pfr_transition) + # fr = open("okay.dat", "w") + # fz = open("border.dat", "w") + # for i ∈ 1:(npol-1) + # for j ∈ 1:(nlvl-1) + # print(fr, okay[i, j], " ") + # print(fz, border[i, j], " ") + # end + # println(fr, "") + # println(fz, "") + # end record_regular_mesh!(grid_ggd, space, mesh_r, mesh_z, pfr_transition) + trim_ext_mesh_with_wall!(dd; grid_ggd_idx=grid_ggd_idx, space_idx=space_idx) return mesh_r, mesh_z, pfr_transition end From 5fa825cb43233645487927a8debd79906a8d0f6a Mon Sep 17 00:00:00 2001 From: David Eldon Date: Fri, 2 Feb 2024 10:23:48 -0800 Subject: [PATCH 9/9] wip: working on mesh trimming --- src/supersize_profile.jl | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 8c5fca7..d8b3bfc 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -673,6 +673,21 @@ function modify_mesh_ext_near_x!( end end +""" + itentify_cutoff_cells() + +Finds cells in the extended mesh that are cut off by the wall. Some cells may +intersect the wall and be divided, and some may be outside the wall completely. +Returns two flags: One is "okay" for cells where the cell area matches the area +of the intersection between the cell and the tokamak interior (to within +some small tolerance), and "border" for cells where the area of intersection is +less than this but more than 0. +dd: the data dictionary +r: the R coordinates of the mesh cells to check +z: the Z coordinates of the mesh cells to check +pfr_transition: the index of the mesh row that should not be connected because + it spans the transition from the PFR to the common flux region of the SOL. +""" function identify_cutoff_cells( dd::OMAS.dd, r::Matrix{Float64}, @@ -907,6 +922,7 @@ function trim_ext_mesh_with_wall!( all_xedges_sub = get_grid_subset_with_index(grid_ggd, 3) all_yedges_sub = get_grid_subset_with_index(grid_ggd, 4) all_cells_sub = get_grid_subset_with_index(grid_ggd, 5) + all_yedge_indices = [element.object[1].index for element in all_yedges_sub.element] deleted_nodes = [] deleted_edges = [] @@ -925,11 +941,12 @@ function trim_ext_mesh_with_wall!( if abs(inter_area - cell_area) > tolerance # print("cell ", i, " has problem: inter=", inter_area, ",cell=",cell_area) # Problem! This cell needs modification. + edges = o2.object[ext_cells_sub.element[i].object[1].index].boundary + edges = [edge.index for edge in edges] + yedges = [index for index in edges if index in all_yedge_indices] if inter_area == 0 # println(", inter area is 0") # Easy: just delete this guy AND all his nodes - edges = o2.object[ext_cells_sub.element[i].object[1].index].boundary - edges = [edge.index for edge in edges] deleted_edges = [deleted_edges; edges] deleted_nodes = [deleted_nodes; nodes] @@ -952,7 +969,23 @@ function trim_ext_mesh_with_wall!( # Create new cell # Reference the new boundary and nodes # p = LibGEOS.Point() + + # -------- okay, maybe that would work, but wouldn't it be easier + # to find nodes that are out of bounds and move them along their + # edges until they are on the wall? Then some of the edges would + # be fixed and there would be no invalid nodes. A complicated wall + # shape might require more nodes and edges to be added. inter = GeoInterface.coordinates(inter)[1] + for yedge in yedges + print(yedge, ": ") + yedge_nodes = o1.object[yedge].nodes + print(yedge_nodes, ": ") + edge_node_coords = [LibGEOS.Point(o0.object[yedge_nodes[1]].geometry[1], o0.object[yedge_nodes[1]].geometry[2]), LibGEOS.Point(o0.object[yedge_nodes[2]].geometry[1], o0.object[yedge_nodes[2]].geometry[2])] + println(edge_node_coords) + edge_line = LibGEOS.LineString([edge_node_coords]) + intersection_point = LibGEOS.intersection(wall_poly, edge_line) + println(intersection_point) + end end end