diff --git a/Project.toml b/Project.toml index 76afd8c7..a87bfbf9 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [weakdeps] diff --git a/docs/Project.toml b/docs/Project.toml index b3e06ef2..31a744c4 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,7 +2,12 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" +Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" + +[compat] +Documenter = "1" +julia = "1.9" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index b497ffde..c78897c8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,40 +1,38 @@ -ENV["MPLBACKEND"]="agg" -using Documenter, ExtendableGrids, Literate, GridVisualize, SimplexGridFactory +ENV["MPLBACKEND"] = "agg" +using Documenter, ExtendableGrids, Literate, GridVisualize, SimplexGridFactory, Gmsh import CairoMakie, Triangulate -CairoMakie.activate!(type="svg",visible=false) +CairoMakie.activate!(; type = "svg", visible = false) -example_md_dir = joinpath(@__DIR__,"src","examples") +example_md_dir = joinpath(@__DIR__, "src", "examples") - - -examples1d=joinpath(@__DIR__,"..","examples","examples1d.jl") +examples1d = joinpath(@__DIR__, "..", "examples", "examples1d.jl") include(examples1d) -examples2d=joinpath(@__DIR__,"..","examples","examples2d.jl") +examples2d = joinpath(@__DIR__, "..", "examples", "examples2d.jl") include(examples2d) -examples3d=joinpath(@__DIR__,"..","examples","examples3d.jl") +examples3d = joinpath(@__DIR__, "..", "examples", "examples3d.jl") include(examples3d) include("makeplots.jl") function mkdocs() + Literate.markdown(examples1d, example_md_dir; documenter = false, info = false) + Literate.markdown(examples2d, example_md_dir; documenter = false, info = false) + Literate.markdown(examples3d, example_md_dir; documenter = false, info = false) + + generated_examples = joinpath.("examples", filter(x -> endswith(x, ".md"), readdir(example_md_dir))) - Literate.markdown(examples1d, example_md_dir, documenter=false,info=false) - Literate.markdown(examples2d, example_md_dir, documenter=false,info=false) - Literate.markdown(examples3d, example_md_dir, documenter=false,info=false) - - generated_examples=joinpath.("examples",filter(x->endswith(x, ".md"),readdir(example_md_dir))) - - makeplots(example_md_dir, Plotter=CairoMakie) - - makedocs(sitename="ExtendableGrids.jl", - modules = [ExtendableGrids], - clean = false, + makeplots(example_md_dir; Plotter = CairoMakie) + + makedocs(; sitename = "ExtendableGrids.jl", + modules = [ExtendableGrids, Base.get_extension(ExtendableGrids, :ExtendableGridsGmshExt)], + clean = false, + warnonly = true, doctest = true, authors = "J. Fuhrmann, Ch. Merdon", - repo="https://github.com/j-fu/ExtendableGrids.jl", - pages=[ - "Home"=>"index.md", + repo = "https://github.com/j-fu/ExtendableGrids.jl", + pages = [ + "Home" => "index.md", "adjacency.md", "vectorofconstants.md", "typehierarchy.md", @@ -52,12 +50,12 @@ function mkdocs() "refinement.md", "regionedit.md", "tokenstream.md", + "gmsh.md", "allindex.md", - "Examples" => generated_examples + "Examples" => generated_examples, ]) end mkdocs() -deploydocs(repo = "github.com/j-fu/ExtendableGrids.jl.git") - +deploydocs(; repo = "github.com/j-fu/ExtendableGrids.jl.git") diff --git a/docs/src/gmsh.md b/docs/src/gmsh.md new file mode 100644 index 00000000..c4ad6bb6 --- /dev/null +++ b/docs/src/gmsh.md @@ -0,0 +1,22 @@ +# Gmsh interoperability +This functionality is in beta stage. +Breaking changes for this API are considered non-breaking for the package. +Therefore, these functions are not exported yet. + + + +## API +These methods become available via a package extension which is loaded together with +[Gmsh.jl](https://github.com/JuliaFEM/Gmsh.jl). +See the general [gmsh documentation](https://gmsh.info/), the [Gmsh reference manual](https://gmsh.info/doc/texinfo/gmsh.html) +and the [Gmsh Julia API source code](https://gitlab.onelab.info/gmsh/gmsh/blob/master/api/gmsh.jl) for information. + + +```@docs +ExtendableGrids.simplexgrid_from_gmsh +ExtendableGrids.simplexgrid_to_gmsh +ExtendableGrids.mixedgrid_from_gmsh +ExtendableGrids.mixedgrid_to_gmsh +``` + + diff --git a/ext/ExtendableGridsGmshExt.jl b/ext/ExtendableGridsGmshExt.jl index 6996c91b..093c4375 100644 --- a/ext/ExtendableGridsGmshExt.jl +++ b/ext/ExtendableGridsGmshExt.jl @@ -6,42 +6,204 @@ else import ..Gmsh: gmsh end -import ExtendableGrids: ExtendableGrid, simplexgrid -import ExtendableGrids: Coordinates, CellNodes, CellRegions, BFaceNodes, BFaceRegions -import ExtendableGrids: simplexgrid_from_gmsh, write_gmsh +import ExtendableGrids: simplexgrid_from_gmsh, simplexgrid_to_gmsh +import ExtendableGrids: mixedgrid_from_gmsh, mixedgrid_to_gmsh +import ExtendableGrids: ExtendableGrid, simplexgrid, VariableTargetAdjacency, num_sources +import ExtendableGrids: Coordinates, CellNodes, CellRegions, BFaceNodes, BFaceRegions, CellGeometries, BFaceGeometries +import ExtendableGrids: Edge1D, Triangle2D, Quadrilateral2D, Tetrahedron3D, Hexahedron3D, Prism3D, VectorOfConstants, + ElementGeometries, Cartesian2D, Cartesian3D, CoordinateSystem, num_cells #!!! Make a license warning at initialization ? Gmsh is GPL - mention this in the readme. ###??? Do we really need this dependency here ? I would rather like to live without, the more that it seems to make some problems. +#using ExtendableGrids using StatsBase: countmap using Bijections +using UUIDs: uuid1 -###!!! I am still thinking about the architecture here. May be we should end up with grid_from_msh +""" +```` +simplexgrid_from_gmsh(filename::String; incomplete=false, Tc=Float32, Ti=Int32) +```` + +The msh file is read and a SimplexGrid is created. +The mesh can also contain an incomplete grid. For this, the function has to be called with ````incomplete=true````. +'incomplete' means that the grid only consists of nodes and cells, it does not have a boundary. We also do not try to read the physical groups for those grids. +`Tc` is the type of coordinates, `Ti` is the index type. + +""" +function simplexgrid_from_gmsh(filename::String; incomplete = false, Tc = Float32, Ti = Int32) + gmshfile_to_simplexgrid(filename; incomplete, Tc, Ti) +end + +""" +```` +mixedgrid_from_gmsh(filename::String; Tc=Float32, Ti=Int32) +```` + +The msh file is read and an ExtendableGrid is created. +This only works for dim=2 grids and the orientation may be wrong. +`Tc` is the type of coordinates, `Ti` is the index type. + +""" +function mixedgrid_from_gmsh(filename::String; Tc = Float32, Ti = Int32) + gmshfile_to_mixedgrid(filename, Tc, Ti) +end + +""" +```` +simplexgrid_from_gmsh(mod::Module; incomplete=false, Tc=Float32, Ti=Int32) +```` + +The mesh contained in the gmsh module is converted to a SimplexGrid. +The mesh can also contain an incomplete grid. For this, the function has to be called with ````incomplete=true````. +'incomplete' means that the grid only consists of nodes and cells, it does not have a boundary. We also do not try to read the physical groups for those grids. +`Tc` is the type of coordinates, `Ti` is the index type. + +""" +function simplexgrid_from_gmsh(mod::Module; incomplete = false, Tc = Float32, Ti = Int32) + if !incomplete + mod_to_simplexgrid(mod, Tc, Ti) + else + incomplete_mod_to_simplexgrid(mod, Tc, Ti) + end +end + +""" +```` +mixedgrid_from_gmsh(mod::Module; Tc=Float32, Ti=Int32) +```` + +The mesh contained in the gmsh module is converted to an ExtendableGrid. +`Tc` is the type of coordinates, `Ti` is the index type. + +""" +function mixedgrid_from_gmsh(mod::Module; Tc = Float32, Ti = Int32) + mod_to_mixedgrid(mod, Tc, Ti) +end + +""" +```` +simplexgrid_to_gmsh(g::ExtendableGrid; filename::String="") +```` + +The SimplexGrid 'g' is loaded into a gmsh module. +If a string (not "") is passed via 'filename', the mesh is written into this file. + + +""" +function simplexgrid_to_gmsh(g::ExtendableGrid; filename::String = "") + simplexgrid_to_gmshfile(g; filename) +end + +""" +```` +mixedgrid_to_gmsh(g::ExtendableGrid; filename::String="") +```` -# ExtendableGrids method extensions -function simplexgrid_from_gmsh(filename::String) - test_gmsh_init() - gmshfile_to_grid(filename) +The ExtendableGrid 'g' is loaded into a gmsh module. +If a string (not "") is passed via 'filename', the mesh is written into this file. + +""" +function mixedgrid_to_gmsh(g::ExtendableGrid; filename::String = "") + mixedgrid_to_gmshfile(g; filename) end +#------------------------------------------------------------------------------------------- -# ExtendableGrids method extensions -function simplexgrid_from_gmsh(mod::Module) - test_gmsh_init() - mod_to_grid(mod) +""" +```` +gmshfile_to_simplexgrid(filename::String, Tc, Ti) +```` + +This function just reads an msh file, and creates a gmsh.model and then calls the 'mod_to_simplexgrid' function +(This function has to be called with an initialized gmsh environment) +This function is called in 'simplexgrid_from_gmsh' +`Tc` is the type of coordinates, `Ti` is the index type. +""" +function gmshfile_to_simplexgrid(filename::String; incomplete = false, Tc, Ti) + gmsh.open(filename) + if !incomplete + return mod_to_simplexgrid(gmsh.model, Tc, Ti) + else + return incomplete_mod_to_simplexgrid(gmsh.model, Tc, Ti) + end end -###!!! maybe grid_to_msh, not sure yet -function write_gmsh(filename, g) - test_gmsh_init() - grid_to_gmshfile(g, filename) +""" +```` +gmshfile_to_mixedgrid(filename::String, Tc, Ti) +```` + +This function just reads an msh file, and creates a gmsh.model and then calls the 'mod_to_mixedgrid' function +(This function has to be called with an initialized gmsh environment) +This function is called in 'mixedgrid_from_gmsh' +`Tc` is the type of coordinates, `Ti` is the index type. +""" +function gmshfile_to_mixedgrid(filename::String, Tc, Ti) + gmsh.open(filename) + return mod_to_mixedgrid(gmsh.model, Tc, Ti) +end + +""" +```` +function simplexgrid_to_gmshfile(grid::ExtendableGrid, filename::String) +```` + +This function takes a simplexgrid, uses 'grid_to_mod' to create a corresponding gmsh module +Then it writes the module to a file +(this function has to be called with an initialized gmsh environment) +This function is called in 'write_gmsh' +""" +function simplexgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "") + if VERSION >= v"1.9" + # This possibility is new in 1.9, see + # https://github.com/JuliaLang/julia/blob/release-1.9/NEWS.md#new-language-features + gmsh.model = simplexgrid_to_mod(grid) + mod = gmsh.model + else + mod = simplexgrid_to_mod(grid) + end + + if filename != "" + gmsh.write(filename) + end + return mod end +""" +```` +mixedgrid_to_gmshfile(grid::ExtendableGrid, filename::String) +```` -###!!! all the rest is (modulo formatting) the untouched stuff from #29 +This function takes a mixed grid, uses 'grid_to_mod' to create a corresponding gmsh module +Then it writes the module to a file +(this function has to be called with an initialized gmsh environment) +This function is called in 'write_gmsh' + +grid[CellNodes] must be a VariableTargetAdjacency structure +""" +function mixedgrid_to_gmshfile(grid::ExtendableGrid; filename::String = "") + if VERSION >= v"1.9" + # This possibility is new in 1.9, see + # https://github.com/JuliaLang/julia/blob/release-1.9/NEWS.md#new-language-features + gmsh.model = mixedgrid_to_mod(grid) + mod = gmsh.model + else + mod = mixedgrid_to_mod(grid) + end + + if filename != "" + gmsh.write(filename) + end + return mod +end + +#--------------------------------------------------------------------------------------------- #= +#old names:!!!! this file contains the 2 main functions: a) "mod_to_grid": takes a gmsh.module and creates an ExtendableGrid from it b) "grid_to_mod": takes an ExtendableGrid and creates a gmsh.module from it @@ -50,43 +212,37 @@ for "mod_to_grid" there also exists "gmshfile_to_grid" which loads the content o for "grid_to_mod" there also exists "grid_to_gmshfile" which loads the content of a grid into a gmsh module and then calls "grid_to_mod" =# - ### general support functions +""" +```` +test_gmsh_init() +```` -function test_gmsh_init() +Very primitive function to test, via a try-catch-block, whether gmsh is already initialized. +If not, it will be initialized. +""" +function test_gmsh_init() try - #@warn "try to use gsmh" - #gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 1) - catch e - @warn "gmsh may not have been initialized. but is initialized now!" - gmsh.initialize() - end - #path = "/home/johannes/Nextcloud/Documents/Uni/VIII/WIAS/ExtendableGrids.jl_in_src/test/" - - #grid1 = ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh") - - #ExtendableGrids.grid_to_gmshfile(grid1, path*"testfile.msh") - - #grid2 = ExtendableGrids.simplexgrid_from_gmsh(path*"testfile.msh") - - #gmsh.finalize() - - #@test seemingly_equal2(grid1, grid2;confidence=:low) - #@test seemingly_equal2(grid1, grid2;confidence=:medium) + gmsh.option.setNumber("General.Terminal", 1) + catch e + @warn "gmsh may not have been initialized. but is initialized now!" + gmsh.initialize() + end end - - """ - take_second(x) - x is a list of 2-tuples, with an Int as second entry - an array of the second entries is returned +```` +take_second(x) +```` + +x is a list of 2-tuples, with an Int as second entry +an array of the second entries is returned """ function take_second(x) - y = zeros(Int64, length(x)) + a, b = x[1] + y = zeros(typeof(b), length(x)) for i = 1:length(x) _, t = x[i] y[i] = t @@ -94,45 +250,81 @@ function take_second(x) return y end - - """ - multiply_indices(indices, n) - for n=3: - [i, j, ..., k], 3 -> [3*i-2, 3*i-1, 3*i, 3*j-1, 3*j-2, 3*j, ..., 3*k-2, 3*k-1, 3*k] - in general: - [i, j, ..., k], n -> [n*i-(n-1), n*i-(n-2), ..., n*i, n*j-(n-1), ...] - This function can be used, if you have the indices of cells, and you want to get all their nodes, but the nodes are stored in one list for all cells: - [node_1_of_cell1, node_2_of_cell1, ... node_n_of_cell1, node_1_of_cell2, ...] +```` +multiply_indices(indices, n) +```` + +for n=3: +[i, j, ..., k], 3 -> [3*i-2, 3*i-1, 3*i, 3*j-1, 3*j-2, 3*j, ..., 3*k-2, 3*k-1, 3*k] +in general: +[i, j, ..., k], n -> [n*i-(n-1), n*i-(n-2), ..., n*i, n*j-(n-1), ...] +This function can be used, if you have the indices of cells, and you want to get all their nodes, but the nodes are stored in one list for all cells: +[node_1_of_cell1, node_2_of_cell1, ... node_n_of_cell1, node_1_of_cell2, ...] """ function multiply_indices(indices, n) m = length(indices) - ind_new = zeros(Int64, n * m) + ind_new = zeros(typeof(indices[1]), n * m) for i = 1:n - ind_new[(i-1)*m+1:i*m] = n * indices .- (n - i) + ind_new[((i - 1) * m + 1):(i * m)] = n * indices .- (n - i) end return sort(ind_new) end +""" +```` +use_vta(VTA, col_ids, num) +```` + +If VTA were a matrix, the result would be equivalent to VTA[:, col_ids]. +Each column of the VTA contains the nodes of one cell. +""" +function use_vta(VTA, col_ids, num) #note + result = zeros(Int64, num * length(col_ids)) + count = 1 + for j in col_ids + for i = 1:num + result[count] = VTA[i, j] + count += 1 + end + end + return result +end +""" +```` +use_geoms(cellgeoms, ids) +```` +If cellgeoms would just be an array/vector, the result would be equivalent to cellgeoms[ids]. +""" +function use_geoms(cellgeoms, ids) + res_cellgeoms = [] + for id in ids + push!(res_cellgeoms, cellgeoms[id]) + end + return res_cellgeoms +end """ - mod_to_grid(model) - function that tries to create an (simplex-) extendable grid from an gmsh.model - model has to be a gmsh.model - (this function has to be called with an initialized gmsh environment) +```` +mod_to_grid(model::Module, Tc, Ti) +```` -""" -function mod_to_grid(model::Module) +Function that tries to create an (simplex-) ExtendableGrid from a gmsh.model. +Model has to be a gmsh.model. +(This function has to be called with an initialized gmsh environment). +This function is called in 'simplexgrid_from_gmsh'. +`Tc` is the type of coordinates, `Ti` is the index type. +""" +function mod_to_simplexgrid(model::Module, Tc, Ti) dim = model.getDimension() node_names, coords, _ = model.mesh.getNodes() cell_types, element_names_cells, base_nodes_cells = model.mesh.getElements(dim, -1) face_types, element_names_faces, base_nodes_faces = model.mesh.getElements(dim - 1, -1) - #check whether cells are tetrahedrons in 3d or triangles in 2d: if dim == 3 @@ -155,9 +347,9 @@ function mod_to_grid(model::Module) if dim == 3 coords_new = reshape(coords, (3, Int(length(coords) / 3))) else - coords_new = zeros(Float64, 2, Int(length(coords) / 3)) + coords_new = zeros(Tc, 2, Int(length(coords) / 3)) for i = 1:Int(length(coords) / 3) - coords_new[:, i] = coords[3*i-2:3*i-1] + coords_new[:, i] = coords[(3 * i - 2):(3 * i - 1)] end end @@ -166,14 +358,13 @@ function mod_to_grid(model::Module) #the nodes making up the cells is stored in "base_nodes_cells", just in the wrong format simplices = reshape(base_nodes_cells[1], (dim + 1, m)) - #the physicalnames are currently unused - cellregion_to_physicalname = Bijection{Int64,String}() + cellregion_to_physicalname = Bijection{Ti, String}() pgnum_to_physcialname = Dict() cr_count = 1 #the cellregions correspond to the physical groups in which the cells are - cellregions = ones(Int64, m) + cellregions = ones(Ti, m) pgs = take_second(model.getPhysicalGroups(dim)) for pg in pgs @@ -186,7 +377,6 @@ function mod_to_grid(model::Module) cr_count += 1 end - for i = 1:m _, _, _, entitytag = model.mesh.getElement(element_names_cells[1][i]) for pg in pgs @@ -197,19 +387,17 @@ function mod_to_grid(model::Module) end end - - # assemble the boundary faces, just reads the faces stored in the msh file # for incomplete boundaries, there will be a function to seal them - + k = length(element_names_faces[1]) bfaces = reshape(base_nodes_faces[1], (dim, k)) # physical groups for bfaces - bfaceregions = ones(Int64, k) + bfaceregions = ones(Ti, k) pgs = take_second(model.getPhysicalGroups(dim - 1)) - bfaceregion_to_physicalname = Bijection{Int64,String}() + bfaceregion_to_physicalname = Bijection{Ti, String}() pgnum_to_physcialname = Dict() fr_count = 1 @@ -227,43 +415,192 @@ function mod_to_grid(model::Module) _, _, _, entitytag = model.mesh.getElement(element_names_faces[1][i]) for pg in pgs if entitytag in model.getEntitiesForPhysicalGroup(dim - 1, pg) - bfaceregions[i] = bfaceregion_to_physicalname(pgnum_to_physcialname[pg]) - break + bfaceregions[i] = bfaceregion_to_physicalname(pgnum_to_physcialname[pg]) + break end end end - - + return simplexgrid(coords_new, simplices, cellregions, bfaces, bfaceregions) - end """ - gmshfile_to_grid(filename::String) - this function just reads an msh file, and creates a gmsh.model and then calls the 'mod_to_grid' function - (this function has to be called with an initialized gmsh environment) +```` +incomplete_mod_to_simplexgrid(model::Module, Tc, Ti) +```` + +Loads an incomplete mesh from a msh file. +Then converts into an ExtendableGrids. +'incomplete' in this context means the boundary is missing. +With the 'ExtendableGrids.seal!(grid::ExtendableGrid)' the boundary can be added. +`Tc` is the type of coordinates, `Ti` is the index type. """ -function gmshfile_to_grid(filename::String) - - gmsh.open(filename) +function incomplete_mod_to_simplexgrid(model::Module, Tc, Ti) + dim = gmsh.model.getDimension() + + node_names, coords, _ = gmsh.model.mesh.getNodes() + cell_types, element_names_cells, base_nodes_cells = gmsh.model.mesh.getElements(dim, -1) + face_types, element_names_faces, base_nodes_faces = gmsh.model.mesh.getElements(dim - 1, -1) + + #check whether cells are tetrahedrons in 3d or triangles in 2d: + if dim == 3 + if cell_types[1] != 4 + @warn "3-dim file, but not (just) tetrahedrons as cells!!!" + return + end + elseif dim == 2 + if cell_types[1] != 2 + @warn "2-dim file, but not (just) triangles as cells!!!" + return + end + else + @warn "dim is neither 3 nor 2" + return + end + + #if dim=3, the coordinates (of the nodes) just have to be reshaped + #for dim=2, the z-coordinate has to be deleted + if dim == 3 + coords_new = reshape(coords, (3, Int(length(coords) / 3))) + else + coords_new = zeros(Tc, 2, Int(length(coords) / 3)) + for i = 1:Int(length(coords) / 3) + coords_new[:, i] = coords[(3 * i - 2):(3 * i - 1)] + end + end + + #number of cells + m = length(element_names_cells[1]) + #the nodes making up the cells is stored in "base_nodes_cells", just in the wrong format + simplices = reshape(base_nodes_cells[1], (dim + 1, m)) - grid = mod_to_grid(gmsh.model) + grid = ExtendableGrid{Tc, Ti}() + grid[Coordinates] = convert(Matrix{Tc}, coords_new) + grid[CellNodes] = convert(Matrix{Ti}, simplices) + grid[CellRegions] = ones(Ti, size(simplices)[2]) #parse(Ti, lines[elems1+1])) + + if dim == 2 + grid[CoordinateSystem] = Cartesian2D + grid[CellGeometries] = VectorOfConstants{ElementGeometries, Ti}(Triangle2D, num_cells(grid)) + else + grid[CoordinateSystem] = Cartesian3D + grid[CellGeometries] = VectorOfConstants{ElementGeometries, Ti}(Tetrahedron3D, num_cells(grid)) + end return grid end +""" +```` +mod_to_mixedgrid(model::Module, Tc, Ti) +```` + +Function that tries to create a (mixed-) ExtendableGrid from a gmsh.model. +Model has to be a gmsh.model. +(This function has to be called with an initialized gmsh environment). +This function is called in 'mixedgrid_from_gmsh'. +`Tc` is the type of coordinates, `Ti` is the index type. +""" +function mod_to_mixedgrid(model::Module, Tc, Ti) + elementtypes_nn_2d = Dict(2 => 3, 3 => 4) #key=elementtype id, val=num nodes + elementtypes_na_2d = Dict(2 => Triangle2D, 3 => Quadrilateral2D) + #elementtypes_na_2d = Dict(2=>ExtendableGrids.Triangle2D, 3=>ExtendableGrids.Quadrangle2D) + #elementtypes2d = [2, 3] + #numnodes_2d = [3, 4] + + elementtypes_nn_3d = Dict(4 => 4, 5 => 8, 6 => 6) #, 7=>5) #key=elementtype id, val=num nodes + elementtypes_na_3d = Dict(4 => Tetrahedron3D, 5 => Hexahedron3D, 6 => Prism3D) + #elementtypes_na_3d = Dict(4=>ExtendableGrids.Tetrahedron3D, 5=>ExtendableGrids.Hexahedron3D, 6=>ExtendableGrids.Prism3D) + #elementtypes3d = [4, 5, 6, 7] + #numnodes_3d = [4, 8, 6, 5] + + dim = model.getDimension() + + node_names, coords, _ = model.mesh.getNodes() + cell_types, element_names_cells, base_nodes_cells = model.mesh.getElements(dim, -1) + face_types, element_names_faces, base_nodes_faces = model.mesh.getElements(dim - 1, -1) + + grid = ExtendableGrid{Tc, Ti}() + + VTA = VariableTargetAdjacency(Ti) + #cellgeom = 0 #[] + cellgeom = VectorOfConstants{ElementGeometries, Ti}(Triangle2D, 0) + + if dim == 3 + @warn "dim=3 is not supported yet" + elseif dim == 2 + #@warn "dim=2" + #cells + for (ti, cell_type) in enumerate(cell_types) + #global cellgeom + + #@warn ti, cellgeom + nn = elementtypes_nn_2d[cell_type] + na = elementtypes_na_2d[cell_type] + + temp_cellgeom = VectorOfConstants{ElementGeometries, Ti}(na, length(element_names_cells[ti])) + + #if ti==1 + # cellgeom = temp_cellgeom + #else + cellgeom = vcat(cellgeom, temp_cellgeom) + #end + + #@warn nn, na + for (ci, cell) in enumerate(element_names_cells[ti]) + Base.append!(VTA, base_nodes_cells[ti][(nn * (ci - 1) + 1):(nn * ci)]) + #push!(cellgeom, na) + end + end + + coords_new = zeros(Tc, 2, Int(length(coords) / 3)) + for i = 1:Int(length(coords) / 3) + coords_new[:, i] = coords[(3 * i - 2):(3 * i - 1)] + end + grid[Coordinates] = convert(Matrix{Tc}, coords_new) + end + + #@warn cellgeom + + grid[CellGeometries] = cellgeom + grid[CellNodes] = VTA + grid[CellRegions] = ones(Ti, num_sources(VTA)) + + grid[BFaceGeometries] = VectorOfConstants{ElementGeometries, Ti}(Edge1D, length(element_names_faces[1])) + #@warn element_names_faces + + k = length(element_names_faces[1]) + grid[BFaceNodes] = convert(Matrix{Ti}, reshape(base_nodes_faces[1], (dim, k))) + grid[BFaceRegions] = ones(Ti, k) + + if dim == 2 + grid[CoordinateSystem] = Cartesian2D + #grid[CellGeometries] = VectorOfConstants{ElementGeometries,Int}(Triangle2D, num_cells(grid)) + else + grid[CoordinateSystem] = Cartesian3D + #grid[CellGeometries] = VectorOfConstants{ElementGeometries,Int}(Tetrahedron3D, num_cells(grid)) + end + + return grid +end """ - grid_to_mod(grid::ExtendableGrid) - this function writes an ExtendableGrid into a gmsh module - (this function has to be called with an initialized gmsh environment) +```` +grid_to_mod(grid::ExtendableGrid) +```` + +This function writes an ExtendableGrid into a gmsh module. +(This function has to be called with an initialized gmsh environment) +At the moment, this function can only be used from the outside via 'write_gmsh', where the newly created gmsh module is written into a msh file. """ -function grid_to_mod(grid::ExtendableGrid) +function simplexgrid_to_mod(grid::ExtendableGrid) # gmsh.initialize() gmsh.option.setNumber("General.Terminal", 1) #(fbase,fext)=splitext(filename) - gmsh.model.add("fbase") + gmsh.model.add("model" * string(uuid1())[1:8]) + Tc = typeof(grid[Coordinates][1, 1]) + Ti = typeof(grid[CellNodes][1, 1]) # formatting the coordinates correctly coords = grid[Coordinates] @@ -276,7 +613,7 @@ function grid_to_mod(grid::ExtendableGrid) elementtype_cell = 4 #tetrahedron elementtype_face = 2 #triangle else - coords3d = vcat(coords, zeros(Float64, (1, num_nodes))) + coords3d = vcat(coords, zeros(Tc, (1, num_nodes))) coords3d = reshape(coords3d, 3 * num_nodes) elementtype_cell = 2 #triangle elementtype_face = 1 #line @@ -288,12 +625,12 @@ function grid_to_mod(grid::ExtendableGrid) entitycount = 2 #Cells - cells = grid[CellNodes] - num_cells = size(cells, 2) - cellregions = grid[CellRegions] + cells = grid[CellNodes] + num_cells = size(cells, 2) + cellregions = grid[CellRegions] cm_cellregions = countmap(cellregions) - celltags0 = collect(num_nodes+1:num_nodes+num_cells) #there is a counter for all elements (nodes, cells, faces) added, since the nodes are already added, we have to start at num_nodes+1 - nodetags0 = reshape(cells, (dim + 1) * num_cells) #vector of nodes contained in the cells: [node_1_of_cell1, node_2_of_cell1, ... node_n_of_cell1, node_1_of_cell2, ...] + celltags0 = collect((num_nodes + 1):(num_nodes + num_cells)) #there is a counter for all elements (nodes, cells, faces) added, since the nodes are already added, we have to start at num_nodes+1 + nodetags0 = reshape(cells, (dim + 1) * num_cells) #vector of nodes contained in the cells: [node_1_of_cell1, node_2_of_cell1, ... node_n_of_cell1, node_1_of_cell2, ...] #we iterate over all cellregions. for each cellregion, we create an entity and add the cells of this cellregion to the entity. Then we add a PhysicalGroup containing the entity for cr_dict_entry in cm_cellregions @@ -311,14 +648,13 @@ function grid_to_mod(grid::ExtendableGrid) entitycount += 1 end - #Faces (basically the same as for the cells) - bfaces = grid[BFaceNodes] - num_bfaces = size(bfaces, 2) - bfaceregions = grid[BFaceRegions] + bfaces = grid[BFaceNodes] + num_bfaces = size(bfaces, 2) + bfaceregions = grid[BFaceRegions] cm_bfaceregions = countmap(bfaceregions) - facetags0 = collect(num_cells+num_nodes+1:num_cells+num_nodes+num_bfaces) - nodetags0 = reshape(bfaces, dim * num_bfaces) + facetags0 = collect((num_cells + num_nodes + 1):(num_cells + num_nodes + num_bfaces)) + nodetags0 = reshape(bfaces, dim * num_bfaces) for fr_dict_entry in cm_bfaceregions gmsh.model.addDiscreteEntity(dim - 1, entitycount) @@ -335,31 +671,154 @@ function grid_to_mod(grid::ExtendableGrid) entitycount += 1 end - return gmsh.model - end +function mixedgrid_to_mod(grid::ExtendableGrid) + elementtypes_nn_2d = Dict(2 => 3, 3 => 4) #key=elementtype id, val=num nodes + elementtypes_na_2d = Dict(Triangle2D => 2, Quadrilateral2D => 3) + #elementtypes_na_2d = Dict(2=>ExtendableGrids.Triangle2D, 3=>ExtendableGrids.Quadrangle2D) + #elementtypes2d = [2, 3] + #numnodes_2d = [3, 4] -""" - grid_to_gmshfile(grid::ExtendableGrid, filename::String) - this function takes a grid, uses 'grid_to_mod' to create a corresponding gmsh module - then it writes the module to a file - (this function has to be called with an initialized gmsh environment) -""" -function grid_to_gmshfile(grid::ExtendableGrid, filename::String) + elementtypes_nn_3d = Dict(4 => 4, 5 => 8, 6 => 6) #, 7=>5) #key=elementtype id, val=num nodes + elementtypes_na_3d = Dict(Tetrahedron3D => 4, Hexahedron3D => 5, Prism3D => 6) - if VERSION>=v"1.9" - # This possibility is new in 1.9, see - # https://github.com/JuliaLang/julia/blob/release-1.9/NEWS.md#new-language-features - gmsh.model = grid_to_mod(grid) + # gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 1) + #(fbase,fext)=splitext(filename) + gmsh.model.add("model" * string(uuid1())[1:8]) + + # formatting the coordinates correctly + coords = grid[Coordinates] + dim = size(coords, 1) + num_nodes = size(coords, 2) + nodetags = collect(1:num_nodes) + #types are taken from https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format-version-1-_0028Legacy_0029 + if dim == 3 + coords3d = reshape(coords, 3 * num_nodes) #vec(reshape(coords, (1,3*num_nodes))) + + elementtypes_nn_face = elementtypes_nn_2d #Dict(2=>3, 3=>4) + elementtypes_na_face = elementtypes_na_2d #Dict(Triangle2D=>2, Quadrilateral2D=>3) + + elementtypes_nn_cell = elementtypes_nn_3d #Dict(4=>4, 5=>8, 6=>6) + elementtypes_na_cell = elementtypes_na_3d #Dict(Tetrahedron3D=>4, Hexahedron3D=>5, Prism3D=>6) + + #elementtype_cell = 4 #tetrahedron + #elementtype_face = 2 #triangle else - grid_to_mod(grid) + coords3d = vcat(coords, zeros(Float64, (1, num_nodes))) + coords3d = reshape(coords3d, 3 * num_nodes) + + elementtypes_nn_face = Dict(1 => 2) + elementtypes_na_face = Dict(Edge1D => 1) + + elementtypes_nn_cell = elementtypes_nn_2d #Dict(2=>3, 3=>4) + elementtypes_na_cell = elementtypes_na_2d #Dict(Triangle2D=>2, Quadrilateral2D=>3) + + #elementtype_cell = 2 #triangle + elementtype_face = 1 #line end - gmsh.write(filename) - -end + #all nodes are added (to the model and to an entity of dimension 'dim' with tag 1) + gmsh.model.addDiscreteEntity(dim, 1) + gmsh.model.mesh.addNodes(dim, 1, nodetags, coords3d, []) + entitycount = 2 + #@warn grid[CellGeometries] + #Cells + cellgeoms = grid[CellGeometries] + VTA = grid[CellNodes] + cellregions = grid[CellRegions] + cm_cellregions = countmap(cellregions) + num_cells = 0 + #@warn cellgeoms + #@warn cellgeoms[1] + #@warn cellgeoms[1,2,3] + + #we iterate over all cellregions. for each cellregion, we create an entity and add the cells of this cellregion to the entity. Then we add a PhysicalGroup containing the entity + for cr_dict_entry in cm_cellregions + gmsh.model.addDiscreteEntity(dim, entitycount) + cr, num_elements = cr_dict_entry + array_with_element_ids = findall(z -> z == cr, cellregions) + #@warn array_with_element_ids + types_in_region = use_geoms(cellgeoms, array_with_element_ids) # cellgeoms[array_with_element_ids] + cm_types = countmap(types_in_region) + + for (type_EG_name, num) in cm_types + type_GMSH_name = elementtypes_na_cell[type_EG_name] + #@warn type_EG_name, type_GMSH_name + array_el_ids_in_reg_of_type = findall(z -> z == type_EG_name, use_geoms(cellgeoms, array_with_element_ids)) + + #@warn "done" + #@warn array_el_ids_in_reg_of_type + celltags = array_with_element_ids[array_el_ids_in_reg_of_type] #.+num_nodes + + #@warn celltags + #@warn VTA[celltags] + + nodetags = use_vta(VTA, celltags, elementtypes_nn_cell[type_GMSH_name]) #vcat(VTA[celltags]...) + num_cells += length(celltags) + celltags .+= num_nodes + + gmsh.model.mesh.addElementsByType(entitycount, type_GMSH_name, celltags, nodetags) + end + + #gmsh.model.mesh.addElementsByType(entitycount, elementtype_cell, celltags, nodetags) + gmsh.model.addPhysicalGroup(dim, [entitycount], cr) + + entitycount += 1 + end + + #---------------------------------------------------- + #= + cells = grid[CellNodes] + num_cells = size(cells, 2) + cellregions = grid[CellRegions] + cm_cellregions = countmap(cellregions) + celltags0 = collect(num_nodes+1:num_nodes+num_cells) #there is a counter for all elements (nodes, cells, faces) added, since the nodes are already added, we have to start at num_nodes+1 + nodetags0 = reshape(cells, (dim + 1) * num_cells) #vector of nodes contained in the cells: [node_1_of_cell1, node_2_of_cell1, ... node_n_of_cell1, node_1_of_cell2, ...] + + #we iterate over all cellregions. for each cellregion, we create an entity and add the cells of this cellregion to the entity. Then we add a PhysicalGroup containing the entity + for cr_dict_entry in cm_cellregions + gmsh.model.addDiscreteEntity(dim, entitycount) + cr, num_elements = cr_dict_entry + array_with_element_ids = findall(z -> z == cr, cellregions) + + #only select those cells which have the right cellregionnumber + celltags = celltags0[array_with_element_ids] + nodetags = nodetags0[multiply_indices(array_with_element_ids, dim + 1)] + + gmsh.model.mesh.addElementsByType(entitycount, elementtype_cell, celltags, nodetags) + gmsh.model.addPhysicalGroup(dim, [entitycount], cr) + + entitycount += 1 + end + =# + + #Faces (basically the same as for the cells) + bfaces = grid[BFaceNodes] + num_bfaces = size(bfaces, 2) + bfaceregions = grid[BFaceRegions] + cm_bfaceregions = countmap(bfaceregions) + facetags0 = collect((num_cells + num_nodes + 1):(num_cells + num_nodes + num_bfaces)) + nodetags0 = reshape(bfaces, dim * num_bfaces) + + for fr_dict_entry in cm_bfaceregions + gmsh.model.addDiscreteEntity(dim - 1, entitycount) + fr, num_elements = fr_dict_entry + array_with_element_ids = findall(z -> z == fr, bfaceregions) + + #only select those cells which have the right bfaceregionnumber + facetags = facetags0[array_with_element_ids] + nodetags = nodetags0[multiply_indices(array_with_element_ids, dim)] + gmsh.model.mesh.addElementsByType(entitycount, elementtype_face, facetags, nodetags) + gmsh.model.addPhysicalGroup(dim - 1, [entitycount], fr) + + entitycount += 1 + end + + return gmsh.model +end end diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl index 94027a83..81075e01 100644 --- a/src/ExtendableGrids.jl +++ b/src/ExtendableGrids.jl @@ -13,16 +13,14 @@ using Random using Dates using LinearAlgebra - -if !isdefined(Base, :get_extension) +if !isdefined(Base, :get_extension) using Requires end - include("adjacency.jl") -export Adjacency,VariableTargetAdjacency,FixedTargetAdjacency -export atranspose,num_targets,num_sources,num_links,append!, max_num_targets_per_source -export asparse,tryfix,makevar +export Adjacency, VariableTargetAdjacency, FixedTargetAdjacency +export atranspose, num_targets, num_sources, num_links, append!, max_num_targets_per_source +export asparse, tryfix, makevar include("serialadjacency.jl") export SerialVariableTargetAdjacency @@ -45,50 +43,45 @@ export AbstractElementGeometry1D export Edge1D export AbstractElementGeometry2D -export Polygon2D,Triangle2D,Quadrilateral2D,Pentagon2D,Hexagon2D,Parallelogram2D,Rectangle2D,Circle2D +export Polygon2D, Triangle2D, Quadrilateral2D, Pentagon2D, Hexagon2D, Parallelogram2D, Rectangle2D, Circle2D export AbstractElementGeometry3D -export Polyhedron3D,Tetrahedron3D, Hexahedron3D,Parallelepiped3D,RectangularCuboid3D,Prism3D,TrianglePrism3D,Sphere3D +export Polyhedron3D, Tetrahedron3D, Hexahedron3D, Parallelepiped3D, RectangularCuboid3D, Prism3D, TrianglePrism3D, Sphere3D export AbstractElementGeometry4D -export Polychoron4D,HyperCube4D +export Polychoron4D, HyperCube4D export dim_element - include("coordinatesystem.jl") -export coordinatesystems,CoordinateSystems +export coordinatesystems, CoordinateSystems export AbstractCoordinateSystem -export Cartesian1D,Cartesian2D,Cartesian3D -export Cylindrical2D,Cylindrical3D -export Polar2D,Polar1D ,Spherical3D,Spherical1D - - +export Cartesian1D, Cartesian2D, Cartesian3D +export Cylindrical2D, Cylindrical3D +export Polar2D, Polar1D, Spherical3D, Spherical1D include("extendablegrid.jl") export ExtendableGrid export instantiate, veryform export AbstractGridComponent -export AbstractGridAdjacency,AbstractElementGeometries,AbstractElementRegions +export AbstractGridAdjacency, AbstractElementGeometries, AbstractElementRegions export Coordinates, CellNodes, BFaceNodes export CellGeometries, BFaceGeometries export CellRegions, BFaceRegions, BEdgeRegions export NumCellRegions, NumBFaceRegions, NumBEdgeRegions export CoordinateSystem -export AbstractGridFloatArray1D,AbstractGridFloatArray2D -export AbstractGridIntegerArray1D,AbstractGridIntegerArray2D +export AbstractGridFloatArray1D, AbstractGridFloatArray2D +export AbstractGridIntegerArray1D, AbstractGridIntegerArray2D export index_type, coord_type export dim_space, dim_grid -export num_nodes, num_cells, num_bfaces, num_bedges +export num_nodes, num_cells, num_bfaces, num_bedges export num_cellregions, num_bfaceregions, num_bedgeregions export gridcomponents -export seemingly_equal +export seemingly_equal include("subgrid.jl") export subgrid - - include("shape_specs.jl") export refcoords_for_geometry export num_nodes @@ -103,7 +96,6 @@ export Tangent4ElemType! export xrefFACE2xrefCELL export xrefFACE2xrefOFACE - include("derived.jl") export Coordinates export CellVolumes, CellFaces, CellEdges, CellFaceSigns, CellFaceOrientations, CellEdgeSigns @@ -126,7 +118,7 @@ include("more.jl") export BFaceCells, BFaceNormals, BFaceEdges, BEdgeNodes include("voronoi.jl") -export tricircumcenter!,VoronoiFaceCenters +export tricircumcenter!, VoronoiFaceCenters include("meshrefinements.jl") export split_grid_into @@ -138,9 +130,7 @@ include("adaptive_meshrefinements.jl") export bulk_mark export RGB_refine - - -include("assemblytypes.jl"); +include("assemblytypes.jl") export AssemblyType export AT_NODES, ON_CELLS, ON_FACES, ON_IFACES, ON_BFACES, ON_EDGES, ON_BEDGES export ItemType4AssemblyType @@ -151,13 +141,12 @@ export GridComponentRegions4AssemblyType export GridComponentUniqueGeometries4AssemblyType export GridComponentAssemblyGroups4AssemblyType -include("l2gtransformations.jl"); +include("l2gtransformations.jl") export L2GTransformer, update_trafo!, eval_trafo!, mapderiv! include("cellfinder.jl") export CellFinder -export gFindLocal!, gFindBruteForce!, interpolate,interpolate! - +export gFindLocal!, gFindBruteForce!, interpolate, interpolate! include("commongrids.jl") export reference_domain @@ -167,13 +156,11 @@ export grid_unitsquare, grid_unitsquare_mixedgeometries export grid_triangle export ringsector - include("regionedit.jl") export cellmask!, bfacemask!, bedgemask!, rect! - include("arraytools.jl") -export glue,geomspace,linspace +export glue, geomspace, linspace include("simplexgrid.jl") export simplexgrid, geomspace, glue @@ -181,20 +168,19 @@ export XCoordinates, YCoordinates, ZCoordinates export writefile include("tokenstream.jl") -export TokenStream, gettoken, expecttoken,trytoken +export TokenStream, gettoken, expecttoken, trytoken include("io.jl") export writeVTK +include("seal.jl") -@static if !isdefined(Base, :get_extension) +@static if !isdefined(Base, :get_extension) function __init__() - @require Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" begin + @require Gmsh="705231aa-382f-11e9-3f0c-b7cb4346fdeb" begin include("../ext/ExtendableGridsGmshExt.jl") end end end - - end # module diff --git a/src/extendablegrid.jl b/src/extendablegrid.jl index 861e0320..42373db7 100644 --- a/src/extendablegrid.jl +++ b/src/extendablegrid.jl @@ -12,9 +12,9 @@ $(TYPEDEF) Grid type wrapping Dict """ -mutable struct ExtendableGrid{Tc,Ti} - components::Dict{Type{<:AbstractGridComponent},Any} - ExtendableGrid{Tc,Ti}() where{Tc,Ti} =new(Dict{Type{<:AbstractGridComponent},Any}()) +mutable struct ExtendableGrid{Tc, Ti} + components::Dict{Type{<:AbstractGridComponent}, Any} + ExtendableGrid{Tc, Ti}() where {Tc, Ti} = new(Dict{Type{<:AbstractGridComponent}, Any}()) end """ @@ -30,8 +30,7 @@ are created on demand. Due to the fact that components are stored as Any the return value triggers type instability. To prevent this, specialized methods must be (and are) defined. """ -Base.getindex( grid::ExtendableGrid, T::Type{<:AbstractGridComponent} ) = get!(grid,T) - +Base.getindex(grid::ExtendableGrid, T::Type{<:AbstractGridComponent}) = get!(grid, T) """ ```` @@ -42,8 +41,7 @@ Union type for element information arrays. If all elements have the same information, it can be stored in an economical form as a [`VectorOfConstants`](@ref). """ -const ElementInfo{T}=Union{Vector{T},VectorOfConstants{T}} - +const ElementInfo{T} = Union{Vector{T}, VectorOfConstants{T}} """ $(TYPEDEF) @@ -51,8 +49,9 @@ $(TYPEDEF) """ abstract type AbstractGridFloatArray1D <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridFloatArray1D} )::Array{Tc,1} where{Tc,Ti} get!(grid,T) end - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridFloatArray1D})::Array{Tc, 1} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -60,7 +59,9 @@ $(TYPEDEF) """ abstract type AbstractGridFloatArray2D <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridFloatArray2D} )::Array{Tc,2} where{Tc,Ti} get!(grid,T) end +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridFloatArray2D})::Array{Tc, 2} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -68,9 +69,9 @@ $(TYPEDEF) """ abstract type AbstractGridIntegerArray1D <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridIntegerArray1D} )::Array{Ti,1} where{Tc,Ti} get!(grid,T) end - - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridIntegerArray1D})::Array{Ti, 1} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -78,8 +79,9 @@ $(TYPEDEF) """ abstract type AbstractGridIntegerArray2D <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridIntegerArray2D} )::Array{Ti,1} where{Tc,Ti} get!(grid,T) end - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridIntegerArray2D})::Array{Ti, 1} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -87,8 +89,9 @@ Integer number """ abstract type AbstractGridIntegerConstant <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridIntegerConstant})::Ti where{Tc,Ti} get!(grid,T) end - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridIntegerConstant})::Ti where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -96,7 +99,9 @@ Floating point number """ abstract type AbstractGridFloatConstant <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridFloatConstant} )::Tc where{Tc,Ti} get!(grid,T) end +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridFloatConstant})::Tc where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -105,7 +110,9 @@ Any kind of adjacency between grid components """ abstract type AbstractGridAdjacency <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractGridAdjacency} )::Adjacency{Ti} where{Tc,Ti} get!(grid,T) end +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractGridAdjacency})::Adjacency{Ti} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -114,10 +121,10 @@ Array of element geometry information. """ abstract type AbstractElementGeometries <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractElementGeometries} )::ElementInfo{ElementGeometries} where{Tc,Ti} get!(grid,T) end - - - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, + T::Type{<:AbstractElementGeometries})::ElementInfo{ElementGeometries} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -126,8 +133,9 @@ Array of element region number information. """ abstract type AbstractElementRegions <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{<:AbstractElementRegions} )::ElementInfo{Ti} where{Tc,Ti} get!(grid,T) end - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{<:AbstractElementRegions})::ElementInfo{Ti} where {Tc, Ti} + get!(grid, T) +end """ $(TYPEDEF) @@ -136,9 +144,9 @@ Coordinate system """ abstract type CoordinateSystem <: AbstractGridComponent end -function Base.getindex( grid::ExtendableGrid{Tc,Ti}, T::Type{CoordinateSystem} )::CoordinateSystems where{Tc,Ti} get!(grid,T) end - - +function Base.getindex(grid::ExtendableGrid{Tc, Ti}, T::Type{CoordinateSystem})::CoordinateSystems where {Tc, Ti} + get!(grid, T) +end ################################################################## # Component key types: for these, access is type stable @@ -164,7 +172,6 @@ Adjacency describing nodes per grid boundary face """ abstract type BFaceNodes <: AbstractGridAdjacency end - """ $(TYPEDEF) @@ -221,8 +228,6 @@ Number of boundary edge regions """ abstract type NumBEdgeRegions <: AbstractGridIntegerConstant end - - ############################################################ # Generic component methods @@ -235,24 +240,21 @@ and possibly transform components to be added to the grid via `setindex!`. The default method just passes data through. """ -veryform(grid::ExtendableGrid,v,::Type{<:AbstractGridComponent})=v - +veryform(grid::ExtendableGrid, v, ::Type{<:AbstractGridComponent}) = v """ $(TYPEDSIGNATURES) Set new grid component """ -Base.setindex!(grid::ExtendableGrid,v,T::Type{<:AbstractGridComponent})= grid.components[T]=veryform(grid,v,T) +Base.setindex!(grid::ExtendableGrid, v, T::Type{<:AbstractGridComponent}) = grid.components[T] = veryform(grid, v, T) """ $(TYPEDSIGNATURES) Remove grid component """ -Base.delete!(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})= delete!(grid.components,T) - - +Base.delete!(grid::ExtendableGrid, T::Type{<:AbstractGridComponent}) = delete!(grid.components, T) """ $(TYPEDSIGNATURES) @@ -260,8 +262,8 @@ $(TYPEDSIGNATURES) To be called by getindex. This triggers lazy creation of non-existing gridcomponents """ -Base.get!(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})= get!( ()->veryform(grid,instantiate(grid,T),T), grid.components,T) - +Base.get!(grid::ExtendableGrid, T::Type{<:AbstractGridComponent}) = get!(() -> veryform(grid, instantiate(grid, T), T), + grid.components, T) """ $(TYPEDSIGNATURES) @@ -270,9 +272,6 @@ $(TYPEDSIGNATURES) """ function instantiate end - - - ############################################################ # Verfication of inserted data @@ -283,14 +282,8 @@ veryform(grid::ExtendableGrid{Tc,Ti},v,T::Type{<:AbstractGridAdjacency}) where{T Check proper type of adjacencies upon insertion """ -veryform(grid::ExtendableGrid{Tc,Ti},v,T::Type{<:AbstractGridAdjacency}) where{Tc,Ti}= typeof(v)<:Adjacency{Ti} ? v : throw("Type mismatch") - - - - - - - +veryform(grid::ExtendableGrid{Tc, Ti}, v, T::Type{<:AbstractGridAdjacency}) where {Tc, Ti} = typeof(v) <: Adjacency{Ti} ? v : + throw("Type mismatch") ############################################################ # Instantiation methods @@ -300,32 +293,29 @@ $(TYPEDSIGNATURES) Instantiate number of cell regions """ -instantiate(grid, ::Type{NumCellRegions})=maximum(grid[CellRegions]) +instantiate(grid, ::Type{NumCellRegions}) = maximum(grid[CellRegions]) """ $(TYPEDSIGNATURES) Instantiate number of bface regions """ -instantiate(grid, ::Type{NumBFaceRegions})=maximum(grid[BFaceRegions]) - +instantiate(grid, ::Type{NumBFaceRegions}) = maximum(grid[BFaceRegions]) """ $(TYPEDSIGNATURES) Instantiate number of boundary edge regions """ -instantiate(grid, ::Type{NumBEdgeRegions})=maximum(grid[BEdgeRegions]) +instantiate(grid, ::Type{NumBEdgeRegions}) = maximum(grid[BEdgeRegions]) function prepare_bedgeregions!(grid::ExtendableGrid) - bedges = grid[BEdgeNodes] - bedgeregions = zeros(Int32, num_bedges(grid)) + bedges = grid[BEdgeNodes] + bedgeregions = zeros(Int32, num_bedges(grid)) grid[BEdgeRegions] = bedgeregions end - -instantiate(grid, ::Type{BEdgeRegions})=prepare_bedgeregions!(grid) - +instantiate(grid, ::Type{BEdgeRegions}) = prepare_bedgeregions!(grid) ############################################################# # General methods @@ -335,112 +325,97 @@ $(TYPEDSIGNATURES) Print the hierarchy of grid component key types (subtypes of [`AbstractGridComponent`](@ref). This includes additionally user defined subptypes. """ -gridcomponents()=AbstractTrees.print_tree(AbstractGridComponent) - - +gridcomponents() = AbstractTrees.print_tree(AbstractGridComponent) """ $(TYPEDSIGNATURES) Keys in grid """ -Base.keys(g::ExtendableGrid)=Base.keys(g.components) - +Base.keys(g::ExtendableGrid) = Base.keys(g.components) """ $(TYPEDSIGNATURES) Check if key is in grid """ -Base.haskey(g::ExtendableGrid,k) = Base.haskey(g.components,k) - +Base.haskey(g::ExtendableGrid, k) = Base.haskey(g.components, k) """ $(SIGNATURES) Space dimension of grid """ -dim_space(grid::ExtendableGrid)= size(grid[Coordinates],1) - +dim_space(grid::ExtendableGrid) = size(grid[Coordinates], 1) """ $(SIGNATURES) Grid dimension dimension of grid (larges element dimension) """ -dim_grid(grid::ExtendableGrid)= dim_element(grid[CellGeometries][1]) - +dim_grid(grid::ExtendableGrid) = dim_element(grid[CellGeometries][1]) """ $(SIGNATURES) Number of nodes in grid """ -num_nodes(grid::ExtendableGrid)::Int= size(grid[Coordinates],2) - +num_nodes(grid::ExtendableGrid)::Int = size(grid[Coordinates], 2) """ $(TYPEDSIGNATURES) Number of cells in grid """ -num_cells(grid::ExtendableGrid)= num_sources(grid[CellNodes]) - +num_cells(grid::ExtendableGrid) = num_sources(grid[CellNodes]) """ $(TYPEDSIGNATURES) Number of boundary faces in grid. """ -num_bfaces(grid::ExtendableGrid)= haskey(grid,BFaceNodes) ? num_sources(grid[BFaceNodes]) : 0 +num_bfaces(grid::ExtendableGrid) = haskey(grid, BFaceNodes) ? num_sources(grid[BFaceNodes]) : 0 """ $(TYPEDSIGNATURES) Number of boundary edges in grid. """ -num_bedges(grid::ExtendableGrid)= haskey(grid,BEdgeNodes) ? num_sources(grid[BEdgeNodes]) : 0 - - +num_bedges(grid::ExtendableGrid) = haskey(grid, BEdgeNodes) ? num_sources(grid[BEdgeNodes]) : 0 """ $(TYPEDSIGNATURES) Maximum cell region number """ -num_cellregions(grid::ExtendableGrid)=grid[NumCellRegions] - +num_cellregions(grid::ExtendableGrid) = grid[NumCellRegions] """ $(TYPEDSIGNATURES) Maximum boundary face region numbers """ -num_bfaceregions(grid::ExtendableGrid)=grid[NumBFaceRegions] - +num_bfaceregions(grid::ExtendableGrid) = grid[NumBFaceRegions] """ $(TYPEDSIGNATURES) Maximum boundary edge region numbers """ -num_bedgeregions(grid::ExtendableGrid)=grid[NumBEdgeRegions] - - - +num_bedgeregions(grid::ExtendableGrid) = grid[NumBEdgeRegions] """ $(SIGNATURES) Type of coordinates in grid """ -coord_type(grid::ExtendableGrid{Tc, Ti}) where {Tc,Ti}=Tc +coord_type(grid::ExtendableGrid{Tc, Ti}) where {Tc, Ti} = Tc """ $(SIGNATURES) Type of indices """ -index_type(grid::ExtendableGrid{Tc, Ti}) where {Tc,Ti}=Ti +index_type(grid::ExtendableGrid{Tc, Ti}) where {Tc, Ti} = Ti """ $(TYPEDSIGNATURES) @@ -448,45 +423,116 @@ $(TYPEDSIGNATURES) Map a function onto node coordinates of grid """ function Base.map(f::Function, grid::ExtendableGrid) - coord=grid[Coordinates] - dim=dim_space(grid) - if dim==1 - @views Base.map(f,coord[1,:]) - elseif dim==2 - @views Base.map(f,coord[1,:], coord[2,:]) + coord = grid[Coordinates] + dim = dim_space(grid) + if dim == 1 + @views Base.map(f, coord[1, :]) + elseif dim == 2 + @views Base.map(f, coord[1, :], coord[2, :]) else - @views Base.map(f,coord[1,:], coord[2,:], coord[3,:]) + @views Base.map(f, coord[1, :], coord[2, :], coord[3, :]) end end - -function Base.show(io::IO,grid::ExtendableGrid) - if num_edges(grid)>0 - str=@sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d, edges: %d\n", - typeof(grid),dim_space(grid),num_nodes(grid), num_cells(grid), num_bfaces(grid), num_edges(grid)) +function Base.show(io::IO, grid::ExtendableGrid) + if num_edges(grid) > 0 + str = @sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d, edges: %d\n", + typeof(grid), dim_space(grid), num_nodes(grid), num_cells(grid), num_bfaces(grid), num_edges(grid)) else - str=@sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d\n", - typeof(grid),dim_space(grid),num_nodes(grid), num_cells(grid), num_bfaces(grid)) + str = @sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d\n", + typeof(grid), dim_space(grid), num_nodes(grid), num_cells(grid), num_bfaces(grid)) end - println(io,str) + println(io, str) +end + +### Tests for the gmsh extension: + +function multidimsort(A) + i1 = sortperm(A[1, :]) + A1 = A[:, i1] + for j = 2:size(A, 1) + cm = countmap(A1[j - 1, :]) + for (key, val) in cm + if val > 1 #if there only is one entry with this key, the reordering is not necessary + inds = findall(z -> z == key, A1[j - 1, :]) #indices of that key in A1 + sorted_inds = sortperm(A1[j, inds]) + A1[:, inds] = A1[:, inds[sorted_inds]] #reorder + end + end + end + return A1 end """ -$(SIGNATURES) + seemingly_equal(grid1, grid2; sort=false, confidence=:full Recursively check seeming equality of two grids. Seemingly means -that long arrays are only compared via random samples -""" -function seemingly_equal(grid1::ExtendableGrid, grid2::ExtendableGrid) - for key in keys(grid1) - if !haskey(grid2,key) - return false +that long arrays are only compared via random samples. + +Keyword args: + +- `sort`: if true, sort grid points +- `confidence`: Confidence level: + - :low : Point numbers etc are the same + - :full : all arrays are equal (besides the coordinate array, the arrays only have to be equal up to permutations) +""" +function seemingly_equal(grid1::ExtendableGrid, grid2::ExtendableGrid; sort = false, confidence = :full) + if !sort + for key in keys(grid1) + if !haskey(grid2, key) + return false + end + if !seemingly_equal(grid1[key], grid2[key]) + return false + end end - if !seemingly_equal(grid1[key],grid2[key]) - return false + return true + end + + if confidence == :full + for key in keys(grid1) + if !haskey(grid2, key) + return false + end + if !isa(grid1[key], AbstractArray) + !seemingly_equal(grid1[key], grid2[key]) && return false + continue + end + + s1 = size(grid1[key]) + s2 = size(grid2[key]) + + if length(s1) == 0 + if length(s1) == 1 + if eltype(grid1[key]) <: Number + ind1 = sortperm(grid1[key]) + ind2 = sortperm(grid2[key]) + sa1 = grid1[key][ind1] + sa2 = grid2[key][ind2] + !seemingly_equal(sa1, sa2) && return false + else + !seemingly_equal(grid1[key], grid2[key]) && return false + end + else + if eltype(grid1[key]) <: Number && key != Coordinates + sa1 = multidimsort(sort(grid1[key]; dims = 1)) + sa2 = multidimsort(sort(grid2[key]; dims = 1)) + !seemingly_equal(sa1, sa2) && return false + else + !seemingly_equal(grid1[key], grid2[key]) && return false + end + end + end end + return true + elseif confidence == :low + grid1_data = (num_nodes(grid1), num_cells(grid1), num_bfaces(grid1)) + grid2_data = (num_nodes(grid2), num_cells(grid2), num_bfaces(grid2)) + return grid1_data == grid2_data + else + error("Confidence level $(confidence) not implemented") + return false end - return true end """ @@ -495,29 +541,30 @@ $(SIGNATURES) Check for seeming equality of two arrays by random sample. """ function seemingly_equal(array1::AbstractArray, array2::AbstractArray) - if size(array1)!=size(array2) + if size(array1) != size(array2) return false end - l=length(array1) - ntests=Float64(min(l,50)) - p=min(0.5,ntests/l) - for i in randsubseq( 1:l,p) - if !seemingly_equal(array1[i],array2[i]) + l = length(array1) + ntests = Float64(min(l, 50)) + p = min(0.5, ntests / l) + for i in randsubseq(1:l, p) + if !seemingly_equal(array1[i], array2[i]) return false end end true end -seemingly_equal(a1::VariableTargetAdjacency, a2::VariableTargetAdjacency)=seemingly_equal(a1.colentries,a2.colentries)&& seemingly_equal(a1.colstart, a2.colstart) -seemingly_equal(x1::Type, x2::Type)=(x1==x2) -seemingly_equal(x1::Number, x2::Number)=(x1≈x2) -seemingly_equal(x1::Any, x2::Any)=(x1==x2) - +function seemingly_equal(a1::VariableTargetAdjacency, a2::VariableTargetAdjacency) + seemingly_equal(a1.colentries, a2.colentries) && seemingly_equal(a1.colstart, a2.colstart) +end +seemingly_equal(x1::Type, x2::Type) = (x1 == x2) +seemingly_equal(x1::Number, x2::Number) = (x1 ≈ x2) +seemingly_equal(x1::Any, x2::Any) = (x1 == x2) -Base.extrema(grid::ExtendableGrid)=Base.extrema(grid[Coordinates],dims=2) +Base.extrema(grid::ExtendableGrid) = Base.extrema(grid[Coordinates]; dims = 2) function bbox(grid) - e=extrema(grid) - map(a->a[1],e),map(a->a[2],e) + e = extrema(grid) + map(a -> a[1], e), map(a -> a[2], e) end diff --git a/src/io.jl b/src/io.jl index 914edbc9..c12485bf 100644 --- a/src/io.jl +++ b/src/io.jl @@ -2,10 +2,10 @@ using WriteVTK # conversion from AbstractElementGeometry to WriteVTK.VTKCellTypes VTKCellType(::Type{<:AbstractElementGeometry1D}) = VTKCellTypes.VTK_LINE -VTKCellType(::Type{<:Triangle2D}) = VTKCellTypes.VTK_TRIANGLE -VTKCellType(::Type{<:Quadrilateral2D}) = VTKCellTypes.VTK_QUAD -VTKCellType(::Type{<:Tetrahedron3D}) = VTKCellTypes.VTK_TETRA -VTKCellType(::Type{<:Hexahedron3D}) = VTKCellTypes.VTK_HEXAHEDRON +VTKCellType(::Type{<:Triangle2D}) = VTKCellTypes.VTK_TRIANGLE +VTKCellType(::Type{<:Quadrilateral2D}) = VTKCellTypes.VTK_QUAD +VTKCellType(::Type{<:Tetrahedron3D}) = VTKCellTypes.VTK_TETRA +VTKCellType(::Type{<:Hexahedron3D}) = VTKCellTypes.VTK_HEXAHEDRON """ $(TYPEDSIGNATURES) @@ -18,25 +18,181 @@ exports grid and optional provided data as a vtk file Each '(key, value)' pair adds another data entry to the vtk file via WriteVTK functionality. """ -function writeVTK(filename::String, grid::ExtendableGrid{Tc,Ti}; kwargs...) where {Tc,Ti} - - ncells = num_cells(grid) # get number of cells in grid - coords = grid[Coordinates] # get coordinates - cells = grid[CellNodes] # get cell-node list +function writeVTK(filename::String, grid::ExtendableGrid{Tc, Ti}; kwargs...) where {Tc, Ti} + ncells = num_cells(grid) # get number of cells in grid + coords = grid[Coordinates] # get coordinates + cells = grid[CellNodes] # get cell-node list geo_dim = size(coords, 1) - cell_geo = grid[CellGeometries] + cell_geo = grid[CellGeometries] - vtk_cells = Array{MeshCell,1}(undef, ncells) + vtk_cells = Array{MeshCell, 1}(undef, ncells) - for icell in 1:ncells + for icell = 1:ncells vtk_cells[icell] = MeshCell(VTKCellType(cell_geo[icell]), view(cells, :, icell)) - end + end vtk_grid(filename, coords, vtk_cells) do vtk for (key, value) in kwargs vtk[String(key)] = value end end +end + +""" +$(TYPEDSIGNATURES) + +Write grid to file. Currently for pdelib sg format only +""" +function Base.write(fname::String, g::ExtendableGrid; format = "") + (fbase, fext) = splitext(fname) + if format == "" + format = fext[2:end] + end + @assert format == "sg" + + dim_g = dim_grid(g) + dim_s = dim_space(g) + nn = num_nodes(g) + nc = num_cells(g) + nbf = num_bfaces(g) + coord = g[Coordinates] + cellnodes = g[CellNodes] + bfacenodes = g[BFaceNodes] + cellregions = g[CellRegions] + bfaceregions = g[BFaceRegions] + + # TODO: replace @sprintf by someting non-allocating + open(fname, "w") do file + write(file, @sprintf("SimplexGrid")) + write(file, @sprintf(" ")) + write(file, @sprintf("2.1\n")) + write(file, @sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n")) + write(file, @sprintf("#mailto:{fuhrmann|streckenbach}@wias-berlin.de\n")) + write(file, @sprintf("#%s\n", Dates.format(Dates.now(), "yyyy-mm-ddTHH-mm-SS"))) + + write(file, @sprintf("DIMENSION\n%d\n", dim_g)) + write(file, @sprintf("NODES\n%d %d\n", nn, dim_s)) + + for inode = 1:nn + for idim = 1:dim_s + write(file, @sprintf("%.20e ", coord[idim, inode])) + write(file, @sprintf("\n")) + end + end + + write(file, @sprintf("CELLS\n%d\n", nc)) + for icell = 1:nc + for inode = 1:(dim_g + 1) + write(file, @sprintf("%d ", cellnodes[inode, icell])) + end + write(file, @sprintf("%d\n", cellregions[icell])) + end + + write(file, @sprintf("FACES\n%d\n", nbf)) + for ibface = 1:nbf + for inode = 1:dim_g + write(file, @sprintf("%d ", bfacenodes[inode, ibface])) + end + write(file, @sprintf("%d\n", bfaceregions[ibface])) + end + write(file, @sprintf("END\n")) + flush(file) + flush(file) + end + nothing +end + +###################################################### +""" +$(TYPEDSIGNATURES) + +Read grid from file. Currently for pdelib sg format only. +""" +function simplexgrid(file::String; format = "") + Ti = Cint + (fbase, fext) = splitext(file) + if format == "" + format = fext[2:end] + end + @assert format == "sg" + + tks = TokenStream(file) + expecttoken(tks, "SimplexGrid") + version = parse(Float64, gettoken(tks)) + version20 = false + + if (version == 2.0) + version20 = true + elseif (version == 2.1) + version20 = false + else + error("Read grid: wrong format version: $(version)") + end + dim::Ti = 0 + coord = Array{Float64, 2}(undef, 0, 0) + cells = Array{Ti, 2}(undef, 0, 0) + regions = Array{Ti, 1}(undef, 0) + faces = Array{Ti, 2}(undef, 0, 0) + bregions = Array{Ti, 1}(undef, 0) + while (true) + if (trytoken(tks, "DIMENSION")) + dim = parse(Ti, gettoken(tks)) + elseif (trytoken(tks, "NODES")) + nnodes = parse(Ti, gettoken(tks)) + embdim = parse(Ti, gettoken(tks)) + if (dim != embdim) + error("Dimension error (DIMENSION $(dim)) in section NODES") + end + coord = Array{Float64, 2}(undef, dim, nnodes) + for inode = 1:nnodes + for idim = 1:embdim + coord[idim, inode] = parse(Float64, gettoken(tks)) + end + end + elseif (trytoken(tks, "CELLS")) + ncells = parse(Ti, gettoken(tks)) + cells = Array{Ti, 2}(undef, dim + 1, ncells) + regions = Array{Ti, 1}(undef, ncells) + for icell = 1:ncells + for inode = 1:(dim + 1) + cells[inode, icell] = parse(Ti, gettoken(tks)) + end + regions[icell] = parse(Ti, gettoken(tks)) + if version20 + for j = 1:(dim + 1) + gettoken(tks) # skip file format garbage + end + end + end + elseif (trytoken(tks, "FACES")) + nfaces = parse(Ti, gettoken(tks)) + faces = Array{Ti, 2}(undef, dim, nfaces) + bregions = Array{Ti, 1}(undef, nfaces) + for iface = 1:nfaces + for inode = 1:dim + faces[inode, iface] = parse(Ti, gettoken(tks)) + end + bregions[iface] = parse(Ti, gettoken(tks)) + if (version20) + for j = 1:(dim + 2) + gettoken(tks) #skip file format garbage + end + end + end + else + expecttoken(tks, "END") + break + end + end + simplexgrid(coord, cells, regions, faces, bregions) end + +function simplexgrid_from_gmsh end + +function simplexgrid_to_gmsh end + +function mixedgrid_from_gmsh end + +function mixedgrid_to_gmsh end diff --git a/src/seal.jl b/src/seal.jl new file mode 100644 index 00000000..47726f75 --- /dev/null +++ b/src/seal.jl @@ -0,0 +1,242 @@ +using StatsBase: countmap + +""" +```` +function encode(x::Vector, nn::Integer) +```` + +Encode th vector `x` into an Int `y`. +The en/-decoding is similar to using the base-`nn` number system. +Example: +`` [x₁, x₂, x₃] → (x₁-1) + (x₂-1)*nn + (x₃-1)*nn²```` +""" +function encode(x::Vector, nn::Integer) + y = 0*nn + for i = 1:length(x) + y += (x[i]-1) * nn^(i - 1) + end + return y +end + +""" +```` +function faces_of_ndim_simplex(x::Vector, dim::Integer, nn::Integer) +```` + +Return all faces of a n-dim simplex. +The orientation is not guaranteed to be right. +`x` contains the nodes of the simplex. +`nn` is the total number of nodes. +The faces (=the nodes contained in the face), are encoded to Integers (of `nn`'s type). + + +""" +function faces_of_ndim_simplex(x::Vector, dim::Integer, nn::Integer) + # for a given array x with n entries. for n=4: x1, x2, x3, x4 + # i want all subsets with length n-1 (and distinct entries) + if dim == 3 # =faces_of_tetrahedron + y = [ + encode(sort([x[1], x[2], x[3]]), nn), + encode(sort([x[1], x[2], x[4]]), nn), + encode(sort([x[1], x[3], x[4]]), nn), + encode(sort([x[2], x[3], x[4]]), nn), + ] + elseif dim == 2 + y = [ + encode(sort([x[1], x[2]]), nn), + encode(sort([x[1], x[3]]), nn), + encode(sort([x[2], x[3]]), nn), + ] + end + return y + +end + +""" +```` +function faces_of_ndim_simplex(x::Vector, dim::Integer, nn::Integer) +```` + +Return all faces of a n-dim simplex. +The orientation is not guaranteed to be right. +`x` contains the nodes of the simplex. +`nn` is the total number of nodes. +The faces (=the nodes contained in the face), are not encoded to Integers. +""" +function faces_of_ndim_simplex_direct(x::Vector) + # for a given array x with n entries. for n=4: x1, x2, x3, x4 + # i want all subsets with length n-1 (and distinct entries) + if length(x) == 4 # =faces_of_tetrahedron + y = [ + sort([x[1], x[2], x[3]]), + sort([x[1], x[2], x[4]]), + sort([x[1], x[3], x[4]]), + sort([x[2], x[3], x[4]]), + ] + elseif length(x) == 3 + y = [ + sort([x[1], x[2]]), + sort([x[1], x[3]]), + sort([x[2], x[3]]), + ] + end + return y +end + +""" +```` +function decode(y::Integer, nn::Integer, dim::Integer) +```` + +Decode `y` to the vector `x`. +`x` has the length `dim`. +The en/-decoding is similar to using the base-`nn` number system. +For details of the encoding, see the documentation of the function `encode`. + +""" +function decode(y::Integer, nn::Integer, dim::Integer) + x = zeros(typeof(y), dim) + x[1] = y % nn + 1 + x[2] = typeof(y).((y ÷ nn) % nn) + 1 + if dim == 3 + x[3] = typeof(y).(y ÷ nn^2) + 1 + end + return x +end + +""" +```` +function assemble_bfaces(simplices, dim, nn, Ti) +```` + +Assemble the BoundaryFaces corresponding to the simplices passed. +In this function, the faces are encoded for performance reasons. If a large grid with many nodes is used, `Ti` has to be chosen accordingly (e.g. `Int128`), or `encode=false` has to be passed to `seal!`. +`simplices` is a ``(dim+1) x 'number cells'`` matrix and `nn` is the total number of nodes. +We can not guarantee, that the orientation of the BoundaryFaces is correct. +""" +function assemble_bfaces(simplices, dim, nn, Ti) + m = size(simplices, 2) + poss_faces = zeros(Ti, (dim + 1) * m) + for i = 1:m + poss_faces[(dim+1)*i-dim:(dim+1)*i] = + faces_of_ndim_simplex(simplices[:, i], dim, nn) + end + dict = countmap(poss_faces) + + k = 0 + for d in dict + (a, b) = d + if b == 1 + k += 1 + #push!(unicats, a) + end + end + + bfaces = zeros(Ti, (dim, k)) + k = 1 + for d in dict + (a, b) = d + if b == 1 + bfaces[:, k] = decode(a, nn, dim) + k += 1 + end + end + + if dim == 3 + @warn "bfaces may not be oriented correctly" + end + + return bfaces +end + + +""" +```` +function assemble_bfaces_direct(simplices, dim, Ti) +```` + +Assemble the BoundaryFaces corresponding to the simplices passed. +In this function, the faces are not encoded. This may make sense for grids with many nodes. +For smaller grids it can lead to performance losses. +`simplices` is a ``(dim+1) x 'number cells'`` matrix and nn is the total number of nodes. +We can not guarantee, that the orientation of the BoundaryFaces is correct. +""" +function assemble_bfaces_direct(simplices, dim, Ti) + m = size(simplices, 2) + poss_faces = fill(zeros(Ti, dim), (dim+1)*m)#zeros(Int64, (dim, (dim + 1)*m)) + for i = 1:m + poss_faces[(dim+1)*i-dim:(dim+1)*i] = faces_of_ndim_simplex_direct(simplices[:, i]) + end + dict = countmap(poss_faces) + + k = 0 + for d in dict + (a, b) = d + if b == 1 + k += 1 + #push!(unicats, a) + end + end + + bfaces = zeros(Ti, (dim, k)) + k = 1 + for d in dict + (a, b) = d + if b == 1 + bfaces[:, k] = a #decode(a, nn, dim) + k += 1 + end + end + + if dim == 3 + @warn "bfaces may not be oriented correctly" + end + + return bfaces +end + + + +""" +```` +function seal!(grid::ExtendableGrid; bfaceregions=[], encode=true, Ti=Int64) +```` + +Take an (simplex-) ExtendableGrid and compute and add the BoundaryFaces. +A so called incomplete ExtendableGrid can e.g. be read from an msh file using the Gmsh.jl-extension of the ExtendableGrids package and the function ````simplexgrid_from_gmsh(filename::String; incomplete=true)````. +If a non empty vector is passed as bfaceregions, this vector is used for the 'BFaceRegions'. +If bfaceregions is empty, all BoundaryFaces get the region number 1. + +For performance reasons, the faces (=the nodes contained in the face) can be encoded (see the function ````encode(x::Vector, nn::Integer)````) to Integers `encoding_type`. To do this, `encode=true` is used. +But for each `encoding_type` there is a limit on the number of nodes: \n + - For Int64 and a 2d grid: 3*10^9 nodes + - For Int64 and a 3d grid: 2*10^6 nodes + - For Int128 and a 2d grid: 1.3*10^19 nodes + - For Int128 and a 3d grid: 5.5*10^12 nodes + +If `encode=false` is passed, there is no limit (besides the MaxValue of the Integer type used). + +""" +function seal!(grid::ExtendableGrid; bfaceregions=[], encode=true, encoding_type=Int64) + dim = size(grid[Coordinates])[1] + Ti2 = typeof(grid[CellNodes][1,1]) + if encode + grid[BFaceNodes] = convert(Matrix{Ti2}, assemble_bfaces(grid[CellNodes], dim, size(grid[Coordinates])[2], encoding_type)) + else + grid[BFaceNodes] = convert(Matrix{Ti2}, assemble_bfaces_direct(grid[CellNodes], dim, encoding_type)) + end + + if bfaceregions==[] + grid[BFaceRegions] = ones(Ti2, size(grid[BFaceNodes])[2]) + else + grid[BFaceRegions] = bfaceregions #Int64.(collect(1:size(grid[BFaceNodes])[2])) + end + + if dim == 2 + grid[BFaceGeometries] = VectorOfConstants{ElementGeometries,Ti2}(Edge1D, size(grid[BFaceNodes])[2]) + else + grid[BFaceGeometries] = VectorOfConstants{ElementGeometries,Ti2}(Triangle2D, size(grid[BFaceNodes])[2]) + end + + return grid +end diff --git a/src/simplexgrid.jl b/src/simplexgrid.jl index d7537006..ad8fe4bb 100644 --- a/src/simplexgrid.jl +++ b/src/simplexgrid.jl @@ -644,163 +644,6 @@ function simplexgrid(grid2::ExtendableGrid, coordZ; bot_offset=0,cell_offset=0,t end -###################################################### -""" -$(TYPEDSIGNATURES) - -Read grid from file. Currently for pdelib sg format only. -""" -function simplexgrid(file::String;format="") - Ti=Cint - (fbase,fext)=splitext(file) - if format=="" - format=fext[2:end] - end - @assert format=="sg" - - tks=TokenStream(file) - expecttoken(tks,"SimplexGrid") - version=parse(Float64,gettoken(tks)) - version20=false; - - if (version==2.0) - version20=true; - elseif (version==2.1) - version20=false; - else - error("Read grid: wrong format version: $(version)") - end - - dim::Ti=0 - coord=Array{Float64,2}(undef,0,0) - cells=Array{Ti,2}(undef,0,0) - regions=Array{Ti,1}(undef,0) - faces=Array{Ti,2}(undef,0,0) - bregions=Array{Ti,1}(undef,0) - while(true) - if (trytoken(tks,"DIMENSION")) - dim=parse(Ti,gettoken(tks)); - elseif (trytoken(tks,"NODES")) - nnodes=parse(Ti,gettoken(tks)); - embdim=parse(Ti,gettoken(tks)); - if(dim!=embdim) - error("Dimension error (DIMENSION $(dim)) in section NODES") - end - coord=Array{Float64,2}(undef,dim,nnodes) - for inode=1:nnodes - for idim=1:embdim - coord[idim,inode]=parse(Float64,gettoken(tks)) - end - end - elseif (trytoken(tks,"CELLS")) - ncells=parse(Ti,gettoken(tks)); - cells=Array{Ti,2}(undef,dim+1,ncells) - regions=Array{Ti,1}(undef,ncells) - for icell=1:ncells - for inode=1:dim+1 - cells[inode,icell]=parse(Ti,gettoken(tks)); - end - regions[icell]=parse(Ti,gettoken(tks)); - if version20 - for j=1:dim+1 - gettoken(tks); # skip file format garbage - end - end - end - elseif (trytoken(tks,"FACES")) - nfaces=parse(Ti,gettoken(tks)); - faces=Array{Ti,2}(undef,dim,nfaces) - bregions=Array{Ti,1}(undef,nfaces) - for iface=1:nfaces - for inode=1:dim - faces[inode,iface]=parse(Ti,gettoken(tks)); - end - bregions[iface]=parse(Ti,gettoken(tks)); - if (version20) - for j=1:dim+2 - gettoken(tks); #skip file format garbage - end - end - end - else - expecttoken(tks,"END") - break - end - end - simplexgrid(coord,cells,regions,faces,bregions); -end - - -# Implementations in Gmsh extension -function simplexgrid_from_gmsh end - -function write_gmsh end - - -""" -$(TYPEDSIGNATURES) - -Write grid to file. Currently for pdelib sg format only -""" -function Base.write(fname::String, g::ExtendableGrid; format="") - (fbase,fext)=splitext(fname) - if format=="" - format=fext[2:end] - end - @assert format=="sg" - - dim_g=dim_grid(g) - dim_s=dim_space(g) - nn=num_nodes(g) - nc=num_cells(g) - nbf=num_bfaces(g) - coord=g[Coordinates] - cellnodes=g[CellNodes] - bfacenodes=g[BFaceNodes] - cellregions=g[CellRegions] - bfaceregions=g[BFaceRegions] - - # TODO: replace @sprintf by someting non-allocating - open(fname, "w") do file - write(file,@sprintf("SimplexGrid")) - write(file,@sprintf(" ")) - write(file,@sprintf("2.1\n")) - write(file,@sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n")) - write(file,@sprintf("#mailto:{fuhrmann|streckenbach}@wias-berlin.de\n")) - write(file,@sprintf("#%s\n",Dates.format(Dates.now(),"yyyy-mm-ddTHH-mm-SS"))) - - write(file,@sprintf("DIMENSION\n%d\n",dim_g)) - write(file,@sprintf("NODES\n%d %d\n",nn,dim_s)) - - for inode=1:nn - for idim=1:dim_s - write(file,@sprintf("%.20e ",coord[idim,inode])) - write(file,@sprintf("\n")) - end - end - - write(file,@sprintf("CELLS\n%d\n",nc)) - for icell=1:nc - for inode=1:dim_g+1 - write(file,@sprintf("%d ", cellnodes[inode,icell])) - end - write(file,@sprintf("%d\n",cellregions[icell])) - end - - write(file,@sprintf("FACES\n%d\n",nbf)) - for ibface=1:nbf - for inode=1:dim_g - write(file,@sprintf("%d ", bfacenodes[inode,ibface])) - end - write(file,@sprintf("%d\n",bfaceregions[ibface])) - end - write(file,@sprintf("END\n")) - flush(file) - flush(file) - end - nothing -end - """ diff --git a/test/gmsh.jl b/test/gmsh.jl new file mode 100644 index 00000000..13733634 --- /dev/null +++ b/test/gmsh.jl @@ -0,0 +1,139 @@ +using Gmsh: gmsh +gmsh.initialize() + +@testset "access gmsh" begin + gmsh.option.setNumber("General.Terminal", 1) + gmsh.model.add("t1") + + lc = 1e-2 + gmsh.model.geo.addPoint(0, 0, 0, lc, 1) + gmsh.model.geo.addPoint(0.1, 0, 0, lc, 2) + gmsh.model.geo.addPoint(0.1, 0.3, 0, lc, 3) + + p4 = gmsh.model.geo.addPoint(0, 0.3, 0, lc) + + gmsh.model.geo.addLine(1, 2, 1) + gmsh.model.geo.addLine(3, 2, 2) + gmsh.model.geo.addLine(3, p4, 3) + gmsh.model.geo.addLine(4, 1, p4) + + gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1) + gmsh.model.geo.addPlaneSurface([1], 1) + + gmsh.model.geo.synchronize() + + gmsh.model.addPhysicalGroup(0, [1, 2], 1) + gmsh.model.addPhysicalGroup(1, [1, 2], 2) + gmsh.model.addPhysicalGroup(2, [1], 6) + + gmsh.model.setPhysicalName(2, 6, "My surface") + + gmsh.model.mesh.generate(2) + grid = ExtendableGrids.simplexgrid_from_gmsh(gmsh.model) + + @test num_nodes(grid) > 0 && num_cells(grid) > 0 && num_bfaces(grid) > 0 + + # @test testgrid(grid, (404, 726, 80)) gmsh generates differently on windows + gmsh.clear() +end + +@testset "Read/write simplex gmsh 2d / 3d" begin + gmsh.option.setNumber("General.Terminal", 0) + gmsh.option.setNumber("General.Verbosity", 0) + + path = "" + + X = collect(0:0.02:2) + Y = collect(0:2:4) + grid1 = simplexgrid(X, Y) #ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh") + ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh") + grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float32, Ti = Int64) + #gmsh.finalize() + @test seemingly_equal(grid2, grid1; sort = true, confidence = :low) + @test seemingly_equal(grid2, grid1; sort = true, confidence = :full) + + gmsh.clear() + + grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_2d.msh"; Tc = Float64, Ti = Int64) + gmsh.clear() + ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh") + grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float64, Ti = Int64) + #gmsh.finalize() + + @test seemingly_equal(grid1, grid2; sort = true, confidence = :low) + @test seemingly_equal(grid1, grid2; sort = true, confidence = :full) + + gmsh.clear() + + grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_3d.msh"; Tc = Float32, Ti = Int64) + gmsh.clear() + ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = path * "testfile.msh") + grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "testfile.msh"; Tc = Float64, Ti = Int32) + #gmsh.finalize() + + @test seemingly_equal(grid1, grid2; sort = true, confidence = :low) + @test seemingly_equal(grid1, grid2; sort = true, confidence = :full) + + gmsh.clear() + + grid1 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_2d.msh") + gmsh.clear() + #grid2 = + #simplexgrid([0, 1, 2], [3, 4, 5]) + grid2 = ExtendableGrids.simplexgrid_from_gmsh(path * "sto_3d.msh"; Tc = Float32, Ti = Int32) + #gmsh.finalize() + + @test !seemingly_equal(grid1, grid2; sort = true, confidence = :low) + @test !seemingly_equal(grid1, grid2; sort = true, confidence = :full) + + gmsh.clear() + + grid1 = ExtendableGrids.simplexgrid_from_gmsh("testmesh.gmsh"; incomplete = true) + ExtendableGrids.seal!(grid1; encode = false) + gmsh.clear() + ExtendableGrids.simplexgrid_to_gmsh(grid1; filename = "completed_testfile.msh") + grid2 = ExtendableGrids.simplexgrid_from_gmsh("completed_testfile.msh") + + gmsh.clear() + + grid3 = ExtendableGrids.simplexgrid_from_gmsh("testmesh.gmsh"; incomplete = true) + ExtendableGrids.seal!(grid3; encode = true) + + gmsh.clear() + + @test seemingly_equal(grid1, grid2; sort = true, confidence = :low) + @test seemingly_equal(grid1, grid2; sort = true, confidence = :full) + @test seemingly_equal(grid1, grid3; sort = true, confidence = :low) + @test seemingly_equal(grid1, grid3; sort = true, confidence = :full) + + x = collect(LinRange(0, 1, 50)) + grid1 = simplexgrid(x, x) + grid1[BFaceRegions] = ones(Int32, length(grid1[BFaceRegions])) #num_faces(grid1)) + grid2 = simplexgrid(x, x) + grid3 = simplexgrid(x, x) + ExtendableGrids.seal!(grid2) + ExtendableGrids.seal!(grid3; encode = false) + + gmsh.finalize() + + @test seemingly_equal(grid2, grid1; sort = true, confidence = :low) + @test seemingly_equal(grid2, grid1; sort = true, confidence = :full) + @test seemingly_equal(grid3, grid1; sort = true, confidence = :low) + @test seemingly_equal(grid3, grid1; sort = true, confidence = :full) +end + +@testset "Read/write mixed gmsh 2d" begin + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 0) + gmsh.option.setNumber("General.Verbosity", 0) + + path = "" + grid1 = ExtendableGrids.mixedgrid_from_gmsh(path * "mixedgrid_2d.msh"; Tc = Float64, Ti = Int64) + gmsh.clear() + ExtendableGrids.mixedgrid_to_gmsh(grid1; filename = path * "testfile.msh") + grid2 = ExtendableGrids.mixedgrid_from_gmsh(path * "testfile.msh"; Tc = Float32, Ti = UInt64) + gmsh.finalize() + + @test seemingly_equal(grid1, grid2; sort = true, confidence = :low) + @test seemingly_equal(grid1, grid2; sort = true, confidence = :full) +end diff --git a/test/mixedgrid_2d.msh b/test/mixedgrid_2d.msh new file mode 100644 index 00000000..4ff16ef2 --- /dev/null +++ b/test/mixedgrid_2d.msh @@ -0,0 +1,2339 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +11 10 2 0 +1 0 0 0 0 +2 0.1 0 0 0 +3 0.1 0.3 0 0 +4 0 0.3 0 0 +5 -0.05 0.05 0 0 +6 -0.05 0.1 0 0 +7 0.2 0.2 0 0 +8 0.2 0.1 0 0 +9 0 0.3 0 0 +10 0.25 0.2 0 0 +11 0.3 0.1 0 0 +1 0 0 0 0.1 0 0 0 2 1 -2 +2 0.1 0 0 0.1 0.3 0 0 2 3 -2 +3 0 0.3 0 0.1 0.3 0 0 2 3 -4 +4 -0.05 0 0 0 0.05 0 0 2 1 -5 +5 -0.05 0.05 0 -0.05 0.1 0 0 2 5 -6 +6 -0.05 0.1 0 0 0.3 0 0 2 6 -4 +10 0.2 0.1 0 0.3 0.1 0 0 2 8 -11 +11 0.25 0.1 0 0.3 0.2 0 0 2 11 -10 +12 0.2 0.2 0 0.25 0.2 0 0 2 10 -7 +13 0.2 0.1 0 0.2 0.2 0 0 2 7 -8 +1 -0.05000000000000001 0 0 0.1 0.3 0 0 6 3 -6 -5 -4 1 -2 +15 0.2 0.1 0 0.3 0.2 0 0 4 13 10 11 12 +$EndEntities +$Nodes +23 701 1 701 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +0.1 0 0 +0 3 0 1 +3 +0.1 0.3 0 +0 4 0 1 +4 +0 0.3 0 +0 5 0 1 +5 +-0.05 0.05 0 +0 6 0 1 +6 +-0.05 0.1 0 +0 7 0 1 +7 +0.2 0.2 0 +0 8 0 1 +8 +0.2 0.1 0 +0 9 0 1 +9 +0 0.3 0 +0 10 0 1 +10 +0.25 0.2 0 +0 11 0 1 +11 +0.3 0.1 0 +1 1 0 28 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +0.01675134899749656 0 0 +0.03071080655893446 0 0 +0.0423436878825755 0 0 +0.05203775563819656 0 0 +0.06011614539899403 0 0 +0.06684813691251688 0 0 +0.07245812988783501 0 0 +0.07713312407354976 0 0 +0.08102895247307268 0 0 +0.08427547609907809 0 0 +0.08698091248825451 0 0 +0.0892354428974928 0 0 +0.09111421828692554 0 0 +0.09267986447704872 0 0 +0.09398456962083367 0 0 +0.0950718239445128 0 0 +0.09597786913611506 0 0 +0.09673290670057005 0 0 +0.09736210469696704 0 0 +0.09788643639360275 0 0 +0.09832337953171119 0 0 +0.09868749882758426 0 0 +0.09899093151065048 0 0 +0.09924379207211487 0 0 +0.09945450923783336 0 0 +0.09963010681159246 0 0 +0.09977643815222635 0 0 +0.0998983809566317 0 0 +1 2 0 18 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +0.1 0.2842105263158026 0 +0.1 0.268421052631616 0 +0.1 0.2526315789474643 0 +0.1 0.236842105263308 0 +0.1 0.2210526315791431 0 +0.1 0.2052631578949782 0 +0.1 0.1894736842108125 0 +0.1 0.1736842105266377 0 +0.1 0.1578947368424728 0 +0.1 0.142105263158266 0 +0.1 0.1263157894740112 0 +0.1 0.1105263157897576 0 +0.1 0.09473684210550881 0 +0.1 0.07894736842125294 0 +0.1 0.06315789473700048 0 +0.1 0.04736842105275713 0 +0.1 0.03157894736850836 0 +0.1 0.01578947368424877 0 +1 3 0 28 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +0.09989838095666016 0.3 0 +0.09977643815229478 0.3 0 +0.09963010681171565 0.3 0 +0.09945450923803045 0.3 0 +0.09924379207241064 0.3 0 +0.0989909315110763 0.3 0 +0.09868749882818036 0.3 0 +0.09832337953252829 0.3 0 +0.09788643639470557 0.3 0 +0.09736210469843745 0.3 0 +0.09673290670250986 0.3 0 +0.09597786913865393 0.3 0 +0.09507182394781449 0.3 0 +0.09398456962510122 0.3 0 +0.0926798644825377 0.3 0 +0.09111421829465347 0.3 0 +0.08923544290644983 0.3 0 +0.08698091249963096 0.3 0 +0.08427547611348578 0.3 0 +0.08102895249127755 0.3 0 +0.07713312409649188 0.3 0 +0.07245812991955933 0.3 0 +0.06684813694870104 0.3 0 +0.06011614544430099 0.3 0 +0.05203775569482921 0.3 0 +0.04234368796022012 0.3 0 +0.03071080664675478 0.3 0 +0.0167513491066254 0.3 0 +1 4 0 4 +86 +87 +88 +89 +-0.009999999999978011 0.009999999999978011 0 +-0.01999999999994795 0.01999999999994795 0 +-0.02999999999994715 0.02999999999994715 0 +-0.03999999999997372 0.03999999999997372 0 +1 5 0 4 +90 +91 +92 +93 +-0.05 0.05999999999997289 0 +-0.05 0.06999999999994574 0 +-0.05 0.07999999999994248 0 +-0.05 0.08999999999997124 0 +1 6 0 8 +94 +95 +96 +97 +98 +99 +100 +101 +-0.04444444444445491 0.1222222222221804 0 +-0.0388888888889152 0.1444444444443392 0 +-0.03333333333337599 0.166666666666496 0 +-0.02777777777783678 0.1888888888886529 0 +-0.02222222222228109 0.2111111111108757 0 +-0.0166666666667102 0.2333333333331592 0 +-0.01111111111114187 0.2555555555554325 0 +-0.005555555555571662 0.2777777777777133 0 +1 10 0 8 +102 +103 +104 +105 +106 +107 +108 +109 +0.211111111111087 0.1 0 +0.222222222222174 0.1 0 +0.2333333333332502 0.1 0 +0.2444444444443263 0.1 0 +0.2555555555554458 0.1 0 +0.2666666666665871 0.1 0 +0.2777777777777283 0.1 0 +0.2888888888888672 0.1 0 +1 11 0 8 +110 +111 +112 +113 +114 +115 +116 +117 +0.2944444444444635 0.1111111111110731 0 +0.2888888888889183 0.1222222222221634 0 +0.2833333333333687 0.1333333333332625 0 +0.2777777777778195 0.1444444444443611 0 +0.2722222222222627 0.1555555555554746 0 +0.2666666666666876 0.1666666666666247 0 +0.2611111111111217 0.1777777777777565 0 +0.2555555555555609 0.1888888888888783 0 +1 12 0 8 +118 +119 +120 +121 +122 +123 +124 +125 +0.2444444444444607 0.2 0 +0.2388888888889268 0.2 0 +0.233333333333393 0.2 0 +0.2277777777778266 0.2 0 +0.2222222222222602 0.2 0 +0.2166666666666938 0.2 0 +0.2111111111111298 0.2 0 +0.2055555555555676 0.2 0 +1 13 0 8 +126 +127 +128 +129 +130 +131 +132 +133 +0.2 0.1888888888889136 0 +0.2 0.1777777777778405 0 +0.2 0.1666666666667607 0 +0.2 0.1555555555556809 0 +0.2 0.1444444444445626 0 +0.2 0.1333333333334256 0 +0.2 0.122222222222284 0 +0.2 0.1111111111111427 0 +2 1 0 504 +134 +135 +136 +137 +138 +139 +140 +141 +142 +143 +144 +145 +146 +147 +148 +149 +150 +151 +152 +153 +154 +155 +156 +157 +158 +159 +160 +161 +162 +163 +164 +165 +166 +167 +168 +169 +170 +171 +172 +173 +174 +175 +176 +177 +178 +179 +180 +181 +182 +183 +184 +185 +186 +187 +188 +189 +190 +191 +192 +193 +194 +195 +196 +197 +198 +199 +200 +201 +202 +203 +204 +205 +206 +207 +208 +209 +210 +211 +212 +213 +214 +215 +216 +217 +218 +219 +220 +221 +222 +223 +224 +225 +226 +227 +228 +229 +230 +231 +232 +233 +234 +235 +236 +237 +238 +239 +240 +241 +242 +243 +244 +245 +246 +247 +248 +249 +250 +251 +252 +253 +254 +255 +256 +257 +258 +259 +260 +261 +262 +263 +264 +265 +266 +267 +268 +269 +270 +271 +272 +273 +274 +275 +276 +277 +278 +279 +280 +281 +282 +283 +284 +285 +286 +287 +288 +289 +290 +291 +292 +293 +294 +295 +296 +297 +298 +299 +300 +301 +302 +303 +304 +305 +306 +307 +308 +309 +310 +311 +312 +313 +314 +315 +316 +317 +318 +319 +320 +321 +322 +323 +324 +325 +326 +327 +328 +329 +330 +331 +332 +333 +334 +335 +336 +337 +338 +339 +340 +341 +342 +343 +344 +345 +346 +347 +348 +349 +350 +351 +352 +353 +354 +355 +356 +357 +358 +359 +360 +361 +362 +363 +364 +365 +366 +367 +368 +369 +370 +371 +372 +373 +374 +375 +376 +377 +378 +379 +380 +381 +382 +383 +384 +385 +386 +387 +388 +389 +390 +391 +392 +393 +394 +395 +396 +397 +398 +399 +400 +401 +402 +403 +404 +405 +406 +407 +408 +409 +410 +411 +412 +413 +414 +415 +416 +417 +418 +419 +420 +421 +422 +423 +424 +425 +426 +427 +428 +429 +430 +431 +432 +433 +434 +435 +436 +437 +438 +439 +440 +441 +442 +443 +444 +445 +446 +447 +448 +449 +450 +451 +452 +453 +454 +455 +456 +457 +458 +459 +460 +461 +462 +463 +464 +465 +466 +467 +468 +469 +470 +471 +472 +473 +474 +475 +476 +477 +478 +479 +480 +481 +482 +483 +484 +485 +486 +487 +488 +489 +490 +491 +492 +493 +494 +495 +496 +497 +498 +499 +500 +501 +502 +503 +504 +505 +506 +507 +508 +509 +510 +511 +512 +513 +514 +515 +516 +517 +518 +519 +520 +521 +522 +523 +524 +525 +526 +527 +528 +529 +530 +531 +532 +533 +534 +535 +536 +537 +538 +539 +540 +541 +542 +543 +544 +545 +546 +547 +548 +549 +550 +551 +552 +553 +554 +555 +556 +557 +558 +559 +560 +561 +562 +563 +564 +565 +566 +567 +568 +569 +570 +571 +572 +573 +574 +575 +576 +577 +578 +579 +580 +581 +582 +583 +584 +585 +586 +587 +588 +589 +590 +591 +592 +593 +594 +595 +596 +597 +598 +599 +600 +601 +602 +603 +604 +605 +606 +607 +608 +609 +610 +611 +612 +613 +614 +615 +616 +617 +618 +619 +620 +621 +622 +623 +624 +625 +626 +627 +628 +629 +630 +631 +632 +633 +634 +635 +636 +637 +0.002179647562980582 0.01118005220911409 0 +-0.00800298242437353 0.02170195833185488 0 +-0.01691807708238987 0.03211374670039643 0 +-0.0248237622169816 0.04251124252348602 0 +-0.0313654530998219 0.05290575550097789 0 +-0.03484458046915367 0.06330783482289694 0 +-0.03662070357167233 0.07384006190227238 0 +-0.03750901707896619 0.08468206036292211 0 +-0.03767684093268259 0.09629313598115957 0 +-0.03655621002275887 0.1103602252490067 0 +-0.03310596515313692 0.1295155064560043 0 +-0.02898027398901567 0.1498020195056517 0 +-0.02440033012969737 0.1707667663956995 0 +-0.01953243876019263 0.1921025241925427 0 +-0.01444928021950475 0.2136467894633583 0 +-0.009131184376161674 0.235317761099278 0 +-0.003409598875143667 0.257068529279962 0 +0.003426242529007883 0.2788255638099427 0 +0.01390097468043122 0.01193810001986666 0 +0.002958032631695102 0.02301774031000928 0 +-0.005739964885226468 0.03384989169935685 0 +-0.0129987472304506 0.04463430300076483 0 +-0.01878551562893036 0.05543890868024167 0 +-0.02268922867761609 0.06632608042650186 0 +-0.02507835654394779 0.07745648277312746 0 +-0.02637114049314996 0.08911712432740891 0 +-0.02666525836281414 0.1018939772724562 0 +-0.02568722065136964 0.1168950762441364 0 +-0.02321055420382989 0.1348559513061863 0 +-0.01989073103670397 0.1541478286958064 0 +-0.01596653919172093 0.1742933946289651 0 +-0.01161009568028275 0.194955750591378 0 +-0.006881511436553083 0.215935472743143 0 +-0.00169006265842702 0.2371191531509951 0 +0.004373678703254783 0.2584248126042158 0 +0.01289824639751419 0.279669624931562 0 +0.02434420541830997 0.01248618438075972 0 +0.01286183901301944 0.02405681602747239 0 +0.004161286760525513 0.03529778965791809 0 +-0.002828086268413696 0.04645995549868516 0 +-0.008329480345055314 0.05766084915353481 0 +-0.01229910008065249 0.06900578002222946 0 +-0.01491756229680492 0.08067854622483406 0 +-0.01641195394714708 0.09299025859227886 0 +-0.01684810221562147 0.1064583142669834 0 +-0.01611784282628983 0.1217639253279337 0 +-0.01419348059470044 0.1391520652482346 0 +-0.01140619655246725 0.1578260075774027 0 +-0.007961036409564781 0.1773782833384425 0 +-0.003998145008277111 0.1975023012794415 0 +0.0004638998653294917 0.2179993590414715 0 +0.005610318859617582 0.2387434866465135 0 +0.01210255642729585 0.2596200805167619 0 +0.02210667990667416 0.2803493900898861 0 +0.0334395223703497 0.01291187942416287 0 +0.0216949429688003 0.02490136804298277 0 +0.01299803854389844 0.0365213283370897 0 +0.006170309226767003 0.04804212318237827 0 +0.0008388028111182838 0.05961449897637104 0 +-0.003103169553144076 0.07137275721404282 0 +-0.005796405618980849 0.08351104041475857 0 +-0.007388027870994789 0.09632611528902815 0 +-0.00793387002764951 0.1102380038573932 0 +-0.00740521269375061 0.1256837054584643 0 +-0.005835143271244866 0.1427554911032567 0 +-0.003428157136594011 0.1610051598632215 0 +-0.0003485064417049438 0.1801024225230048 0 +0.003303999549049702 0.1997817282189889 0 +0.007561835354348187 0.2198574137022944 0 +0.01270704245709834 0.2401985714357359 0 +0.01960370115778175 0.2606614193842658 0 +0.03065802563826879 0.2809010611361263 0 +0.04131358104631582 0.0132568837725304 0 +0.02953181632993396 0.02560512712243303 0 +0.02091105761117688 0.03756926725151683 0 +0.01423678666381627 0.04942411369063969 0 +0.00905723492861365 0.06133985443988713 0 +0.00519176105551392 0.07346803680088999 0 +0.002504768021727786 0.08600371053960369 0 +0.0008818993848456548 0.09921636478194852 0 +0.0002679988426741984 0.113441495866599 0 +0.0006513433910077077 0.1289742841121188 0 +0.001976910578381826 0.1458513765504801 0 +0.004099085119288068 0.1637897332262052 0 +0.006890517824193866 0.1825235644590882 0 +0.01029104919035112 0.2018259364707103 0 +0.01438371624845215 0.2215273431738179 0 +0.01953397340772282 0.2414952759604224 0 +0.02673418042601322 0.2615632170574396 0 +0.03840091607153021 0.2813529135333104 0 +0.04812875874917764 0.01354421340019939 0 +0.03648088479829401 0.02620335617209428 0 +0.02801825708341231 0.03847834182874154 0 +0.02152481049192824 0.05064145685424133 0 +0.01651162433620201 0.06287275506287031 0 +0.01275959217524276 0.07533238195578945 0 +0.01013032872995907 0.08821027066514159 0 +0.008522264158659705 0.1017462958884348 0 +0.007877170100683076 0.116209370277157 0 +0.008160873980118236 0.1318076494391748 0 +0.009309781358216642 0.1485539662092931 0 +0.01121077251364322 0.1662520140292546 0 +0.01376766169336346 0.1846860185676463 0 +0.01695583976944366 0.2036620558550011 0 +0.02090260677796698 0.22302645403607 0 +0.02603510851226126 0.2426471128415678 0 +0.03340376275974218 0.2623428558294936 0 +0.04531416893099043 0.2817265106267423 0 +0.05403951535165688 0.01378810346090621 0 +0.04265391722623377 0.026720086343659 0 +0.03442272676435045 0.03927623657853038 0 +0.02814675899314787 0.0517229780683675 0 +0.02332391583056972 0.06424391593730808 0 +0.01971572879639976 0.07700175965244151 0 +0.01718022972283139 0.09017776933509727 0 +0.01562012053853164 0.1039837094494521 0 +0.01497290582651994 0.1186369645909062 0 +0.01519104001283466 0.1342890579718985 0 +0.01620793834098983 0.1509407387957862 0 +0.01793142765430717 0.1684454213790176 0 +0.0202910580925041 0.1866253570556461 0 +0.02329226263876477 0.2053138964697393 0 +0.02709737171398965 0.2243717853323951 0 +0.03217184520902541 0.2436690565829768 0 +0.03957097390077857 0.2630178228592762 0 +0.0514401618984137 0.282038326745054 0 +0.05918181728061671 0.0139981001251248 0 +0.04815421406153512 0.02717235183684442 0 +0.04021445246654045 0.03998400171004208 0 +0.03419118142523277 0.05269187337514096 0 +0.02958303288794185 0.06547898013177181 0 +0.02614325866287524 0.07850648526164129 0 +0.02372740805552581 0.09194501795926409 0 +0.0222389563285111 0.1059808850582408 0 +0.02161118399370873 0.1207913329877215 0 +0.02178830557141464 0.1364890457813579 0 +0.02270353214415203 0.1530676174630632 0 +0.02427989890136046 0.1704112001903821 0 +0.02646793449253029 0.1883710009098085 0 +0.02929717643485207 0.206802684303339 0 +0.03295494708007218 0.2255798544952074 0 +0.03792420248992177 0.2445763176252503 0 +0.04523038139232315 0.2636041991684693 0 +0.05684887058868479 0.2823010143207051 0 +0.06367216057728574 0.01418098333201086 0 +0.05307318788855375 0.02757264183448337 0 +0.04547150546407001 0.04061783752061166 0 +0.03973063511141941 0.0535667633759256 0 +0.03535811414997617 0.06659907751070084 0 +0.03210560572175121 0.07987163022445375 0 +0.02982761109879924 0.09354346155471104 0 +0.02842717021311673 0.1077782189200799 0 +0.02783399278260859 0.1227214467039944 0 +0.02798730877601328 0.1384580157171881 0 +0.02882172077969002 0.1549767851856251 0 +0.03027219006788581 0.1721821963755943 0 +0.03230610821196495 0.1899477560207938 0 +0.0349713431029579 0.2081474251904074 0 +0.0384707626515637 0.2266662958269346 0 +0.04328835747119597 0.2453834234783151 0 +0.05039985126776312 0.2641160628020894 0 +0.06161905770870218 0.2825243470710828 0 +0.06760933513656484 0.01434177512988451 0 +0.05749022036543426 0.02793037256642249 0 +0.05026136049053762 0.04119037227015435 0 +0.04482589483846511 0.05436262056394485 0 +0.04070512669317138 0.06762152859409153 0 +0.03765328217937038 0.08111778737638474 0 +0.03552515336524921 0.09499850982712867 0 +0.03422322245002422 0.1094072894729831 0 +0.03367457031353994 0.1244642589634003 0 +0.03381571415169072 0.140233688891929 0 +0.03458364190286956 0.1567011132311764 0 +0.03592308107040325 0.1737851591094349 0 +0.03781474291287574 0.1913767879126175 0 +0.04031957804419506 0.2093651002344152 0 +0.04364788864276779 0.2276455582122074 0 +0.04827276703893633 0.2461036734016291 0 +0.05511062473183996 0.264565399005587 0 +0.0658287448771855 0.2827159237782923 0 +0.07107680801112973 0.01448431662397343 0 +0.06147368193844674 0.02825281086669599 0 +0.05464225048210236 0.04171157465214 0 +0.04952856011604147 0.05509153551148489 0 +0.0456706011780657 0.06856052538176097 0 +0.04282781482784033 0.0822618603491802 0 +0.04085645340935404 0.09633077360319953 0 +0.03965884972119518 0.1108931307088338 0 +0.03916051202427036 0.1260484239422059 0 +0.03929694059814247 0.1418452793323791 0 +0.04000814939957088 0.1582668606541631 0 +0.04124703875938802 0.1752422354415168 0 +0.04300461581794527 0.1926762864202287 0 +0.04535038253069652 0.2104708054490987 0 +0.04849555460493853 0.2285307261824146 0 +0.05289428591019443 0.2467488991373271 0 +0.05940045018327453 0.264962253579371 0 +0.06955068530951734 0.2828816947074458 0 +0.07414505203957571 0.01461162109016303 0 +0.06508227825485134 0.0285456775953768 0 +0.0586644511313126 0.04218941074816866 0 +0.0538828597730055 0.05576333049539217 0 +0.05029397190519033 0.0694277284231495 0 +0.04766421215519098 0.08331776707113417 0 +0.04585232124011818 0.09755708691040829 0 +0.04476115116392081 0.1122558738421002 0 +0.0443157084189387 0.1274966725517684 0 +0.04445179356616132 0.1433160166437154 0 +0.04511285618134295 0.159695419705736 0 +0.04625870220083009 0.1765719955467914 0 +0.04788808847594328 0.1938619551399599 0 +0.05007531395333965 0.2114778885182339 0 +0.05302753151380696 0.2293334642711903 0 +0.05717492087435765 0.2473294318631995 0 +0.06330917455210645 0.2653149799631227 0 +0.0728505151283121 0.2830263553240891 0 +0.07687361964291338 0.01472610169816973 0 +0.06836640121174795 0.02881355589521004 0 +0.06237144630101973 0.04263032311316621 0 +0.05792700643262803 0.05638604461772036 0 +0.05460917788333696 0.07023276830038794 0 +0.05219262771041623 0.08429703297834042 0 +0.05053951535205446 0.09869132107841484 0 +0.04955399932534097 0.1135119378030687 0 +0.04916161018769126 0.1288273942804095 0 +0.04929948985626259 0.1446647687139985 0 +0.0499147713134313 0.1610045019478238 0 +0.05097311700498732 0.1777901715615471 0 +0.05247892422578227 0.1949473977405415 0 +0.05450828636648374 0.2123980998321736 0 +0.05726066466945016 0.2300640564106744 0 +0.06113938665064205 0.2478541887943178 0 +0.06687609214225598 0.2656305032049472 0 +0.07578631494322785 0.283153640770454 0 +0.07931292361123543 0.01482972461169786 0 +0.07136936878508018 0.02906017577703262 0 +0.06580096270733071 0.04303958534141909 0 +0.06169427232346432 0.056966316378779 0 +0.05864584032784078 0.07098365890051574 0 +0.05643955303586495 0.08520928002815549 0 +0.05494185380621478 0.09974502682885543 0 +0.05405903235262347 0.1146749052555472 0 +0.05371808747851473 0.1300557313645117 0 +0.05385832589307252 0.1459071408203093 0 +0.0544306941918009 0.1622089792508595 0 +0.05540582825769948 0.1789102130217786 0 +0.0567920494909977 0.1959444387650779 0 +0.0586649206560574 0.2132417569996159 0 +0.06121368286540781 0.2307315065605207 0 +0.06481340063068394 0.2483308180974517 0 +0.07013847140572352 0.2659145654771687 0 +0.07840888562738206 0.2832665464738539 0 +0.08150574197942721 0.01492411550527088 0 +0.07412852086059568 0.02928861774335308 0 +0.06898587944360796 0.0434215686143468 0 +0.06521387679859815 0.05750968644362352 0 +0.06243017781663195 0.07168713560458198 0 +0.06042871719019466 0.08606262640693435 0 +0.05908103885925836 0.10072794109023 0 +0.05829637457647781 0.1157561783934571 0 +0.05800403204362055 0.1311943612906091 0 +0.05814612621415768 0.1470562521906564 0 +0.05867746081313071 0.1633215016640572 0 +0.059572899525187 0.1799437187558456 0 +0.06084331740780476 0.1968633980658963 0 +0.0625620028522675 0.2140179152632791 0 +0.06490630712431333 0.2313436716907903 0 +0.06822257315889405 0.2487658636686033 0 +0.07313083801640152 0.2661719398725393 0 +0.08076235683133474 0.2833674940069694 0 +0.08348848083808109 0.01501063586703202 0 +0.07667617122346015 0.02950146194569996 0 +0.07195502394978325 0.0437799452984377 0 +0.06851173617538665 0.05802083885617296 0 +0.06598574522380653 0.07234893170457564 0 +0.06418179001045027 0.08686401249164576 0 +0.0629772866203881 0.1016483888239241 0 +0.06228517182550272 0.116765477920241 0 +0.06203778814552639 0.1322540722149533 0 +0.06218055009419075 0.1481233035431115 0 +0.06267209838562937 0.1643529645344251 0 +0.06349089744708861 0.1809007810895661 0 +0.06464930335553784 0.1977133285590181 0 +0.06621706732356121 0.2147345353916014 0 +0.06835862929407818 0.231907406162594 0 +0.0713917418915901 0.2491649279766508 0 +0.07588472599967096 0.2664066101209802 0 +0.08288492143529914 0.2834584564511845 0 +0.08529223136238168 0.01509043908223047 0 +0.07904042792863872 0.02970089983336297 0 +0.07473386687670872 0.04411784690575466 0 +0.07161110546659667 0.05850379470289757 0 +0.06933404392352209 0.07297400456394369 0 +0.06771894394548469 0.08761946661168725 0 +0.06664981557332271 0.1025136031214573 0 +0.0660439935384757 0.1177112278921832 0 +0.06583746357215119 0.1332441955931286 0 +0.06597930224759103 0.1491180028678257 0 +0.06643192149145753 0.1653128694443468 0 +0.06717686313774388 0.1817902646270524 0 +0.06822714252283575 0.1985022224580199 0 +0.06964809707155654 0.2153986422501784 0 +0.0715907097979336 0.2324287039290729 0 +0.07434461949333758 0.2495328219537908 0 +0.07842870675340376 0.2666219187201334 0 +0.08480958662713545 0.2835410537880573 0 +0.08694365465394893 0.0151645126294419 0 +0.08124589756615076 0.02988881942009314 0 +0.07734512704694452 0.0444379884312248 0 +0.07453312982524721 0.05896206868972203 0 +0.07249503022380939 0.07356672081534776 0 +0.07105930644173161 0.08833432160138008 0 +0.07011722679187712 0.1033299809986393 0 +0.06959113423803064 0.1186008555110026 0 +0.06942115262374371 0.1341729368871379 0 +0.06956027548823038 0.1500488910740899 0 +0.06997458914660702 0.1662096074260965 0 +0.0706482801298126 0.1826200340395445 0 +0.07159440801581606 0.1992371891607003 0 +0.07287332148999685 0.2160164695139022 0 +0.07462234065867968 0.2329128306451688 0 +0.07710365190329735 0.249873697340604 0 +0.08078857338859603 0.2668206877186082 0 +0.08656488676360412 0.2836166260488554 0 +0.0884657237090324 0.01523371052675922 0 +0.0833142886284326 0.03006687161326848 0 +0.07980929587726822 0.04474276743763095 0 +0.07729731682942742 0.059398796341171 0 +0.07548753799887847 0.07413100786004587 0 +0.07422132251162304 0.08901339081383318 0 +0.07339779694559891 0.1041032878954707 0 +0.07294483573203984 0.119441026058213 0 +0.07280709169265359 0.1350476311984484 0 +0.07294164274850373 0.1509235943115593 0 +0.07331813379610687 0.1670506833000172 0 +0.07392304183839027 0.1833971406988227 0 +0.07476902207345099 0.1999246071931738 0 +0.07591108818908733 0.2165935880297158 0 +0.07747292523963477 0.2333644416805635 0 +0.07969001057149373 0.2501911604576898 0 +0.08298760281096558 0.2670053166459223 0 +0.08817553339213485 0.2836862898872534 0 +0.08987834740677571 0.01529877880446726 0 +0.08526492757253148 0.03023652266567063 0 +0.08214508891539281 0.04503434374973565 0 +0.07992193702271759 0.05981683748330902 0 +0.07832962577039361 0.07467047730079274 0 +0.07722304056679219 0.08966111040875717 0 +0.07650969867137454 0.104838820685452 0 +0.0761234449127343 0.1202378273154894 0 +0.07601376143338434 0.1358749420540018 0 +0.07614190967814401 0.1517490209051104 0 +0.07648096840507733 0.1678428940917961 0 +0.07701941793817925 0.1841279752315086 0 +0.07776919077029133 0.2005702523762469 0 +0.07877978632105631 0.2171350169142498 0 +0.0801614346111192 0.233787684529341 0 +0.08212366353770946 0.2504883681953621 0 +0.08504684734541625 0.2671778616989354 0 +0.08966299433166734 0.2837509827314612 0 +0.0911988974502224 0.01536037591270382 0 +0.08711519928738337 0.03039909624012595 0 +0.0843698308653692 0.04531470390252633 0 +0.08242435828879253 0.06021886021388595 0 +0.08103885575614725 0.07518852472797335 0 +0.08008233126836738 0.09028165333859252 0 +0.07947115927746957 0.1055415367249938 0 +0.0791455181877517 0.1209969139355373 0 +0.07905994552626244 0.1366610161347093 0 +0.07917993535458176 0.1525315153601081 0 +0.07948187505129073 0.168592470779113 0 +0.0799560174860858 0.1848183913653676 0 +0.08061335265733641 0.2011794043359312 0 +0.08149780193359793 0.2176453175678687 0 +0.08270640933227187 0.2341862856151269 0 +0.08442348679766574 0.2507681078345605 0 +0.08698542671570349 0.267340099698603 0 +0.0910460025412464 0.2838114975758676 0 +0.09244265580976904 0.01541908940641182 0 +0.08888092199024068 0.03055580753726568 0 +0.08649977973025716 0.04558571333690979 0 +0.08482131900222842 0.06060740857396933 0 +0.08363251134957322 0.07568840943691157 0 +0.08281704756155936 0.09087901946994238 0 +0.08230056673025325 0.1062161548847307 0 +0.08202988124113403 0.1217236196289866 0 +0.08196475356201545 0.1374116032266708 0 +0.08207492642166385 0.1532769784082377 0 +0.08233997696960124 0.1693051902301176 0 +0.08275174627904701 0.1854738053513247 0 +0.08332013306445513 0.2017569334882289 0 +0.08408348980345794 0.2181286714992752 0 +0.08512598352215313 0.2345636223465755 0 +0.08660738982004798 0.2510328625992306 0 +0.08882080347840333 0.2674935797239399 0 +0.09234099997677175 0.283868510678027 0 +0.09362319748455053 0.01547544987570303 0 +0.09057666527053332 0.03070779125241929 0 +0.08855039546085743 0.04584915859068193 0 +0.08712914500462392 0.06098495650538261 0 +0.08612775932690055 0.07617331710846897 0 +0.08544513336710126 0.09145710560282136 0 +0.08501653069884346 0.1068672334045717 0 +0.08479565132878926 0.1224230432748323 0 +0.08474761398927205 0.138132148363084 0 +0.08484640864155397 0.1539909597904677 0 +0.08507469532297987 0.1699864626221657 0 +0.08542575626771072 0.1860992745280097 0 +0.08590829757920145 0.2023073706329721 0 +0.0865551501718969 0.218588943289847 0 +0.08743791496620225 0.2349227817346361 0 +0.08869243827814163 0.2512848648643959 0 +0.09056903300481302 0.2676396647955613 0 +0.09356252373987986 0.2839226038523612 0 +0.09475272128604816 0.01552994284122036 0 +0.09221601910963374 0.0308561246877655 0 +0.09053655866790813 0.04610678127828931 0 +0.08936391597216028 0.06135395031208159 0 +0.08854176331449165 0.07664640815885532 0 +0.08798468810754345 0.0920197587054884 0 +0.08763790567964326 0.1074992286573757 0 +0.08746222829562576 0.1231001138614564 0 +0.0874282420168952 0.1388278615729913 0 +0.08751417904886083 0.154678728958396 0 +0.08770569175340363 0.1706413985701762 0 +0.08799738550258029 0.186699558106961 0 +0.08839670011813645 0.2028349612774454 0 +0.08893100242337112 0.2190297302748087 0 +0.08965961062653993 0.2352666071053958 0 +0.09069496386783799 0.2515261388550351 0 +0.09224498437958345 0.2677795655436702 0 +0.09472354268948988 0.2839742826362911 0 +0.09584233979648693 0.01558301917084371 0 +0.09381182118178773 0.03100184707539813 0 +0.09247274539857098 0.04636030540346552 0 +0.09154158742016875 0.06171684166145426 0 +0.09089175539287471 0.07711085429036242 0 +0.09045399408262299 0.09257081539795611 0 +0.09018378266750759 0.1081165383859687 0 +0.09004925976021085 0.1237596383373668 0 +0.09002658664555897 0.1395037695943434 0 +0.09009824143085227 0.1553453278459627 0 +0.09025279779832551 0.1712748594676937 0 +0.0904860878609499 0.187279162911964 0 +0.09080422279671835 0.2033437067777373 0 +0.09122915085954528 0.2194544006457369 0 +0.09180814141550189 0.2355977345330429 0 +0.09263065613881426 0.2517585345454954 0 +0.09386253212831394 0.267914367451782 0 +0.09583575294513459 0.2840239913038722 0 +0.0969023386366806 0.01563510447057158 0 +0.09537634948688495 0.0311459760129735 0 +0.094373164554477 0.04661145943840019 0 +0.093678075255792 0.06207611309970934 0 +0.0931950730002027 0.07756986569800234 0 +0.09287151354027108 0.09311413053857694 0 +0.09267345537175012 0.1087235326186742 0 +0.09257658533959093 0.1244063358728891 0 +0.09256276050072069 0.1401647531750084 0 +0.092618727599439 0.1559956091518722 0 +0.09273593242467318 0.1718914940356769 0 +0.09291135257429124 0.1878423765536176 0 +0.09314970638197111 0.2038373943382644 0 +0.09346754029410709 0.2198661217307071 0 +0.09390024361680263 0.2359186206363749 0 +0.09451463502763785 0.2519837543420853 0 +0.09543472109719714 0.2680450530028318 0 +0.09690984076090155 0.284072125484694 0 +0.0979424146482243 0.01568660784160181 0 +0.09692148725479312 0.03128952183716663 0 +0.09625186510901963 0.0468619945820156 0 +0.09578931032668735 0.0624342980661703 0 +0.09546916844623506 0.07802671136760722 0 +0.09525586206257894 0.09365359765599737 0 +0.09512636654283001 0.1093245752155486 0 +0.09506416432281332 0.1250448616231402 0 +0.0950569557809994 0.1408155730668006 0 +0.09509580679315333 0.1566342630709084 0 +0.09517500916174201 0.1724957637094112 0 +0.09529261428983128 0.1883932903266647 0 +0.09545187145001566 0.2043196178513181 0 +0.0956639012671482 0.2202678802128395 0 +0.09595230706219648 0.2362315633616371 0 +0.09636150556640277 0.2522033740139941 0 +0.09697390843102592 0.2681725198499242 0 +0.09795572130203035 0.2841190429886794 0 +0.0989719025361445 0.01573793036658277 0 +0.09845886695523882 0.03143350073210717 0 +0.09812282076373513 0.04711370064251993 0 +0.0978912707583601 0.06279399744163996 0 +0.09773159833260397 0.07848473491226286 0 +0.09762576458540817 0.09419316389632564 0 +0.09756203955795306 0.1099240388128554 0 +0.09753199081115524 0.1256798227809212 0 +0.09752934937492716 0.141460887440157 0 +0.0975495855010083 0.1572658350470997 0 +0.09758983456598229 0.1730919592219104 0 +0.09764915503456141 0.1889358139599983 0 +0.09772923147780814 0.2047937914905861 0 +0.09783568631510411 0.2206624960213323 0 +0.09798035223614308 0.2365387173298136 0 +0.09818539825846269 0.2524188592292663 0 +0.09849188761160503 0.2682975959713053 0 +0.09898276204372107 0.2841650733241296 0 +2 15 0 64 +638 +639 +640 +641 +642 +643 +644 +645 +646 +647 +648 +649 +650 +651 +652 +653 +654 +655 +656 +657 +658 +659 +660 +661 +662 +663 +664 +665 +666 +667 +668 +669 +670 +671 +672 +673 +674 +675 +676 +677 +678 +679 +680 +681 +682 +683 +684 +685 +686 +687 +688 +689 +690 +691 +692 +693 +694 +695 +696 +697 +698 +699 +700 +701 +0.2059961831449814 0.1888888888889119 0 +0.2119988867644767 0.1888888888889087 0 +0.2180155380751619 0.1888888888889049 0 +0.2240556741122001 0.1888888888889008 0 +0.2301328907141563 0.1888888888888968 0 +0.2362683876336542 0.1888888888888928 0 +0.24249787851132 0.1888888888888885 0 +0.2488857471250591 0.1888888888888839 0 +0.2064873066762567 0.1777777777778285 0 +0.2129887407492314 0.1777777777778178 0 +0.2195200819696328 0.1777777777778079 0 +0.2261007501890784 0.1777777777777987 0 +0.2327565152148688 0.1777777777777901 0 +0.2395233584456571 0.1777777777777818 0 +0.2464527348732792 0.1777777777777734 0 +0.2536172870746236 0.177777777777765 0 +0.2070179342666262 0.1666666666667404 0 +0.2140569067281896 0.1666666666667223 0 +0.2211398938273211 0.1666666666667063 0 +0.2282939388956915 0.1666666666666921 0 +0.235552606155097 0.166666666666679 0 +0.2429587071945336 0.1666666666666664 0 +0.2505667826308733 0.1666666666666536 0 +0.258443919548154 0.1666666666666399 0 +0.2075815014094561 0.1555555555556453 0 +0.215189339096983 0.1555555555556191 0 +0.2228517131450438 0.1555555555555981 0 +0.230600581016662 0.1555555555555802 0 +0.2384735335695898 0.1555555555555638 0 +0.2465154298015984 0.1555555555555475 0 +0.2547793164192706 0.1555555555555294 0 +0.2633260102171202 0.155555555555507 0 +0.2081773687301962 0.1444444444445295 0 +0.2163842812841149 0.1444444444445036 0 +0.2246518532201926 0.1444444444444824 0 +0.2330142603332175 0.1444444444444639 0 +0.2415100128052106 0.1444444444444468 0 +0.2501828051141725 0.1444444444444295 0 +0.259081708428393 0.1444444444444103 0 +0.2682606599830627 0.1444444444443878 0 +0.2088111462310094 0.133333333333401 0 +0.2176524152486283 0.13333333333338 0 +0.2265550673137395 0.133333333333362 0 +0.2355524917053182 0.1333333333333458 0 +0.2446809032659029 0.1333333333333306 0 +0.2539796153537501 0.1333333333333151 0 +0.2634908781313161 0.1333333333332985 0 +0.2732595522920556 0.1333333333332807 0 +0.2094953242014855 0.1222222222222679 0 +0.2190177101365892 0.1222222222222532 0 +0.2285948167035521 0.1222222222222401 0 +0.2382553371425895 0.1222222222222281 0 +0.2480291575116575 0.1222222222222166 0 +0.2579472146712042 0.122222222222205 0 +0.2680412154162276 0.1222222222221926 0 +0.2783436510576747 0.122222222222179 0 +0.2102509571807479 0.1111111111111339 0 +0.2205202711324805 0.111111111111126 0 +0.2308263554271822 0.111111111111119 0 +0.2411876294313529 0.1111111111111124 0 +0.2516223151636553 0.1111111111111062 0 +0.2621481529094727 0.1111111111110996 0 +0.2727822613694358 0.1111111111110925 0 +0.2835415566893428 0.1111111111110842 0 +$EndNodes +$Elements +23 856 1 856 +0 1 15 1 +1 1 +0 2 15 1 +2 2 +0 3 15 1 +3 3 +0 4 15 1 +4 4 +0 5 15 1 +5 5 +0 6 15 1 +6 6 +0 7 15 1 +7 7 +0 8 15 1 +8 8 +0 9 15 1 +9 9 +0 10 15 1 +10 10 +0 11 15 1 +11 11 +1 1 1 29 +12 1 12 +13 12 13 +14 13 14 +15 14 15 +16 15 16 +17 16 17 +18 17 18 +19 18 19 +20 19 20 +21 20 21 +22 21 22 +23 22 23 +24 23 24 +25 24 25 +26 25 26 +27 26 27 +28 27 28 +29 28 29 +30 29 30 +31 30 31 +32 31 32 +33 32 33 +34 33 34 +35 34 35 +36 35 36 +37 36 37 +38 37 38 +39 38 39 +40 39 2 +1 2 1 19 +41 3 40 +42 40 41 +43 41 42 +44 42 43 +45 43 44 +46 44 45 +47 45 46 +48 46 47 +49 47 48 +50 48 49 +51 49 50 +52 50 51 +53 51 52 +54 52 53 +55 53 54 +56 54 55 +57 55 56 +58 56 57 +59 57 2 +1 3 1 29 +60 3 58 +61 58 59 +62 59 60 +63 60 61 +64 61 62 +65 62 63 +66 63 64 +67 64 65 +68 65 66 +69 66 67 +70 67 68 +71 68 69 +72 69 70 +73 70 71 +74 71 72 +75 72 73 +76 73 74 +77 74 75 +78 75 76 +79 76 77 +80 77 78 +81 78 79 +82 79 80 +83 80 81 +84 81 82 +85 82 83 +86 83 84 +87 84 85 +88 85 4 +1 4 1 5 +89 1 86 +90 86 87 +91 87 88 +92 88 89 +93 89 5 +1 5 1 5 +94 5 90 +95 90 91 +96 91 92 +97 92 93 +98 93 6 +1 6 1 9 +99 6 94 +100 94 95 +101 95 96 +102 96 97 +103 97 98 +104 98 99 +105 99 100 +106 100 101 +107 101 4 +1 10 1 9 +108 8 102 +109 102 103 +110 103 104 +111 104 105 +112 105 106 +113 106 107 +114 107 108 +115 108 109 +116 109 11 +1 11 1 9 +117 11 110 +118 110 111 +119 111 112 +120 112 113 +121 113 114 +122 114 115 +123 115 116 +124 116 117 +125 117 10 +1 12 1 9 +126 10 118 +127 118 119 +128 119 120 +129 120 121 +130 121 122 +131 122 123 +132 123 124 +133 124 125 +134 125 7 +1 13 1 9 +135 7 126 +136 126 127 +137 127 128 +138 128 129 +139 129 130 +140 130 131 +141 131 132 +142 132 133 +143 133 8 +2 1 3 551 +144 1 12 134 86 +145 86 134 135 87 +146 87 135 136 88 +147 88 136 137 89 +148 89 137 138 5 +149 5 138 139 90 +150 90 139 140 91 +151 91 140 141 92 +152 92 141 142 93 +153 93 142 143 6 +154 6 143 144 94 +155 94 144 145 95 +156 95 145 146 96 +157 96 146 147 97 +158 97 147 148 98 +159 98 148 149 99 +160 99 149 150 100 +161 100 150 151 101 +162 101 151 85 4 +163 12 13 152 134 +164 134 152 153 135 +165 135 153 154 136 +166 136 154 155 137 +167 137 155 156 138 +168 138 156 157 139 +169 139 157 158 140 +170 140 158 159 141 +171 141 159 160 142 +172 142 160 161 143 +173 143 161 162 144 +174 144 162 163 145 +175 145 163 164 146 +176 146 164 165 147 +177 147 165 166 148 +178 148 166 167 149 +179 149 167 168 150 +180 150 168 169 151 +181 151 169 84 85 +182 13 14 170 152 +183 152 170 171 153 +184 153 171 172 154 +185 154 172 173 155 +186 155 173 174 156 +187 156 174 175 157 +188 157 175 176 158 +189 158 176 177 159 +190 159 177 178 160 +191 160 178 179 161 +192 161 179 180 162 +193 162 180 181 163 +194 163 181 182 164 +195 164 182 183 165 +196 165 183 184 166 +197 166 184 185 167 +198 167 185 186 168 +199 168 186 187 169 +200 169 187 83 84 +201 14 15 188 170 +202 170 188 189 171 +203 171 189 190 172 +204 172 190 191 173 +205 173 191 192 174 +206 174 192 193 175 +207 175 193 194 176 +208 176 194 195 177 +209 177 195 196 178 +210 178 196 197 179 +211 179 197 198 180 +212 180 198 199 181 +213 181 199 200 182 +214 182 200 201 183 +215 183 201 202 184 +216 184 202 203 185 +217 185 203 204 186 +218 186 204 205 187 +219 187 205 82 83 +220 15 16 206 188 +221 188 206 207 189 +222 189 207 208 190 +223 190 208 209 191 +224 191 209 210 192 +225 192 210 211 193 +226 193 211 212 194 +227 194 212 213 195 +228 195 213 214 196 +229 196 214 215 197 +230 197 215 216 198 +231 198 216 217 199 +232 199 217 218 200 +233 200 218 219 201 +234 201 219 220 202 +235 202 220 221 203 +236 203 221 222 204 +237 204 222 223 205 +238 205 223 81 82 +239 16 17 224 206 +240 206 224 225 207 +241 207 225 226 208 +242 208 226 227 209 +243 209 227 228 210 +244 210 228 229 211 +245 211 229 230 212 +246 212 230 231 213 +247 213 231 232 214 +248 214 232 233 215 +249 215 233 234 216 +250 216 234 235 217 +251 217 235 236 218 +252 218 236 237 219 +253 219 237 238 220 +254 220 238 239 221 +255 221 239 240 222 +256 222 240 241 223 +257 223 241 80 81 +258 17 18 242 224 +259 224 242 243 225 +260 225 243 244 226 +261 226 244 245 227 +262 227 245 246 228 +263 228 246 247 229 +264 229 247 248 230 +265 230 248 249 231 +266 231 249 250 232 +267 232 250 251 233 +268 233 251 252 234 +269 234 252 253 235 +270 235 253 254 236 +271 236 254 255 237 +272 237 255 256 238 +273 238 256 257 239 +274 239 257 258 240 +275 240 258 259 241 +276 241 259 79 80 +277 18 19 260 242 +278 242 260 261 243 +279 243 261 262 244 +280 244 262 263 245 +281 245 263 264 246 +282 246 264 265 247 +283 247 265 266 248 +284 248 266 267 249 +285 249 267 268 250 +286 250 268 269 251 +287 251 269 270 252 +288 252 270 271 253 +289 253 271 272 254 +290 254 272 273 255 +291 255 273 274 256 +292 256 274 275 257 +293 257 275 276 258 +294 258 276 277 259 +295 259 277 78 79 +296 19 20 278 260 +297 260 278 279 261 +298 261 279 280 262 +299 262 280 281 263 +300 263 281 282 264 +301 264 282 283 265 +302 265 283 284 266 +303 266 284 285 267 +304 267 285 286 268 +305 268 286 287 269 +306 269 287 288 270 +307 270 288 289 271 +308 271 289 290 272 +309 272 290 291 273 +310 273 291 292 274 +311 274 292 293 275 +312 275 293 294 276 +313 276 294 295 277 +314 277 295 77 78 +315 20 21 296 278 +316 278 296 297 279 +317 279 297 298 280 +318 280 298 299 281 +319 281 299 300 282 +320 282 300 301 283 +321 283 301 302 284 +322 284 302 303 285 +323 285 303 304 286 +324 286 304 305 287 +325 287 305 306 288 +326 288 306 307 289 +327 289 307 308 290 +328 290 308 309 291 +329 291 309 310 292 +330 292 310 311 293 +331 293 311 312 294 +332 294 312 313 295 +333 295 313 76 77 +334 21 22 314 296 +335 296 314 315 297 +336 297 315 316 298 +337 298 316 317 299 +338 299 317 318 300 +339 300 318 319 301 +340 301 319 320 302 +341 302 320 321 303 +342 303 321 322 304 +343 304 322 323 305 +344 305 323 324 306 +345 306 324 325 307 +346 307 325 326 308 +347 308 326 327 309 +348 309 327 328 310 +349 310 328 329 311 +350 311 329 330 312 +351 312 330 331 313 +352 313 331 75 76 +353 22 23 332 314 +354 314 332 333 315 +355 315 333 334 316 +356 316 334 335 317 +357 317 335 336 318 +358 318 336 337 319 +359 319 337 338 320 +360 320 338 339 321 +361 321 339 340 322 +362 322 340 341 323 +363 323 341 342 324 +364 324 342 343 325 +365 325 343 344 326 +366 326 344 345 327 +367 327 345 346 328 +368 328 346 347 329 +369 329 347 348 330 +370 330 348 349 331 +371 331 349 74 75 +372 23 24 350 332 +373 332 350 351 333 +374 333 351 352 334 +375 334 352 353 335 +376 335 353 354 336 +377 336 354 355 337 +378 337 355 356 338 +379 338 356 357 339 +380 339 357 358 340 +381 340 358 359 341 +382 341 359 360 342 +383 342 360 361 343 +384 343 361 362 344 +385 344 362 363 345 +386 345 363 364 346 +387 346 364 365 347 +388 347 365 366 348 +389 348 366 367 349 +390 349 367 73 74 +391 24 25 368 350 +392 350 368 369 351 +393 351 369 370 352 +394 352 370 371 353 +395 353 371 372 354 +396 354 372 373 355 +397 355 373 374 356 +398 356 374 375 357 +399 357 375 376 358 +400 358 376 377 359 +401 359 377 378 360 +402 360 378 379 361 +403 361 379 380 362 +404 362 380 381 363 +405 363 381 382 364 +406 364 382 383 365 +407 365 383 384 366 +408 366 384 385 367 +409 367 385 72 73 +410 25 26 386 368 +411 368 386 387 369 +412 369 387 388 370 +413 370 388 389 371 +414 371 389 390 372 +415 372 390 391 373 +416 373 391 392 374 +417 374 392 393 375 +418 375 393 394 376 +419 376 394 395 377 +420 377 395 396 378 +421 378 396 397 379 +422 379 397 398 380 +423 380 398 399 381 +424 381 399 400 382 +425 382 400 401 383 +426 383 401 402 384 +427 384 402 403 385 +428 385 403 71 72 +429 26 27 404 386 +430 386 404 405 387 +431 387 405 406 388 +432 388 406 407 389 +433 389 407 408 390 +434 390 408 409 391 +435 391 409 410 392 +436 392 410 411 393 +437 393 411 412 394 +438 394 412 413 395 +439 395 413 414 396 +440 396 414 415 397 +441 397 415 416 398 +442 398 416 417 399 +443 399 417 418 400 +444 400 418 419 401 +445 401 419 420 402 +446 402 420 421 403 +447 403 421 70 71 +448 27 28 422 404 +449 404 422 423 405 +450 405 423 424 406 +451 406 424 425 407 +452 407 425 426 408 +453 408 426 427 409 +454 409 427 428 410 +455 410 428 429 411 +456 411 429 430 412 +457 412 430 431 413 +458 413 431 432 414 +459 414 432 433 415 +460 415 433 434 416 +461 416 434 435 417 +462 417 435 436 418 +463 418 436 437 419 +464 419 437 438 420 +465 420 438 439 421 +466 421 439 69 70 +467 28 29 440 422 +468 422 440 441 423 +469 423 441 442 424 +470 424 442 443 425 +471 425 443 444 426 +472 426 444 445 427 +473 427 445 446 428 +474 428 446 447 429 +475 429 447 448 430 +476 430 448 449 431 +477 431 449 450 432 +478 432 450 451 433 +479 433 451 452 434 +480 434 452 453 435 +481 435 453 454 436 +482 436 454 455 437 +483 437 455 456 438 +484 438 456 457 439 +485 439 457 68 69 +486 29 30 458 440 +487 440 458 459 441 +488 441 459 460 442 +489 442 460 461 443 +490 443 461 462 444 +491 444 462 463 445 +492 445 463 464 446 +493 446 464 465 447 +494 447 465 466 448 +495 448 466 467 449 +496 449 467 468 450 +497 450 468 469 451 +498 451 469 470 452 +499 452 470 471 453 +500 453 471 472 454 +501 454 472 473 455 +502 455 473 474 456 +503 456 474 475 457 +504 457 475 67 68 +505 30 31 476 458 +506 458 476 477 459 +507 459 477 478 460 +508 460 478 479 461 +509 461 479 480 462 +510 462 480 481 463 +511 463 481 482 464 +512 464 482 483 465 +513 465 483 484 466 +514 466 484 485 467 +515 467 485 486 468 +516 468 486 487 469 +517 469 487 488 470 +518 470 488 489 471 +519 471 489 490 472 +520 472 490 491 473 +521 473 491 492 474 +522 474 492 493 475 +523 475 493 66 67 +524 31 32 494 476 +525 476 494 495 477 +526 477 495 496 478 +527 478 496 497 479 +528 479 497 498 480 +529 480 498 499 481 +530 481 499 500 482 +531 482 500 501 483 +532 483 501 502 484 +533 484 502 503 485 +534 485 503 504 486 +535 486 504 505 487 +536 487 505 506 488 +537 488 506 507 489 +538 489 507 508 490 +539 490 508 509 491 +540 491 509 510 492 +541 492 510 511 493 +542 493 511 65 66 +543 32 33 512 494 +544 494 512 513 495 +545 495 513 514 496 +546 496 514 515 497 +547 497 515 516 498 +548 498 516 517 499 +549 499 517 518 500 +550 500 518 519 501 +551 501 519 520 502 +552 502 520 521 503 +553 503 521 522 504 +554 504 522 523 505 +555 505 523 524 506 +556 506 524 525 507 +557 507 525 526 508 +558 508 526 527 509 +559 509 527 528 510 +560 510 528 529 511 +561 511 529 64 65 +562 33 34 530 512 +563 512 530 531 513 +564 513 531 532 514 +565 514 532 533 515 +566 515 533 534 516 +567 516 534 535 517 +568 517 535 536 518 +569 518 536 537 519 +570 519 537 538 520 +571 520 538 539 521 +572 521 539 540 522 +573 522 540 541 523 +574 523 541 542 524 +575 524 542 543 525 +576 525 543 544 526 +577 526 544 545 527 +578 527 545 546 528 +579 528 546 547 529 +580 529 547 63 64 +581 34 35 548 530 +582 530 548 549 531 +583 531 549 550 532 +584 532 550 551 533 +585 533 551 552 534 +586 534 552 553 535 +587 535 553 554 536 +588 536 554 555 537 +589 537 555 556 538 +590 538 556 557 539 +591 539 557 558 540 +592 540 558 559 541 +593 541 559 560 542 +594 542 560 561 543 +595 543 561 562 544 +596 544 562 563 545 +597 545 563 564 546 +598 546 564 565 547 +599 547 565 62 63 +600 35 36 566 548 +601 548 566 567 549 +602 549 567 568 550 +603 550 568 569 551 +604 551 569 570 552 +605 552 570 571 553 +606 553 571 572 554 +607 554 572 573 555 +608 555 573 574 556 +609 556 574 575 557 +610 557 575 576 558 +611 558 576 577 559 +612 559 577 578 560 +613 560 578 579 561 +614 561 579 580 562 +615 562 580 581 563 +616 563 581 582 564 +617 564 582 583 565 +618 565 583 61 62 +619 36 37 584 566 +620 566 584 585 567 +621 567 585 586 568 +622 568 586 587 569 +623 569 587 588 570 +624 570 588 589 571 +625 571 589 590 572 +626 572 590 591 573 +627 573 591 592 574 +628 574 592 593 575 +629 575 593 594 576 +630 576 594 595 577 +631 577 595 596 578 +632 578 596 597 579 +633 579 597 598 580 +634 580 598 599 581 +635 581 599 600 582 +636 582 600 601 583 +637 583 601 60 61 +638 37 38 602 584 +639 584 602 603 585 +640 585 603 604 586 +641 586 604 605 587 +642 587 605 606 588 +643 588 606 607 589 +644 589 607 608 590 +645 590 608 609 591 +646 591 609 610 592 +647 592 610 611 593 +648 593 611 612 594 +649 594 612 613 595 +650 595 613 614 596 +651 596 614 615 597 +652 597 615 616 598 +653 598 616 617 599 +654 599 617 618 600 +655 600 618 619 601 +656 601 619 59 60 +657 38 39 620 602 +658 602 620 621 603 +659 603 621 622 604 +660 604 622 623 605 +661 605 623 624 606 +662 606 624 625 607 +663 607 625 626 608 +664 608 626 627 609 +665 609 627 628 610 +666 610 628 629 611 +667 611 629 630 612 +668 612 630 631 613 +669 613 631 632 614 +670 614 632 633 615 +671 615 633 634 616 +672 616 634 635 617 +673 617 635 636 618 +674 618 636 637 619 +675 619 637 58 59 +676 39 2 57 620 +677 620 57 56 621 +678 621 56 55 622 +679 622 55 54 623 +680 623 54 53 624 +681 624 53 52 625 +682 625 52 51 626 +683 626 51 50 627 +684 627 50 49 628 +685 628 49 48 629 +686 629 48 47 630 +687 630 47 46 631 +688 631 46 45 632 +689 632 45 44 633 +690 633 44 43 634 +691 634 43 42 635 +692 635 42 41 636 +693 636 41 40 637 +694 637 40 3 58 +2 15 2 162 +695 7 126 125 +696 125 126 638 +697 125 638 124 +698 124 638 639 +699 124 639 123 +700 123 639 640 +701 123 640 122 +702 122 640 641 +703 122 641 121 +704 121 641 642 +705 121 642 120 +706 120 642 643 +707 120 643 119 +708 119 643 644 +709 119 644 118 +710 118 644 645 +711 118 645 10 +712 10 645 117 +713 126 127 638 +714 638 127 646 +715 638 646 639 +716 639 646 647 +717 639 647 640 +718 640 647 648 +719 640 648 641 +720 641 648 649 +721 641 649 642 +722 642 649 650 +723 642 650 643 +724 643 650 651 +725 643 651 644 +726 644 651 652 +727 644 652 645 +728 645 652 653 +729 645 653 117 +730 117 653 116 +731 127 128 646 +732 646 128 654 +733 646 654 647 +734 647 654 655 +735 647 655 648 +736 648 655 656 +737 648 656 649 +738 649 656 657 +739 649 657 650 +740 650 657 658 +741 650 658 651 +742 651 658 659 +743 651 659 652 +744 652 659 660 +745 652 660 653 +746 653 660 661 +747 653 661 116 +748 116 661 115 +749 128 129 654 +750 654 129 662 +751 654 662 655 +752 655 662 663 +753 655 663 656 +754 656 663 664 +755 656 664 657 +756 657 664 665 +757 657 665 658 +758 658 665 666 +759 658 666 659 +760 659 666 667 +761 659 667 660 +762 660 667 668 +763 660 668 661 +764 661 668 669 +765 661 669 115 +766 115 669 114 +767 129 130 662 +768 662 130 670 +769 662 670 663 +770 663 670 671 +771 663 671 664 +772 664 671 672 +773 664 672 665 +774 665 672 673 +775 665 673 666 +776 666 673 674 +777 666 674 667 +778 667 674 675 +779 667 675 668 +780 668 675 676 +781 668 676 669 +782 669 676 677 +783 669 677 114 +784 114 677 113 +785 130 131 670 +786 670 131 678 +787 670 678 671 +788 671 678 679 +789 671 679 672 +790 672 679 680 +791 672 680 673 +792 673 680 681 +793 673 681 674 +794 674 681 682 +795 674 682 675 +796 675 682 683 +797 675 683 676 +798 676 683 684 +799 676 684 677 +800 677 684 685 +801 677 685 113 +802 113 685 112 +803 131 132 678 +804 678 132 686 +805 678 686 679 +806 679 686 687 +807 679 687 680 +808 680 687 688 +809 680 688 681 +810 681 688 689 +811 681 689 682 +812 682 689 690 +813 682 690 683 +814 683 690 691 +815 683 691 684 +816 684 691 692 +817 684 692 685 +818 685 692 693 +819 685 693 112 +820 112 693 111 +821 132 133 686 +822 686 133 694 +823 686 694 687 +824 687 694 695 +825 687 695 688 +826 688 695 696 +827 688 696 689 +828 689 696 697 +829 689 697 690 +830 690 697 698 +831 690 698 691 +832 691 698 699 +833 691 699 692 +834 692 699 700 +835 692 700 693 +836 693 700 701 +837 693 701 111 +838 111 701 110 +839 133 8 694 +840 694 8 102 +841 694 102 695 +842 695 102 103 +843 695 103 696 +844 696 103 104 +845 696 104 697 +846 697 104 105 +847 697 105 698 +848 698 105 106 +849 698 106 699 +850 699 106 107 +851 699 107 700 +852 700 107 108 +853 700 108 701 +854 701 108 109 +855 701 109 110 +856 110 109 11 +$EndElements diff --git a/test/runtests.jl b/test/runtests.jl index b28d3664..4ec308e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,540 +1,345 @@ -ENV["MPLBACKEND"]="agg" - +ENV["MPLBACKEND"] = "agg" using Test, ExtendableGrids, GridVisualize, SHA, SimplexGridFactory, Triangulate import StatsBase -using Gmsh: gmsh using StatsBase: countmap import CairoMakie +CairoMakie.activate!(; type = "svg", visible = false) -CairoMakie.activate!(type="svg",visible=false) - - -### Tests for the gmsh extension: - - - -function multidimsort(A) - i1 = sortperm(A[1,:]) - A1 = A[:,i1] - for j=2:size(A,1) - cm = countmap(A1[j-1,:]) - for (key,val) in cm - if val > 1 #if there only is one entry with this key, the reordering is not necessary - inds = findall(z->z==key, A1[j-1,:]) #indices of that key in A1 - sorted_inds = sortperm(A1[j, inds]) - A1[:,inds] = A1[:,inds[sorted_inds]] #reorder - end - end - end - return A1 - -end - -#= -Confidence level: -- :low : Point numbers etc are the same -- :full : all arrays are equal (besides the coordinate array, the arrays only have to be equal up to permutations) -=# - - -function seemingly_equal_with_sort(grid1::ExtendableGrid, grid2::ExtendableGrid; confidence=:full) - if confidence==:full - for key in keys(grid1) - #@warn key - if !haskey(grid2,key) - return false - end - - if key == Coordinates - #@warn "now checking coords" - if !seemingly_equal(grid1[key],grid2[key]) - return false - end - else - try - #@warn "1" - sorted = false - s1 = size(grid1[key]) - s2 = size(grid2[key]) - - - #@warn "2" - - if length(s1) == 1 - #@warn "l(s) = 1" - ind1 = sortperm(grid1[key]) - ind2 = sortperm(grid2[key]) - sa1 = grid1[key][ind1] - sa2 = grid2[key][ind2] - else - #@warn "l(s) > 1" - sa1 = multidimsort(grid1[key]) - sa2 = multidimsort(grid2[key]) - - end - - - if !seemingly_equal(sa1,sa2) - return false - end - - catch - #@warn "doesn't have a size" #"not a sortable array") - if !seemingly_equal(grid1[key],grid2[key]) - return false - end - end - end - - end - return true - elseif confidence==:low - grid1_data=(num_nodes(grid1),num_cells(grid1), num_bfaces(grid1)) - grid2_data=(num_nodes(grid2),num_cells(grid2), num_bfaces(grid2)) - return grid1_data==grid2_data - else - error("Confidence level $(confidence) not implemented") - return false - end -end - - - - - -@testset "Read/write gmsh simplexgrid" begin - - gmsh.initialize() - - path = "" #"/home/johannes/Nextcloud/Documents/Uni/VIII/WIAS/Julia Packages (tries)/ExtendableGrids.jl_moving_to_ext/test/" - X = collect(0:0.02:2) - Y = collect(0:2:4) - grid1 = simplexgrid(X, Y) #ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh") - - ExtendableGrids.write_gmsh(path*"testfile.msh", grid1) - grid2 = ExtendableGrids.simplexgrid_from_gmsh(path*"testfile.msh") - gmsh.finalize() - - @test seemingly_equal_with_sort(grid2, grid1;confidence=:low) - @test seemingly_equal_with_sort(grid2, grid1;confidence=:full) - -end - -@testset "Read/write gmsh 2d" begin - - gmsh.initialize() - - path = "" - grid1 = ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh") - ExtendableGrids.write_gmsh(path*"testfile.msh", grid1) - grid2 = ExtendableGrids.simplexgrid_from_gmsh(path*"testfile.msh") - gmsh.finalize() - - @test seemingly_equal_with_sort(grid1, grid2;confidence=:low) - @test seemingly_equal_with_sort(grid1, grid2;confidence=:full) - +function testgrid(grid, testdata) + (num_nodes(grid), num_cells(grid), num_bfaces(grid)) == testdata end -@testset "Read/write gmsh 3d" begin - - gmsh.initialize() - - path = "" - grid1 = ExtendableGrids.simplexgrid_from_gmsh(path*"sto_3d.msh") - ExtendableGrids.write_gmsh(path*"testfile.msh", grid1) - grid2 = ExtendableGrids.simplexgrid_from_gmsh(path*"testfile.msh") - gmsh.finalize() - - @test seemingly_equal_with_sort(grid1, grid2;confidence=:low) - @test seemingly_equal_with_sort(grid1, grid2;confidence=:full) - -end - -@testset "gmsh neg test" begin - - gmsh.initialize() - - path = "" - grid1 = ExtendableGrids.simplexgrid_from_gmsh(path*"sto_2d.msh") - #grid2 = - #simplexgrid([0, 1, 2], [3, 4, 5]) - grid2 = ExtendableGrids.simplexgrid_from_gmsh(path*"sto_3d.msh") - gmsh.finalize() - - @test !seemingly_equal_with_sort(grid1, grid2;confidence=:low) - @test !seemingly_equal_with_sort(grid1, grid2;confidence=:full) - -end - - -### Original tests: - - @testset "Geomspace" begin function test_geomspace() - X1=geomspace(2.0,3.0,0.2,0.2) - X2=collect(2:0.2:3) - @assert X1 ≈ X2 - X1=geomspace(2.0,3.0,0.1,0.2) - X2=geomspace(2.0,3.0,0.2,0.1) - - - X1=geomspace(2.0,3.0,1.0e-5,0.2) - X2=geomspace(2.0,3.0,0.2,1.0e-5) - (X1[2]-X1[1]) ≈ (X2[end]-X2[end-1]) && - (X2[2]-X2[1]) ≈ (X1[end]-X1[end-1]) + X1 = geomspace(2.0, 3.0, 0.2, 0.2) + X2 = collect(2:0.2:3) + @assert X1 ≈ X2 + X1 = geomspace(2.0, 3.0, 0.1, 0.2) + X2 = geomspace(2.0, 3.0, 0.2, 0.1) + + X1 = geomspace(2.0, 3.0, 1.0e-5, 0.2) + X2 = geomspace(2.0, 3.0, 0.2, 1.0e-5) + (X1[2] - X1[1]) ≈ (X2[end] - X2[end - 1]) && + (X2[2] - X2[1]) ≈ (X1[end] - X1[end - 1]) end - function test_geomspace1(n,h0,h1) - for i=1:n - X=geomspace(0,1,h0*rand(),h1*rand()) + function test_geomspace1(n, h0, h1) + for i = 1:n + X = geomspace(0, 1, h0 * rand(), h1 * rand()) end true end - - @test test_geomspace() - @test test_geomspace1(100,0.01,0.01) - @test test_geomspace1(100,0.01,0.1) - @test test_geomspace1(100,0.1,0.01) - @test test_geomspace1(100,0.001,0.1) - - - @test length(geomspace(0,1,0.01,0.1))==27 - @test length(geomspace(0,1,0.001,0.1))==47 - @test length(geomspace(0,1,0.0001,0.1))==68 - @test length(geomspace(0,1,0.00001,0.1))==90 + @test test_geomspace() + @test test_geomspace1(100, 0.01, 0.01) + @test test_geomspace1(100, 0.01, 0.1) + @test test_geomspace1(100, 0.1, 0.01) + @test test_geomspace1(100, 0.001, 0.1) + + @test length(geomspace(0, 1, 0.01, 0.1)) == 27 + @test length(geomspace(0, 1, 0.001, 0.1)) == 47 + @test length(geomspace(0, 1, 0.0001, 0.1)) == 68 + @test length(geomspace(0, 1, 0.00001, 0.1)) == 90 end @testset "Basic" begin + function test_prepare_edges() + ## Compared with pdelib; have more of these + nodes = [0.0 1; + 1 0; + 0 1; + 1 1]' + cells = [1 2 3; + 2 3 4]' + cellmat = [1, 1] + bfaces = [1 2; + 1 3; + 2 4]' - - function test_prepare_edges() - ## Compared with pdelib; have more of these - nodes=[ - 0.0 1; - 1 0; - 0 1; - 1 1]' - - cells=[ - 1 2 3; - 2 3 4 - ]' - - cellmat=[1,1] - bfaces=[ 1 2; - 1 3; - 2 4]' - - bfacemat=[1,1] - - grid=simplexgrid(nodes,cells,cellmat,bfaces,bfacemat) - - grid[CellEdges]==Int32[1 2; 2 4; 3 5] && - grid[FaceNodes]==Int32[1 2 3 3 4; 2 3 1 4 2] && - grid[EdgeCells]==Int32[1 1 1 2 2; 0 2 0 0 0] + bfacemat = [1, 1] + grid = simplexgrid(nodes, cells, cellmat, bfaces, bfacemat) + + grid[CellEdges] == Int32[1 2; 2 4; 3 5] && + grid[FaceNodes] == Int32[1 2 3 3 4; 2 3 1 4 2] && + grid[EdgeCells] == Int32[1 1 1 2 2; 0 2 0 0 0] end @test test_prepare_edges() - @test let - g=simplexgrid(0:10) - bfc=g[BFaceCells] - bfc[1,1]==1 && bfc[1,2]==10 + g = simplexgrid(0:10) + bfc = g[BFaceCells] + bfc[1, 1] == 1 && bfc[1, 2] == 10 end @test let - g=simplexgrid(0:10) - bfn=g[BFaceNormals] - bfn==[-1.0 1.0] + g = simplexgrid(0:10) + bfn = g[BFaceNormals] + bfn == [-1.0 1.0] end - @test let - g=simplexgrid(0:2, 0:2) - bfn=g[BFaceNormals] - bfn==[0.0 -1.0; - -1.0 0.0; - 0.0 -1.0; - 1.0 0.0; - 0.0 1.0; - -1.0 0.0; - 1.0 0.0; - 0.0 1.0]' + g = simplexgrid(0:2, 0:2) + bfn = g[BFaceNormals] + bfn == [0.0 -1.0; + -1.0 0.0; + 0.0 -1.0; + 1.0 0.0; + 0.0 1.0; + -1.0 0.0; + 1.0 0.0; + 0.0 1.0]' end @test let - g=simplexgrid(0:2, 0:2) - bfc=g[BFaceCells] - tryfix(bfc)==[1 2 3 3 6 6 7 8;] + g = simplexgrid(0:2, 0:2) + bfc = g[BFaceCells] + tryfix(bfc) == [1 2 3 3 6 6 7 8;] end @test let - X=[0 0; 2 0 ; 1 2; 1 0.5]' - C=[ 1 2 4 ; 2 3 4 ; 3 1 4]' - CR=ones(Int,3) - F=[ 1 2 ; 2 3 ; 3 1]' - FR=[1,2,3] - g=simplexgrid(X,C,CR,F,FR) - bfc=g[BFaceCells] - tryfix(bfc)==[1 2 3;] + X = [0 0; 2 0; 1 2; 1 0.5]' + C = [1 2 4; 2 3 4; 3 1 4]' + CR = ones(Int, 3) + F = [1 2; 2 3; 3 1]' + FR = [1, 2, 3] + g = simplexgrid(X, C, CR, F, FR) + bfc = g[BFaceCells] + tryfix(bfc) == [1 2 3;] end - - end - - - -function testrw(grid,format;confidence=:full) +function testrw(grid, format; confidence = :full) #@warn format - ftmp=tempname()*"."*format - write(ftmp,grid) - grid1=simplexgrid(ftmp) - seemingly_equal_with_sort(grid1,grid;confidence=confidence) + ftmp = tempname() * "." * format + write(ftmp, grid) + grid1 = simplexgrid(ftmp) + seemingly_equal(grid1, grid; confidence = confidence) end @testset "Read/Write sg" begin - X=collect(0:0.05:1) - @test testrw(simplexgrid(X),"sg") - @test testrw(simplexgrid(X,X),"sg") - @test testrw(simplexgrid(X,X,X),"sg") + X = collect(0:0.05:1) + @test testrw(simplexgrid(X), "sg") + @test testrw(simplexgrid(X, X), "sg") + @test testrw(simplexgrid(X, X, X), "sg") end - -function testgrid(grid,testdata) - (num_nodes(grid),num_cells(grid), num_bfaces(grid))==testdata -end - -examples1d=joinpath(@__DIR__,"..","examples","examples1d.jl") +examples1d = joinpath(@__DIR__, "..", "examples", "examples1d.jl") include(examples1d) -examples2d=joinpath(@__DIR__,"..","examples","examples2d.jl") +examples2d = joinpath(@__DIR__, "..", "examples", "examples2d.jl") include(examples2d) -examples3d=joinpath(@__DIR__,"..","examples","examples3d.jl") +examples3d = joinpath(@__DIR__, "..", "examples", "examples3d.jl") include(examples3d) - - - @testset "1D" begin function rect1d() - X=collect(0:0.1:1) - g=simplexgrid(X) - rect!(g, [0.3],[0.6], region=2,bregions=[3,4]) + X = collect(0:0.1:1) + g = simplexgrid(X) + rect!(g, [0.3], [0.6]; region = 2, bregions = [3, 4]) end - - @test testgrid(interval_from_vector(),(21,20,2)) - @test testgrid(interval_localref(), (27,26,2)) - @test testgrid(interval_multiregion(),(21,20,3)) - @test testgrid(interval_subgrid(),(51,50,2)) - @test testgrid(rect1d(),(11,10,4)) + + @test testgrid(interval_from_vector(), (21, 20, 2)) + @test testgrid(interval_localref(), (27, 26, 2)) + @test testgrid(interval_multiregion(), (21, 20, 3)) + @test testgrid(interval_subgrid(), (51, 50, 2)) + @test testgrid(rect1d(), (11, 10, 4)) end @testset "2D" begin function rect2d() - X=collect(0:0.1:1) - g=simplexgrid(X,X) - rect!(g, [0.3, 0.3],[0.6,0.6], region=2,bregions=[3,4,5,6]) + X = collect(0:0.1:1) + g = simplexgrid(X, X) + rect!(g, [0.3, 0.3], [0.6, 0.6]; region = 2, bregions = [3, 4, 5, 6]) end - @test testgrid(rectangle(),(441,800,80)) - @test testgrid(rectangle_localref(),(729, 1352, 104)) - @test testgrid(rectangle_multiregion(),(441,800,100)) - @test testgrid(rectangle_subgrid(),(360, 600, 120)) - @test testgrid(rect2d(),(121,200,52)) - @test testgrid(rect2d_bregion_function(),(79,112,44)) - - g,sg,sf=sorted_subgrid() - @test testgrid(g,(187,306,66)) - @test testgrid(sg,(17,16,0)) - @test issorted(view(sg[Coordinates],1,:)) - + @test testgrid(rectangle(), (441, 800, 80)) + @test testgrid(rectangle_localref(), (729, 1352, 104)) + @test testgrid(rectangle_multiregion(), (441, 800, 100)) + @test testgrid(rectangle_subgrid(), (360, 600, 120)) + @test testgrid(rect2d(), (121, 200, 52)) + @test testgrid(rect2d_bregion_function(), (79, 112, 44)) + + g, sg, sf = sorted_subgrid() + @test testgrid(g, (187, 306, 66)) + @test testgrid(sg, (17, 16, 0)) + @test issorted(view(sg[Coordinates], 1, :)) end @testset "3D" begin - function subgen(;h=0.2) - X=collect(-1:h:2) - Y=collect(0:h:1) - Z=collect(0:h:2) - g=simplexgrid(X,Y,Z) + function subgen(; h = 0.2) + X = collect(-1:h:2) + Y = collect(0:h:1) + Z = collect(0:h:2) + g = simplexgrid(X, Y, Z) # Mark all boundaries as region 1 - bfacemask!(g,[-1,0,0], [2,1,2],1, allow_new=false) - + bfacemask!(g, [-1, 0, 0], [2, 1, 2], 1; allow_new = false) + # Default cell region is 1 # Mark cut-out regions as 2 - cellmask!(g, [-1,0,0], [0,1,1], 2) - cellmask!(g, [1,0,0], [2,1,1], 2) - + cellmask!(g, [-1, 0, 0], [0, 1, 1], 2) + cellmask!(g, [1, 0, 0], [2, 1, 1], 2) + # add new interface elements - bfacemask!(g, [-1,0,1], [2,1,1],3,allow_new=true) - bfacemask!(g, [0,0,0], [0,1,1],3,allow_new=true) - bfacemask!(g, [1,0,0], [1,1,1],3,allow_new=true) - + bfacemask!(g, [-1, 0, 1], [2, 1, 1], 3; allow_new = true) + bfacemask!(g, [0, 0, 0], [0, 1, 1], 3; allow_new = true) + bfacemask!(g, [1, 0, 0], [1, 1, 1], 3; allow_new = true) + # return subgrid of region 1 - subgrid(g,[1]) + subgrid(g, [1]) end function rect3d() - X=collect(0:0.1:1) - g=simplexgrid(X,X,X) - rect!(g, [0.3, 0.3,0.4],[0.6,0.6,0.7], region=2,bregion=8) + X = collect(0:0.1:1) + g = simplexgrid(X, X, X) + rect!(g, [0.3, 0.3, 0.4], [0.6, 0.6, 0.7]; region = 2, bregion = 8) g end - - - @test testgrid(quadrilateral(),(330,1200,440)) + @test testgrid(quadrilateral(), (330, 1200, 440)) @test mask_bedges() - X=collect(0:0.25:1) - gxy=simplexgrid(X,X) - gxyz=simplexgrid(gxy,X) - g=simplexgrid(X,X,X) - @test testgrid(gxyz,(125,384,192)) - @test g[Coordinates]≈gxyz[Coordinates] - @test testgrid(subgen(),(756,3000,950)) - @test testgrid(rect3d(),(1331,6000,1308)) - @test testgrid(cross3d(),(189,480,344)) + X = collect(0:0.25:1) + gxy = simplexgrid(X, X) + gxyz = simplexgrid(gxy, X) + g = simplexgrid(X, X, X) + @test testgrid(gxyz, (125, 384, 192)) + @test g[Coordinates] ≈ gxyz[Coordinates] + @test testgrid(subgen(), (756, 3000, 950)) + @test testgrid(rect3d(), (1331, 6000, 1308)) + @test testgrid(cross3d(), (189, 480, 344)) end - @testset "plotting examples" begin include("../docs/makeplots.jl") - picdir=mktempdir() - - @test makeplot("interval_from_vector",picdir; Plotter=CairoMakie) - @test makeplot("interval_localref",picdir; Plotter=CairoMakie) - @test makeplot("interval_multiregion",picdir; Plotter=CairoMakie) - @test makeplot("interval_subgrid",picdir; Plotter=CairoMakie) - @test makeplot("rectangle",picdir; Plotter=CairoMakie) - @test makeplot("rectangle_localref",picdir; Plotter=CairoMakie) - @test makeplot("rectangle_multiregion",picdir; Plotter=CairoMakie) - @test makeplot("rectangle_subgrid",picdir; Plotter=CairoMakie) - @test makeplot("quadrilateral",picdir; Plotter=CairoMakie) - @test makeplotx("sorted_subgrid",picdir; Plotter=CairoMakie) + picdir = mktempdir() + + @test makeplot("interval_from_vector", picdir; Plotter = CairoMakie) + @test makeplot("interval_localref", picdir; Plotter = CairoMakie) + @test makeplot("interval_multiregion", picdir; Plotter = CairoMakie) + @test makeplot("interval_subgrid", picdir; Plotter = CairoMakie) + @test makeplot("rectangle", picdir; Plotter = CairoMakie) + @test makeplot("rectangle_localref", picdir; Plotter = CairoMakie) + @test makeplot("rectangle_multiregion", picdir; Plotter = CairoMakie) + @test makeplot("rectangle_subgrid", picdir; Plotter = CairoMakie) + @test makeplot("quadrilateral", picdir; Plotter = CairoMakie) + @test makeplotx("sorted_subgrid", picdir; Plotter = CairoMakie) end +function tglue(; dim = 2, breg = 0) + X1 = linspace(0, 1, 5) + Y1 = linspace(0, 1, 5) + Z1 = linspace(0, 1, 5) -function tglue(;dim=2,breg=0) - X1=linspace(0,1,5) - Y1=linspace(0,1,5) - Z1=linspace(0,1,5) + X2 = linspace(1, 2, 5) + Y2 = linspace(0, 0.5, 3) + Z2 = linspace(0, 0.5, 3) - X2=linspace(1,2,5) - Y2=linspace(0,0.5,3) - Z2=linspace(0,0.5,3) - - if dim==1 - g1=simplexgrid(X1) - g2=simplexgrid(X2) + if dim == 1 + g1 = simplexgrid(X1) + g2 = simplexgrid(X2) end - - if dim==2 - g1=simplexgrid(X1,Y1) - g2=simplexgrid(X2,Y2) + + if dim == 2 + g1 = simplexgrid(X1, Y1) + g2 = simplexgrid(X2, Y2) end - - if dim==3 - g1=simplexgrid(X1,Y1,Z1) - g2=simplexgrid(X2,Y2,Z2) + + if dim == 3 + g1 = simplexgrid(X1, Y1, Z1) + g2 = simplexgrid(X2, Y2, Z2) end - - glue(g1,g2,interface=breg) + + glue(g1, g2; interface = breg) end @testset "Glue" begin - @test testgrid(tglue(;dim=1,breg=0),(9,8,2)) - @test testgrid(tglue(;dim=1,breg=1),(9,8,3)) - @test testgrid(tglue(;dim=2,breg=0),(37, 48, 24)) - @test testgrid(tglue(;dim=2,breg=1),(37, 48, 26)) - @test testgrid(tglue(;dim=3,breg=0),(161, 480, 256)) - @test testgrid(tglue(;dim=3,breg=1),(161, 480, 264)) + @test testgrid(tglue(; dim = 1, breg = 0), (9, 8, 2)) + @test testgrid(tglue(; dim = 1, breg = 1), (9, 8, 3)) + @test testgrid(tglue(; dim = 2, breg = 0), (37, 48, 24)) + @test testgrid(tglue(; dim = 2, breg = 1), (37, 48, 26)) + @test testgrid(tglue(; dim = 3, breg = 0), (161, 480, 256)) + @test testgrid(tglue(; dim = 3, breg = 1), (161, 480, 264)) end - - function voronoitest() - g=simplexgrid(0:0.1:1) - @test g[VoronoiFaceCenters]≈[0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95] - g=simplexgrid(0:0.5:1,0:0.5:1) - @test g[VoronoiFaceCenters]≈[0.25 0.5 0.25 0.25 0.0 0.75 1.0 0.75 0.75 0.5 0.25 0.25 0.0 1.0 0.75 0.75; - 0.0 0.25 0.25 0.5 0.25 0.0 0.25 0.25 0.5 0.75 0.75 1.0 0.75 0.75 0.75 1.0] + g = simplexgrid(0:0.1:1) + @test g[VoronoiFaceCenters] ≈ [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95] + g = simplexgrid(0:0.5:1, 0:0.5:1) + @test g[VoronoiFaceCenters] ≈ [0.25 0.5 0.25 0.25 0.0 0.75 1.0 0.75 0.75 0.5 0.25 0.25 0.0 1.0 0.75 0.75; + 0.0 0.25 0.25 0.5 0.25 0.0 0.25 0.25 0.5 0.75 0.75 1.0 0.75 0.75 0.75 1.0] end @testset "Voronoi" begin voronoitest() end - + @testset "Grid Stuff" begin include("test_gridstuff.jl") run_grid_tests() -end +end @testset "Interpolations" begin X_from = collect(0:0.1:1) X_to = collect(0:0.02:1) - grid_from=simplexgrid(X_from) - grid_to=simplexgrid(X_to) - u_from=map((x)->x,grid_from) - u_to_exact=map((x)->x,grid_to) - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - - u_from=Matrix(hcat(u_from,u_from)') - u_to_exact=Matrix(hcat(u_to_exact,u_to_exact)') - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - - grid_from=simplexgrid(X_from,X_from) - grid_to=simplexgrid(X_to,X_to) - u_from=map((x,y)->x+y,grid_from) - u_to_exact=map((x,y)->x+y,grid_to) - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - - u_from=Matrix(hcat(u_from,u_from)') - u_to_exact=Matrix(hcat(u_to_exact,u_to_exact)') - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - - grid_from=simplexgrid(X_from,X_from,X_from) - grid_to=simplexgrid(X_to,X_to,X_to) - u_from=map((x,y,z)->x+y+z,grid_from) - u_to_exact=map((x,y,z)->x+y+z,grid_to) - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - - u_from=Matrix(hcat(u_from,u_from)') - u_to_exact=Matrix(hcat(u_to_exact,u_to_exact)') - u_to=interpolate(grid_to,u_from,grid_from) - @test u_to≈u_to_exact - + grid_from = simplexgrid(X_from) + grid_to = simplexgrid(X_to) + u_from = map((x) -> x, grid_from) + u_to_exact = map((x) -> x, grid_to) + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact + + u_from = Matrix(hcat(u_from, u_from)') + u_to_exact = Matrix(hcat(u_to_exact, u_to_exact)') + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact + + grid_from = simplexgrid(X_from, X_from) + grid_to = simplexgrid(X_to, X_to) + u_from = map((x, y) -> x + y, grid_from) + u_to_exact = map((x, y) -> x + y, grid_to) + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact + + u_from = Matrix(hcat(u_from, u_from)') + u_to_exact = Matrix(hcat(u_to_exact, u_to_exact)') + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact + + grid_from = simplexgrid(X_from, X_from, X_from) + grid_to = simplexgrid(X_to, X_to, X_to) + u_from = map((x, y, z) -> x + y + z, grid_from) + u_to_exact = map((x, y, z) -> x + y + z, grid_to) + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact + + u_from = Matrix(hcat(u_from, u_from)') + u_to_exact = Matrix(hcat(u_to_exact, u_to_exact)') + u_to = interpolate(grid_to, u_from, grid_from) + @test u_to ≈ u_to_exact end @testset "writeVTK" begin X = collect(-1:1.0:1) Y = collect(-1:1.0:1) Z = collect(-2:1.0:2) - g = simplexgrid(X,Y,Z) + g = simplexgrid(X, Y, Z) nx = num_nodes(g) nc = num_cells(g) # ensure calculation of these data is free of roundoff errors - point_data = map((x,y,z) -> (x+y+z), g) + point_data = map((x, y, z) -> (x + y + z), g) field_data = [1.0, 2, 3, 4, 5, 6] - writeVTK("testfile_writevtk.vtu", g; - cellregions = g[CellRegions], - point_data = point_data, - field_data = field_data) + writeVTK("testfile_writevtk.vtu", g; + cellregions = g[CellRegions], + point_data = point_data, + field_data = field_data) sha_code = open("testfile_writevtk.vtu") do f sha256(f) @@ -543,21 +348,4 @@ end @test sha_code == "93a31139ccb3ae3017351d7cef0c2639c5def97c9744699543fe8bc58e1ebcea" end - - - - - - - - - - - - - - - - - - +include("gmsh.jl") diff --git a/test/testmesh.gmsh b/test/testmesh.gmsh new file mode 100644 index 00000000..d99a1580 --- /dev/null +++ b/test/testmesh.gmsh @@ -0,0 +1,366 @@ +$NOD +134 +1 0.000000 0.000000 0.000000 +2 1.000000 0.000000 0.000000 +3 1.000000 1.000000 0.000000 +4 0.000000 1.000000 0.000000 +5 0.100000 0.000000 0.000000 +6 0.200000 0.000000 0.000000 +7 0.300000 0.000000 0.000000 +8 0.400000 0.000000 0.000000 +9 0.500000 0.000000 0.000000 +10 0.600000 0.000000 0.000000 +11 0.700000 0.000000 0.000000 +12 0.800000 0.000000 0.000000 +13 0.900000 0.000000 0.000000 +14 1.000000 0.100000 0.000000 +15 1.000000 0.200000 0.000000 +16 1.000000 0.300000 0.000000 +17 1.000000 0.400000 0.000000 +18 1.000000 0.500000 0.000000 +19 1.000000 0.600000 0.000000 +20 1.000000 0.700000 0.000000 +21 1.000000 0.800000 0.000000 +22 1.000000 0.900000 0.000000 +23 0.900000 1.000000 0.000000 +24 0.800000 1.000000 0.000000 +25 0.700000 1.000000 0.000000 +26 0.600000 1.000000 0.000000 +27 0.500000 1.000000 0.000000 +28 0.400000 1.000000 0.000000 +29 0.300000 1.000000 0.000000 +30 0.200000 1.000000 0.000000 +31 0.100000 1.000000 0.000000 +32 -0.000000 0.900000 0.000000 +33 -0.000000 0.800000 0.000000 +34 -0.000000 0.700000 0.000000 +35 -0.000000 0.600000 0.000000 +36 -0.000000 0.500000 0.000000 +37 -0.000000 0.400000 0.000000 +38 -0.000000 0.300000 0.000000 +39 -0.000000 0.200000 0.000000 +40 -0.000000 0.100000 0.000000 +41 0.117837 0.129850 0.000000 +42 0.243574 0.110825 0.000000 +43 0.350998 0.102095 0.000000 +44 0.454129 0.095354 0.000000 +45 0.554997 0.088968 0.000000 +46 0.653715 0.082751 0.000000 +47 0.749032 0.076370 0.000000 +48 0.837509 0.071561 0.000000 +49 0.920410 0.079540 0.000000 +50 0.927991 0.162777 0.000000 +51 0.922271 0.252667 0.000000 +52 0.915416 0.351415 0.000000 +53 0.911849 0.454465 0.000000 +54 0.911404 0.557822 0.000000 +55 0.911852 0.660590 0.000000 +56 0.912830 0.763439 0.000000 +57 0.915174 0.866557 0.000000 +58 0.839928 0.921280 0.000000 +59 0.756266 0.909852 0.000000 +60 0.656133 0.907278 0.000000 +61 0.554888 0.910581 0.000000 +62 0.455518 0.912667 0.000000 +63 0.356263 0.914260 0.000000 +64 0.258227 0.918651 0.000000 +65 0.166885 0.923882 0.000000 +66 0.081781 0.917043 0.000000 +67 0.073865 0.831164 0.000000 +68 0.077755 0.738308 0.000000 +69 0.079674 0.638035 0.000000 +70 0.076437 0.538384 0.000000 +71 0.070888 0.445230 0.000000 +72 0.074090 0.355127 0.000000 +73 0.089809 0.257243 0.000000 +74 0.301556 0.209516 0.000000 +75 0.193642 0.225458 0.000000 +76 0.511342 0.184048 0.000000 +77 0.610968 0.171580 0.000000 +78 0.707077 0.159684 0.000000 +79 0.795936 0.145790 0.000000 +80 0.866837 0.132908 0.000000 +81 0.834733 0.297796 0.000000 +82 0.852608 0.204302 0.000000 +83 0.822635 0.404790 0.000000 +84 0.821605 0.513708 0.000000 +85 0.823378 0.618963 0.000000 +86 0.823505 0.724413 0.000000 +87 0.827344 0.830430 0.000000 +88 0.607818 0.817696 0.000000 +89 0.722337 0.798546 0.000000 +90 0.509348 0.825004 0.000000 +91 0.412228 0.826061 0.000000 +92 0.311663 0.827908 0.000000 +93 0.214098 0.842268 0.000000 +94 0.137960 0.858972 0.000000 +95 0.163205 0.678049 0.000000 +96 0.152484 0.780664 0.000000 +97 0.158669 0.572556 0.000000 +98 0.146712 0.477177 0.000000 +99 0.131259 0.401887 0.000000 +100 0.158117 0.326945 0.000000 +101 0.407948 0.197268 0.000000 +102 0.667919 0.249807 0.000000 +103 0.570914 0.264625 0.000000 +104 0.764654 0.232637 0.000000 +105 0.727084 0.468991 0.000000 +106 0.734354 0.577583 0.000000 +107 0.725765 0.344068 0.000000 +108 0.737723 0.679433 0.000000 +109 0.467640 0.742943 0.000000 +110 0.556540 0.740430 0.000000 +111 0.647185 0.721379 0.000000 +112 0.360015 0.306223 0.000000 +113 0.253801 0.314595 0.000000 +114 0.372738 0.734792 0.000000 +115 0.257975 0.729853 0.000000 +116 0.232889 0.499030 0.000000 +117 0.205650 0.405718 0.000000 +118 0.249153 0.604779 0.000000 +119 0.469574 0.288072 0.000000 +120 0.667098 0.636540 0.000000 +121 0.626572 0.550185 0.000000 +122 0.629842 0.425602 0.000000 +123 0.514462 0.677606 0.000000 +124 0.347399 0.631459 0.000000 +125 0.587824 0.654606 0.000000 +126 0.301943 0.411464 0.000000 +127 0.330448 0.520993 0.000000 +128 0.440803 0.654925 0.000000 +129 0.627473 0.330712 0.000000 +130 0.423478 0.418713 0.000000 +131 0.541644 0.367871 0.000000 +132 0.528525 0.484918 0.000000 +133 0.431077 0.551945 0.000000 +134 0.520194 0.595510 0.000000 +$ENDNOD +$ELM +226 +1 2 1 1 3 1 5 40 +2 2 1 1 3 5 6 41 +3 2 1 1 3 6 7 42 +4 2 1 1 3 6 42 41 +5 2 1 1 3 7 8 43 +6 2 1 1 3 7 43 42 +7 2 1 1 3 8 9 44 +8 2 1 1 3 8 44 43 +9 2 1 1 3 9 10 45 +10 2 1 1 3 9 45 44 +11 2 1 1 3 10 11 46 +12 2 1 1 3 10 46 45 +13 2 1 1 3 11 12 47 +14 2 1 1 3 11 47 46 +15 2 1 1 3 12 13 48 +16 2 1 1 3 12 48 47 +17 2 1 1 3 2 49 13 +18 2 1 1 3 13 49 48 +19 2 1 1 3 2 14 49 +20 2 1 1 3 14 15 50 +21 2 1 1 3 14 50 49 +22 2 1 1 3 15 16 51 +23 2 1 1 3 15 51 50 +24 2 1 1 3 16 17 52 +25 2 1 1 3 16 52 51 +26 2 1 1 3 17 18 53 +27 2 1 1 3 17 53 52 +28 2 1 1 3 18 19 54 +29 2 1 1 3 18 54 53 +30 2 1 1 3 19 20 55 +31 2 1 1 3 19 55 54 +32 2 1 1 3 20 21 56 +33 2 1 1 3 20 56 55 +34 2 1 1 3 21 22 57 +35 2 1 1 3 21 57 56 +36 2 1 1 3 3 23 22 +37 2 1 1 3 23 24 58 +38 2 1 1 3 24 25 59 +39 2 1 1 3 24 59 58 +40 2 1 1 3 25 26 60 +41 2 1 1 3 25 60 59 +42 2 1 1 3 26 27 61 +43 2 1 1 3 26 61 60 +44 2 1 1 3 27 28 62 +45 2 1 1 3 27 62 61 +46 2 1 1 3 28 29 63 +47 2 1 1 3 28 63 62 +48 2 1 1 3 29 30 64 +49 2 1 1 3 29 64 63 +50 2 1 1 3 30 31 65 +51 2 1 1 3 30 65 64 +52 2 1 1 3 4 66 31 +53 2 1 1 3 31 66 65 +54 2 1 1 3 4 32 66 +55 2 1 1 3 32 33 67 +56 2 1 1 3 32 67 66 +57 2 1 1 3 33 34 68 +58 2 1 1 3 33 68 67 +59 2 1 1 3 34 35 69 +60 2 1 1 3 34 69 68 +61 2 1 1 3 35 36 70 +62 2 1 1 3 35 70 69 +63 2 1 1 3 36 37 71 +64 2 1 1 3 36 71 70 +65 2 1 1 3 37 38 72 +66 2 1 1 3 37 72 71 +67 2 1 1 3 38 39 73 +68 2 1 1 3 38 73 72 +69 2 1 1 3 39 40 41 +70 2 1 1 3 5 41 40 +71 2 1 1 3 22 23 57 +72 2 1 1 3 23 58 57 +73 2 1 1 3 39 41 73 +74 2 1 1 3 42 43 74 +75 2 1 1 3 41 42 75 +76 2 1 1 3 42 74 75 +77 2 1 1 3 44 45 76 +78 2 1 1 3 45 46 77 +79 2 1 1 3 45 77 76 +80 2 1 1 3 46 47 78 +81 2 1 1 3 46 78 77 +82 2 1 1 3 47 48 79 +83 2 1 1 3 47 79 78 +84 2 1 1 3 48 49 80 +85 2 1 1 3 48 80 79 +86 2 1 1 3 49 50 80 +87 2 1 1 3 51 52 81 +88 2 1 1 3 50 51 82 +89 2 1 1 3 50 82 80 +90 2 1 1 3 51 81 82 +91 2 1 1 3 52 53 83 +92 2 1 1 3 52 83 81 +93 2 1 1 3 53 54 84 +94 2 1 1 3 53 84 83 +95 2 1 1 3 54 55 85 +96 2 1 1 3 54 85 84 +97 2 1 1 3 55 56 86 +98 2 1 1 3 55 86 85 +99 2 1 1 3 56 57 87 +100 2 1 1 3 56 87 86 +101 2 1 1 3 57 58 87 +102 2 1 1 3 60 61 88 +103 2 1 1 3 59 60 89 +104 2 1 1 3 60 88 89 +105 2 1 1 3 58 59 87 +106 2 1 1 3 61 62 90 +107 2 1 1 3 61 90 88 +108 2 1 1 3 62 63 91 +109 2 1 1 3 62 91 90 +110 2 1 1 3 63 64 92 +111 2 1 1 3 63 92 91 +112 2 1 1 3 64 65 93 +113 2 1 1 3 64 93 92 +114 2 1 1 3 65 66 94 +115 2 1 1 3 65 94 93 +116 2 1 1 3 66 67 94 +117 2 1 1 3 68 69 95 +118 2 1 1 3 67 68 96 +119 2 1 1 3 67 96 94 +120 2 1 1 3 68 95 96 +121 2 1 1 3 69 70 97 +122 2 1 1 3 69 97 95 +123 2 1 1 3 70 71 98 +124 2 1 1 3 70 98 97 +125 2 1 1 3 71 72 99 +126 2 1 1 3 71 99 98 +127 2 1 1 3 72 73 100 +128 2 1 1 3 72 100 99 +129 2 1 1 3 43 44 101 +130 2 1 1 3 44 76 101 +131 2 1 1 3 59 89 87 +132 2 1 1 3 43 101 74 +133 2 1 1 3 73 75 100 +134 2 1 1 3 41 75 73 +135 2 1 1 3 77 78 102 +136 2 1 1 3 76 77 103 +137 2 1 1 3 77 102 103 +138 2 1 1 3 79 80 82 +139 2 1 1 3 79 82 104 +140 2 1 1 3 83 84 105 +141 2 1 1 3 84 85 106 +142 2 1 1 3 84 106 105 +143 2 1 1 3 81 83 107 +144 2 1 1 3 83 105 107 +145 2 1 1 3 85 86 108 +146 2 1 1 3 85 108 106 +147 2 1 1 3 86 87 89 +148 2 1 1 3 90 91 109 +149 2 1 1 3 88 90 110 +150 2 1 1 3 90 109 110 +151 2 1 1 3 88 111 89 +152 2 1 1 3 88 110 111 +153 2 1 1 3 78 79 104 +154 2 1 1 3 74 101 112 +155 2 1 1 3 74 113 75 +156 2 1 1 3 75 113 100 +157 2 1 1 3 74 112 113 +158 2 1 1 3 91 92 114 +159 2 1 1 3 91 114 109 +160 2 1 1 3 92 93 115 +161 2 1 1 3 92 115 114 +162 2 1 1 3 93 94 96 +163 2 1 1 3 97 98 116 +164 2 1 1 3 98 99 117 +165 2 1 1 3 98 117 116 +166 2 1 1 3 95 97 118 +167 2 1 1 3 97 116 118 +168 2 1 1 3 100 113 117 +169 2 1 1 3 86 89 108 +170 2 1 1 3 76 119 101 +171 2 1 1 3 101 119 112 +172 2 1 1 3 76 103 119 +173 2 1 1 3 78 104 102 +174 2 1 1 3 81 104 82 +175 2 1 1 3 81 107 104 +176 2 1 1 3 89 111 108 +177 2 1 1 3 93 96 115 +178 2 1 1 3 99 100 117 +179 2 1 1 3 95 115 96 +180 2 1 1 3 95 118 115 +181 2 1 1 3 102 104 107 +182 2 1 1 3 108 111 120 +183 2 1 1 3 106 108 120 +184 2 1 1 3 105 106 121 +185 2 1 1 3 106 120 121 +186 2 1 1 3 105 122 107 +187 2 1 1 3 105 121 122 +188 2 1 1 3 109 123 110 +189 2 1 1 3 114 115 124 +190 2 1 1 3 115 118 124 +191 2 1 1 3 110 125 111 +192 2 1 1 3 111 125 120 +193 2 1 1 3 110 123 125 +194 2 1 1 3 112 126 113 +195 2 1 1 3 113 126 117 +196 2 1 1 3 116 117 126 +197 2 1 1 3 116 127 118 +198 2 1 1 3 118 127 124 +199 2 1 1 3 116 126 127 +200 2 1 1 3 109 114 128 +201 2 1 1 3 109 128 123 +202 2 1 1 3 114 124 128 +203 2 1 1 3 102 129 103 +204 2 1 1 3 112 130 126 +205 2 1 1 3 112 119 130 +206 2 1 1 3 107 122 129 +207 2 1 1 3 102 107 129 +208 2 1 1 3 103 131 119 +209 2 1 1 3 119 131 130 +210 2 1 1 3 103 129 131 +211 2 1 1 3 122 132 131 +212 2 1 1 3 122 131 129 +213 2 1 1 3 120 125 121 +214 2 1 1 3 121 132 122 +215 2 1 1 3 124 127 133 +216 2 1 1 3 124 133 128 +217 2 1 1 3 121 125 134 +218 2 1 1 3 121 134 132 +219 2 1 1 3 123 134 125 +220 2 1 1 3 126 130 127 +221 2 1 1 3 123 128 134 +222 2 1 1 3 130 131 132 +223 2 1 1 3 127 130 133 +224 2 1 1 3 128 133 134 +225 2 1 1 3 130 132 133 +226 2 1 1 3 132 134 133 +$ENDELM$