JOcTree Reorganization

src/IO/exportUBC.jl

@@ -0,0 +1,77 @@
+
+export exportUBCOcTreeMesh, exportUBCOcTreeModel
+
+
+"""
+    exportUBCOcTreeMesh(fname, mesh)
+
+    Writes an OcTree mesh to disk in UBC format.
+
+    Input:
+
+        name::AbstractString - File to write
+        mesh::OcTreeMesh     - The mesh to export
+
+"""
+function exportUBCOcTreeMesh(fname::AbstractString, mesh::OcTreeMesh)
+
+    m1,m2,m3     = mesh.n
+    i1,i2,i3,bsz = find3(mesh.S)
+    h1,h2,h3     = mesh.h
+    x1,x2,x3     = mesh.x0
+
+    # Roman's code starts the OcTree at the top corner. Change from bottom
+    # corner.
+    i3 = m3 + 2 .- i3 - bsz
+    x3 = x3 + m3 * h3
+    S = sub2ind( (m1,m2,m3), i1,i2,i3 )
+    p = sortpermFast(S)[1]
+    n = length(bsz)
+
+    # Write OcTree mesh
+    f = open(fname, "w")
+    println(f, m1, " ", m2, " ", m3, " ! # of cells in underlying mesh")
+    println(f, x1, " ", x2, " ", x3, " ! top corner")
+    println(f, h1, " ", h2, " ", h3, " ! cell size")
+    println(f, n, " ! size of octree mesh")
+    for i = 1:n
+        idx = p[i]
+        @printf(f,"%i %i %i %i\n", i1[idx], i2[idx], i3[idx], bsz[idx])
+    end
+    close(f)
+    return
+end
+
+"""
+    exportUBCOcTreeModel(fname, mesh, model)
+
+    Writes an OcTree model to disk in UBC format.
+
+    Input:
+
+        name::AbstractString - File to write
+        mesh::OcTreeMesh     - The mesh corresponding to the model
+        model::Union{Array{Float64,1}, Array{Int64,1}} - The model 
+
+"""
+function exportUBCOcTreeModel(name::AbstractString, mesh::OcTreeMesh, model::Union{Array{Float64,1}, Array{Int64,1}})
+
+    m1,m2,m3     = mesh.n
+    i1,i2,i3,bsz = find3(mesh.S)
+
+    # Roman's code starts the OcTree at the top corner. Change from bottom
+    # corner.
+    i3 = m3 + 2 .- i3 - bsz
+    n = nnz(mesh.S)
+    S = sub2ind( (m1,m2,m3), i1,i2,i3 )
+    p = sortpermFast(S)[1]
+    modelPerm = model[p]
+
+    # Write model vector
+    f = open(name, "w")
+    for i = 1:n
+        println(f, modelPerm[i])
+    end
+    close(f)
+    return
+end

src/IO/importUBC.jl

@@ -0,0 +1,140 @@
+export importUBCOcTreeMesh, importUBCOcTreeModel
+
+"""
+    mesh = importUBCOcTreeMesh(meshfile)
+
+    Reads an OcTree mesh in UBC format from disk.
+
+    Input:
+
+        meshfile::AbstractString - File to read
+
+    Output:
+
+        mesh::OcTreeMesh - The mesh
+
+"""
+function importUBCOcTreeMesh(meshfile::AbstractString)
+
+    # open file (throws error if file doesn't exist)
+    f    = open(meshfile,"r")
+
+    # number of cells of underlying tensor mesh along dimension 1, 2, 3
+    line = split(readline(f))
+    m1 = parse(Int64,line[1])
+    m2 = parse(Int64,line[2])
+    m3 = parse(Int64,line[3])
+
+    # top corner coordinates
+    line = split(readline(f))
+    x1 = parse(Float64,line[1])
+    x2 = parse(Float64,line[2])
+    x3 = parse(Float64,line[3])
+
+    # cell size
+    line = split(readline(f))
+    h1 = parse(Float64,line[1])
+    h2 = parse(Float64,line[2])
+    h3 = parse(Float64,line[3])
+
+    # number of OcTree cells
+    line = split(readline(f))
+    n = parse(Int64,line[1])
+
+    # read rest of file at ones
+    lines = readlines(f)
+
+    # close file
+    close(f)
+
+    # check correct number of lines read
+    if n != length(lines)
+        error("Invalid number of (i,j,k,bsz) lines in file $meshfile.")
+    end
+
+    # allocate space
+    i1  = zeros(Int64, n)
+    i2  = zeros(Int64, n)
+    i3  = zeros(Int64, n)
+    bsz = zeros(Int64, n)
+
+    # convert string array to numbers
+    for i = 1:n
+        line   = split(lines[i])
+
+        i1[i]  = parse(Int64,line[1])
+        i2[i]  = parse(Int64,line[2])
+        i3[i]  = parse(Int64,line[3])
+        bsz[i] = parse(Int64,line[4])
+    end
+
+    # Roman's code starts the OcTree at the top corner. Change to bottom
+    # corner.
+    i3 = m3 + 2 .- i3 - bsz
+    x3 = x3 - m3 * h3
+
+    #S   = sortrows([i3 i2 i1 bsz])
+    S = sub2ind( (m1,m2,m3), i1,i2,i3 )
+    p = sortpermFast(S)[1]
+
+    i1  = i1[p] #S[:,3]
+    i2  = i2[p] #S[:,2]
+    i3  = i3[p] #S[:,1]
+    bsz = bsz[p] #S[:,4]
+
+    # create mesh object
+    S = sparse3(i1,i2,i3,bsz,[m1,m2,m3])
+    mesh = getOcTreeMeshFV(S,[h1,h2,h3];x0=[x1,x2,x3])
+    return mesh
+end
+
+"""
+    model = importUBCOcTreeMesh(modelfile, mesh)
+
+    Reads an OcTree model in UBC format from disk.
+
+    Input:
+
+        modelfile::AbstractString - File to read
+        mesh::OcTreeMeshFV        - The corresponding mesh
+
+    Output:
+
+        model::Array{Float64,1} - The model
+
+"""
+function importUBCOcTreeModel(modelfile::AbstractString, mesh::OcTreeMesh)
+
+    # open file (throws error if file doesn't exist)
+    f    = open(modelfile,"r")
+
+    # read everything
+    s = readlines(f)
+
+    # close
+    close(f)
+
+    # convert to numbers
+    model = Array{Float64}(length(s))
+    for i = 1:length(s)
+        model[i] = parse(Float64,s[i])
+    end
+
+    # check if we have the correct number of cell values
+    if length(model) != mesh.nc
+        error("Incorrect number of cell values")
+    end
+
+    # Roman's code starts the OcTree at the top corner. Here, we start with the bottom corner. Therefore, we need to permute the cells values.
+    m1,m2,m3 = mesh.n
+    i1,i2,i3,bsz = find3(mesh.S)
+    i3 = m3 + 2 .- i3 - bsz
+
+    n = nnz(mesh.S)
+
+    S = sub2ind( (m1,m2,m3), i1,i2,i3 )
+    p = sortpermFast(S)[1]
+    ipermute!(model,p)
+
+    return model
+end

src/JOcTree.jl

@@ -1,75 +1,70 @@
 module JOcTree
 
-# using jInv.Mesh.AbstractMesh
-# using jInv.Mesh.ndgrid
-importall jInv.Mesh
-using jInv.Utils
-export OcTreeMesh
-abstract OcTreeMesh <: AbstractMesh
+    importall jInv.Mesh
+    using jInv.Utils
+    export OcTreeMesh
+    abstract OcTreeMesh <: AbstractMesh
 
-include("sparse3.jl")
-include("OcTreeMeshFV.jl")
-include("OcTreeMeshFEM.jl")
+    include("sparse3.jl")
+    include("OcTreeMeshFV.jl")
+    include("display.jl")
 
-include("findNonRegularBlocks.jl")
-include("findBlocks.jl")
-include("getCellNumbering.jl")
-include("getCellCenteredGrid.jl")
-include("getCurlMatrixRec.jl")
-include("getDivergenceMatrixRec.jl")
-include("getEdgeToCellCenteredMatrix.jl")
-include("getEdgeMassMatrix.jl")
-include("getEdgeMassMatrixAnisotropic.jl")
-include("getEdgeMassMatrixAnisotropicNoWeight.jl")
-include("getEdgeNumbering.jl")
-include("getEdgeGrids.jl")
-include("getEdgeSize.jl")
-include("getEdgeInterpolationMatrix.jl")
-include("getFaceToCellCenteredMatrix.jl")
-include("getFaceMassMatrix.jl")
-include("getFaceGrids.jl")
-include("getFaceSize.jl")
-include("getFaceNumbering.jl")
-include("getFaceInterpolationMatrix.jl")
-include("getInterpolationMatrix.jl")
-include("getNodalConstraints.jl")
-include("getFaceConstraints.jl")
-include("getNodalMassMatrix.jl")
-include("getNodalInterpolationMatrix.jl")
-include("getNodalNumbering.jl")
-include("getNodalGradientRec.jl")
-include("getNodalGrid.jl")
-include("getNodesToCellCenteredMatrix.jl")
-include("getNumberOfNeighbors.jl")
-include("regularizeOcTree.jl")
-include("refineOcTree.jl")
-include("getVolume.jl")
-include("getLength.jl")
-include("uniteOcTrees.jl")
+    include("getNumbering.jl")
+    include("getGrids.jl")
+    include("getEdgeSize.jl")
+    include("getFaceSize.jl")
 
-include("createOcTreeFromBox.jl")
-include("getEdgeIntegralOfPolygonalChain.jl")
-include("createOcTreeFromImage.jl")
-include("initCoarseOcTree.jl")
+    include("getCurlMatrixRec.jl")
+    include("getDivergenceMatrixRec.jl")
+    include("getNodalGradientRec.jl")
 
-include("display.jl")
+    include("getEdgeMassMatrix.jl")
+    include("getEdgeMassMatrixAnisotropic.jl")
+    include("getEdgeMassMatrixAnisotropicNoWeight.jl")
+    include("getFaceMassMatrix.jl")
+    include("getNodalMassMatrix.jl")
 
-include("getEdgeConstraints.jl")
-include("createOcTreeMesh.jl")
-include("OctreeBoxPolygon.jl")
+    include("getEdgeInterpolationMatrix.jl")
+    include("getFaceInterpolationMatrix.jl")
+    include("getNodalInterpolationMatrix.jl")
+
+    include("getInterpolationMatrix.jl")
+    include("getEdgeToCellCenteredMatrix.jl")
+    include("getFaceToCellCenteredMatrix.jl")
+    include("getNodesToCellCenteredMatrix.jl")
 
-#include("createSmallMeshFromTX.jl")
-#include("createOcTreeFromTRX.jl")
+    include("findNonRegularBlocks.jl")
+    include("findBlocks.jl")
 
-include("regularizeOcTree2.jl")
-include("splitCells.jl")
-include("getInterfaceWeights.jl")
+    include("getNodalConstraints.jl")
+    include("getEdgeConstraints.jl")
+    include("getFaceConstraints.jl")
 
-#include("getLocalElementMatrices.jl") 
-#include("getMassMatrixFEM.jl") 
-#include("getDiffMassMatrixFEM.jl")
+    include("regularizeOcTree.jl")
+    include("regularizeOcTree2.jl")
 
-# include("plot.jl")
+    include("getNumberOfNeighbors.jl")
+    include("getVolume.jl")
+    include("getLength.jl")
 
+    include("getEdgeIntegralOfPolygonalChain.jl")
+
+    include("refineOcTree.jl")    
+    include("splitCells.jl")
+    include("uniteOcTrees.jl")
+
+    include("getInterfaceWeights.jl")
+
+    # Mesh Creation
+    include("createOcTree/initCoarseOcTree.jl")
+    include("createOcTree/createOcTreeFromBox.jl")
+    include("createOcTree/createOcTreeFromImage.jl")
+    # include("createOcTree/createOcTreeFromPoints.jl")
+    # include("createOcTree/createOcTreeFromTopo.jl")
+    include("createOcTree/createOcTreeMesh.jl")
+    include("createOcTree/OctreeBoxPolygon.jl")
+
+    include("IO/importUBC.jl")
+    include("IO/exportUBC.jl")
 
 end