diff --git a/Project.toml b/Project.toml index 8a7e8c94..c18a3330 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableGrids" uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8" authors = ["Juergen Fuhrmann ", "Christian Merdon ", "Johannes Taraz ", "Patrick Jaap "] -version = "1.14.0" +version = "1.15.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl index be3d404a..8bf3ecf5 100644 --- a/src/ExtendableGrids.jl +++ b/src/ExtendableGrids.jl @@ -12,6 +12,7 @@ using Printf using Random using Dates using LinearAlgebra +using WriteVTK import Graphs @@ -184,6 +185,7 @@ export TokenStream, gettoken, expecttoken, trytoken include("io.jl") export writeVTK +export BRegionDomCode include("seal.jl") diff --git a/src/extendablegrid.jl b/src/extendablegrid.jl index df0327ce..0200ddfa 100644 --- a/src/extendablegrid.jl +++ b/src/extendablegrid.jl @@ -614,55 +614,83 @@ Keyword args: - :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 - return true +function seemingly_equal(grid1::ExtendableGrid, grid2::ExtendableGrid; sort = false, confidence = :full, skipkeys = [], trimgrids = true, keep = []) + if trimgrids + grid1 = trim(grid1; keep) + grid2 = trim(grid2; keep) 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 + if !sort + for key in keys(grid1) + if key in skipkeys + continue + end + if !haskey(grid2, key) + @warn "missing key $(key) in grid2" + return false + end + if !seemingly_equal(grid1[key], grid2[key]) + @warn "grid1[$(key)] and grid2[$(key)] differ" + return false + end 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 + return true + else + for key in keys(grid1) + if key in skipkeys + continue + end + if !haskey(grid2, key) + @warn "missing key $(key) in grid2" + return false + end + if !isa(grid1[key], AbstractArray) + if !seemingly_equal(grid1[key], grid2[key]) + @warn "grid1[$(key)] and grid2[$(key)] differ" + 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 + 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] + if !seemingly_equal(sa1, sa2) + @warn "grid1[$(key)] and grid2[$(key)] differ" + return false + end + else + if !seemingly_equal(grid1[key], grid2[key]) + @warn "grid1[$(key)] and grid2[$(key)] differ" + return false + end + end else - !seemingly_equal(grid1[key], grid2[key]) && return false + if eltype(grid1[key]) <: Number && key != Coordinates + sa1 = multidimsort(sort(grid1[key]; dims = 1)) + sa2 = multidimsort(sort(grid2[key]; dims = 1)) + if !seemingly_equal(sa1, sa2) + @warn "grid1[$(key)] and grid2[$(key)] differ" + return false + end + else + if !seemingly_equal(grid1[key], grid2[key]) + @warn "grid1[$(key)] and grid2[$(key)] differ" + return false + end + end end end end + return true 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)) diff --git a/src/io.jl b/src/io.jl index 553aa026..68f20f03 100644 --- a/src/io.jl +++ b/src/io.jl @@ -1,5 +1,6 @@ -using WriteVTK - +#################### +## VTK +################### # conversion from AbstractElementGeometry to WriteVTK.VTKCellTypes WriteVTK.VTKCellType(::Type{<:AbstractElementGeometry1D}) = VTKCellTypes.VTK_LINE WriteVTK.VTKCellType(::Type{<:Triangle2D}) = VTKCellTypes.VTK_TRIANGLE @@ -10,15 +11,15 @@ WriteVTK.VTKCellType(::Type{<:Hexahedron3D}) = VTKCellTypes.VTK_HEXAHEDRON """ $(TYPEDSIGNATURES) -exports grid and optional provided data as a vtk file + exports grid and optional provided data as a vtk file -- `filename`: filename of the exported file -- `grid`: grid + - `filename`: filename of the exported file + - `grid`: grid -Each '(key, value)' pair adds another data entry to the vtk file via WriteVTK functionality. + Each '(key, value)' pair adds another data entry to the vtk file via WriteVTK functionality. -For the arguments 'append' and 'compress', see documentation of vtk_grid of WriteVTK. -""" + For the arguments 'append' and 'compress', see documentation of vtk_grid of WriteVTK. + """ function writeVTK(filename::String, grid::ExtendableGrid{Tc, Ti}; append = false, compress = false, kwargs...) where {Tc, Ti} ncells = num_cells(grid) # get number of cells in grid coords = grid[Coordinates] # get coordinates @@ -40,143 +41,68 @@ function writeVTK(filename::String, grid::ExtendableGrid{Tc, Ti}; append = false end end +######################################################## +## Read/Write grids +####################################################### + """ $(TYPEDSIGNATURES) - -Write grid to file. -Supported formats: -- "*.sg": pdelib sg format. See [`simplexgrid(::String;kwargs...)`](@ref) - -""" -function Base.write(fname::String, g::ExtendableGrid; format = "", kwargs...) - (fbase, fext) = splitext(fname) + + Read grid from file. + Supported formats: + - "*.sg": pdelib sg files. Format versions: + - format version 2.0: long version with some unnecessary data + - format version 2.1: shortened version only with cells, cellnodes, cellregions, bfacenodes, bfaceregions + - format version 2.2: like 2.1, but additional info on cell and node partitioning. Edge partitioning + is not stored in the file and may be re-established by [`induce_edge_partitioning!`](@ref). + - "*.geo": gmsh geometry description (requires `using Gmsh`) + - "*.msh": gmsh mesh (requires `using Gmsh`) + - "*.dom": WIAS-TeSCA dom format + - The resulting grid contains `grid[BRegionDomCode]` describing the boundary condition code for + each boundary region + """ +function simplexgrid(file::String; format = "", kwargs...) + (fbase, fext) = splitext(file) if format == "" format = fext[2:end] end return try - writegrid(fname, g, Val{Symbol(format)}; kwargs...) - catch e - throw(ErrorException("Writing $(fext) files not supported")) - end -end - -function writegrid(filename::String, g::ExtendableGrid, ::Type{Val{:msh}}; kwargs...) - return try - simplexgrid_to_gmsh(g; filename) + simplexgrid(file, Val{Symbol(format)}; kwargs...) catch e - throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to write msh files.")) + throw(ErrorException("Reading $(fext) files not supported")) end end -function writegrid(fname::String, g::ExtendableGrid, ::Type{Val{:sg}}; version = v"2.2", kwargs...) - 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 something non-allocating - open(fname, "w") do file - write(file, @sprintf("SimplexGrid")) - write(file, @sprintf(" ")) - write(file, @sprintf("%s\n", "$version"[1:(end - 2)])) - write(file, @sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n")) - write(file, @sprintf("#%s\n", Dates.format(Dates.now(), "yyyy-mm-ddTHH-mm-SS"))) - - write(file, @sprintf("DIMENSION %d\n", dim_g)) - write(file, @sprintf("NODES %d %d\n", nn, dim_s)) - - for inode in 1:nn - for idim in 1:dim_s - write(file, @sprintf("%.20e ", coord[idim, inode])) - write(file, @sprintf("\n")) - end - end - - write(file, @sprintf("CELLS %d\n", nc)) - for icell in 1:nc - for inode in 1:(dim_g + 1) - write(file, @sprintf("%d ", cellnodes[inode, icell])) - end - write(file, @sprintf("%d\n", cellregions[icell])) - end - - write(file, @sprintf("FACES %d\n", nbf)) - for ibface in 1:nbf - for inode in 1:dim_g - write(file, @sprintf("%d ", bfacenodes[inode, ibface])) - end - write(file, @sprintf("%d\n", bfaceregions[ibface])) - end - - if version > v"2.1" - function writeitems(key, label) - data = g[key] - write(file, "$(label) $(length(data))\n") - for d in data - write(file, "$d\n") - end - return - end - writeitems(PColorPartitions, "PCOLORPARTITIONS") - writeitems(PartitionCells, "PARTITIONCELLS") - writeitems(PartitionBFaces, "PARTITIONBFACES") - writeitems(PartitionNodes, "PARTITIONNODES") - end - write(file, @sprintf("END\n")) - flush(file) - flush(file) - end - return nothing -end -###################################################### """ $(TYPEDSIGNATURES) - -Read grid from file. + +Write grid to file. Supported formats: -- "*.sg": pdelib sg files. Format versions: - - `format=v"2.0"`: long version with some unnecessary data - - `format=v"2.1"`: shortened version only with cells, cellnodes, cellregions, bfacenodes, bfaceregions - - `format=v"2.2"`: like 2.1, but additional info on cell and node partitioning. Edge partitioning is not - stored in the file and may be re-established by [`induce_edge_partitioning!`](@ref). -- "*.geo": gmsh geometry description (requires `using Gmsh`) -- "*.msh": gmsh mesh (requires `using Gmsh`) +- "*.sg": pdelib sg format. Keywords: + - `version`: format version +- "*.msh": gmsh grid format +- "*.dom": WIAS-TeSCA dom format. Keywords: + - domcodes: Vector of boundary codes per boundary region. + +See [`simplexgrid(::String;kwargs...)`](@ref) """ -function simplexgrid(file::String; format = "", kwargs...) - (fbase, fext) = splitext(file) +function Base.write(fname::String, g::ExtendableGrid; format = "", kwargs...) + (fbase, fext) = splitext(fname) if format == "" format = fext[2:end] end return try - simplexgrid(file, Val{Symbol(format)}; kwargs...) - catch e - throw(ErrorException("Reading $(fext) files not supported")) - end -end - -function simplexgrid(file::String, ::Type{Val{:msh}}; kwargs...) - return try - simplexgrid_from_gmsh(file) + writegrid(fname, g, Val{Symbol(format)}; kwargs...) catch e - throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to read msh files.")) + throw(ErrorException("Writing $(fext) files not supported")) end end -function simplexgrid(file::String, ::Type{Val{:geo}}; kwargs...) - return try - simplexgrid_from_gmsh(file) - catch e - throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to read geo files.")) - end -end +######################################################## +## PDELIB sg format +####################################################### function simplexgrid(file::String, ::Type{Val{:sg}}; kwargs...) Ti = Cint tks = TokenStream(file) @@ -272,6 +198,219 @@ function simplexgrid(file::String, ::Type{Val{:sg}}; kwargs...) return g end +function writegrid(fname::String, g::ExtendableGrid, ::Type{Val{:sg}}; version = v"2.2", kwargs...) + 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 something non-allocating + open(fname, "w") do file + write(file, @sprintf("SimplexGrid")) + write(file, @sprintf(" ")) + write(file, @sprintf("%s\n", "$version"[1:(end - 2)])) + write(file, @sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n")) + write(file, @sprintf("#%s\n", Dates.format(Dates.now(), "yyyy-mm-ddTHH-mm-SS"))) + + write(file, @sprintf("DIMENSION %d\n", dim_g)) + write(file, @sprintf("NODES %d %d\n", nn, dim_s)) + + for inode in 1:nn + for idim in 1:dim_s + write(file, @sprintf("%.20e ", coord[idim, inode])) + write(file, @sprintf("\n")) + end + end + + write(file, @sprintf("CELLS %d\n", nc)) + for icell in 1:nc + for inode in 1:(dim_g + 1) + write(file, @sprintf("%d ", cellnodes[inode, icell])) + end + write(file, @sprintf("%d\n", cellregions[icell])) + end + + write(file, @sprintf("FACES %d\n", nbf)) + for ibface in 1:nbf + for inode in 1:dim_g + write(file, @sprintf("%d ", bfacenodes[inode, ibface])) + end + write(file, @sprintf("%d\n", bfaceregions[ibface])) + end + + if version > v"2.1" + function writeitems(key, label) + data = g[key] + write(file, "$(label) $(length(data))\n") + for d in data + write(file, "$d\n") + end + return + end + writeitems(PColorPartitions, "PCOLORPARTITIONS") + writeitems(PartitionCells, "PARTITIONCELLS") + writeitems(PartitionBFaces, "PARTITIONBFACES") + writeitems(PartitionNodes, "PARTITIONNODES") + end + write(file, @sprintf("END\n")) + flush(file) + flush(file) + end + return nothing +end + +######################################################## +## WIAS-TeSCA dom format +####################################################### +struct BRegionDomCode <: AbstractGridComponent end + +function simplexgrid( + file::String, ::Type{Val{:dom}}; + kwargs... + ) + Ti = Cint + tks = TokenStream(file) + nn = parse(Ti, gettoken(tks)) + nc = parse(Ti, gettoken(tks)) + + coord = Array{Float64, 2}(undef, 2, nn) + cellnodes = Array{Ti, 2}(undef, 3, nc) + cellregions = Vector{Ti}(undef, nc) + bfacenodes = ElasticArray{Ti, 2}(undef, 2, 0) + bfaceregions = Vector{Ti}(undef, 0) + + for inode in 1:nn + coord[1, inode] = parse(Float64, gettoken(tks)) + coord[2, inode] = parse(Float64, gettoken(tks)) + end + + function add_bface!(bfacenodes, bfaceregions, cn1, cn2, bfr) + n1 = min(cn1, cn2) + n2 = max(cn1, cn2) + found = false + for ibface in 1:size(bfacenodes, 2) + if bfacenodes[1] == n1 && bfacenodes[2] == n2 + found = true + end + end + if !found && bfr != 0 + append!(bfacenodes, (n1, n2)) + push!(bfaceregions, bfr) + end + return + end + + for icell in 1:nc + for inode in 1:3 + cellnodes[inode, icell] = parse(Ti, gettoken(tks)) + end + add_bface!(bfacenodes, bfaceregions, cellnodes[2, icell], cellnodes[3, icell], parse(Ti, gettoken(tks))) + add_bface!(bfacenodes, bfaceregions, cellnodes[1, icell], cellnodes[3, icell], parse(Ti, gettoken(tks))) + add_bface!(bfacenodes, bfaceregions, cellnodes[1, icell], cellnodes[2, icell], parse(Ti, gettoken(tks))) + cellregions[icell] = parse(Ti, gettoken(tks)) + end + domcodes = sort(unique(bfaceregions)) + for i in 1:length(bfaceregions) + bfaceregions[i] = findfirst(r -> r == bfaceregions[i], domcodes) + end + println("Boundary code replacement (see grid[BRegionDomCode]):") + for i in 1:length(domcodes) + @printf("%4d => %d\n", domcodes[i], i) + end + g = simplexgrid(coord, cellnodes, cellregions, bfacenodes, bfaceregions) + g[BRegionDomCode] = domcodes + return g +end + +function writegrid( + fname::String, g::ExtendableGrid, ::Type{Val{:dom}}; + domcodes = collect(1:num_bfaceregions(g)), + kwargs... + ) + println("Boundary code replacement:") + for i in 1:length(domcodes) + @printf("%4d => %d\n", i, domcodes[i]) + end + dim_g = dim_grid(g) + dim_s = dim_space(g) + if dim_g != 2 + error("dom files are defined for 2D only") + end + nn = num_nodes(g) + nc = num_cells(g) + nbf = num_bfaces(g) + nbr = num_bfaceregions(g) + + coord = g[Coordinates] + cellnodes = g[CellNodes] + bfacenodes = g[BFaceNodes] + cellregions = g[CellRegions] + bfaceregions = g[BFaceRegions] + cellfaces = g[CellFaces] + bfacefaces = g[BFaceFaces] + nf = maximum(cellfaces) + + + facemarkers = zeros(Int, nf) + for ibf in 1:nbf + facemarkers[bfacefaces[ibf]] = domcodes[bfaceregions[ibf]] + end + + open(fname, "w") do file + write(file, @sprintf("%d %d\n", nn, nc)) + for inode in 1:nn + for idim in 1:dim_s + write(file, @sprintf("%.20e ", coord[idim, inode])) + end + write(file, @sprintf("\n")) + end + + for icell in 1:nc + for inode in 1:(dim_g + 1) + write(file, @sprintf("%d ", cellnodes[inode, icell])) + end + write(file, @sprintf("%d ", facemarkers[cellfaces[2, icell]])) + write(file, @sprintf("%d ", facemarkers[cellfaces[3, icell]])) + write(file, @sprintf("%d ", facemarkers[cellfaces[1, icell]])) + write(file, @sprintf("%d\n", cellregions[icell])) + end + end + return nothing +end + +######################################################## +## Wrapper to gmsh extension +####################################################### +function simplexgrid(file::String, ::Type{Val{:msh}}; kwargs...) + return try + simplexgrid_from_gmsh(file) + catch e + throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to read msh files.")) + end +end + +function simplexgrid(file::String, ::Type{Val{:geo}}; kwargs...) + return try + simplexgrid_from_gmsh(file) + catch e + throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to read geo files.")) + end +end +function writegrid(filename::String, g::ExtendableGrid, ::Type{Val{:msh}}; kwargs...) + return try + simplexgrid_to_gmsh(g; filename) + catch e + throw(ErrorException("Missing Gmsh extension. Add Gmsh.jl to your environment and import it to write msh files.")) + end +end + + function simplexgrid_from_gmsh end function simplexgrid_to_gmsh end diff --git a/test/runtests.jl b/test/runtests.jl index 41c02f66..c4c79fd6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -222,12 +222,12 @@ end @test map(vfxyz, gxyz) == map(vfv, gxyz) end -function testrw(grid, format; confidence = :full, kwargs...) +function testrw(grid, format; compare_kwargs = (confidence = :full,), kwargs...) #@warn format ftmp = tempname() * "." * format write(ftmp, grid; kwargs...) grid1 = simplexgrid(ftmp) - return seemingly_equal(grid1, grid; confidence = confidence) + return seemingly_equal(grid1, grid; compare_kwargs...) end @testset "Read/Write sg" begin @@ -239,6 +239,12 @@ end end end +@testset "Read/Write dom" begin + X = collect(0:0.05:1) + @test testrw(simplexgrid(X, X), "dom"; compare_kwargs = (sort = true, skipkeys = [XCoordinates, YCoordinates, BRegionDomCode])) +end + + @testset "rectnd" begin function rect1d() X = collect(0:0.1:1)