diff --git a/src/IO/exportUBC.jl b/src/IO/exportUBC.jl new file mode 100644 index 0000000..1704216 --- /dev/null +++ b/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 diff --git a/src/IO/importUBC.jl b/src/IO/importUBC.jl new file mode 100644 index 0000000..31eaacb --- /dev/null +++ b/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 diff --git a/src/JOcTree.jl b/src/JOcTree.jl index dda2c3d..d9f5e94 100644 --- a/src/JOcTree.jl +++ b/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 \ No newline at end of file diff --git a/src/OcTreeMeshFEM.jl b/src/OcTreeMeshFEM.jl deleted file mode 100644 index 094c6a0..0000000 --- a/src/OcTreeMeshFEM.jl +++ /dev/null @@ -1,120 +0,0 @@ -export OcTreeMeshFEM, getOcTreeMeshFEM - -type OcTreeMeshFEM <: OcTreeMesh - S::SparseArray3D - h::Vector{Float64} # (3) cell size - x0::Vector{Float64} - n::Vector{Int64} - nc::Int # number of cells - nf::Vector{Int} # (3) number of faces - ne::Vector{Int} # (3) number of edges - Af::SparseMatrixCSC - Ae::SparseMatrixCSC - V::SparseMatrixCSC - N::SparseMatrixCSC # nullspace matrix - FX::SparseArray3D # X face size - FY::SparseArray3D # Y face size - FZ::SparseArray3D # Z face size - EX::SparseArray3D # X edge size - EY::SparseArray3D # Y edge size - EZ::SparseArray3D # Z edge size - NFX::SparseArray3D - NFY::SparseArray3D - NFZ::SparseArray3D - NEX::SparseArray3D - NEY::SparseArray3D - NEZ::SparseArray3D -end # type OcTreeMeshFEM - - -function getOcTreeMeshFEM(S,h;x0=zeros(3)) - - empt = spzeros(0,0) - - # get number of cells - i, = find3(S) - nc = size(i,1) - - # get number of faces - FX,FY,FZ = getFaceSize(S) - iex, = find3(FX) - iey, = find3(FY) - iez, = find3(FZ) - nf = [size(iex,1); size(iey,1); size(iez,1)] - - # get number of edges - EX,EY,EZ = getEdgeSize(S) - iex, = find3(EX) - iey, = find3(EY) - iez, = find3(EZ) - ne = [size(iex,1); size(iey,1); size(iez,1)] - - empt3 = sparse3([size(S,1),size(S,2),size(S,3)]) - - return OcTreeMeshFEM(S,h,x0, - S.sz,nc,nf,ne, - empt,empt, empt,empt, # no Af,Ae, V,N - FX,FY,FZ, EX,EY,EZ, - empt3,empt3,empt3, # no NFX,NFY,NFZ - empt3,empt3,empt3) # no NEX,NEY,NEZ - -end # function getOcTreeMeshFEM - -import Base.clear! -function clear!(M::OcTreeMeshFEM) - M.S = clear!(M.S) - M.Af = clear!(M.Af ) - M.Ae = clear!(M.Ae ) - M.V = clear!(M.V ) - M.FX = clear!(M.FX ) - M.FY = clear!(M.FY ) - M.FZ = clear!(M.FZ ) - M.N = clear!(M.N ) - M.EX = clear!(M.EX ) - M.EY = clear!(M.EY ) - M.EZ = clear!(M.EZ ) - M.NFX = clear!(M.NFX) - M.NFY = clear!(M.NFY) - M.NFZ = clear!(M.NFZ) - M.NEX = clear!(M.NEX) - M.NEY = clear!(M.NEY) - M.NEZ = clear!(M.NEZ) - -end # function clear - - - -import Base.== -function ==(M1::OcTreeMeshFEM, M2::OcTreeMeshFEM) - isEqual = fill(true,21) - - - # check mandatory fields - isEqual[1] = M1.S==M2.S - isEqual[2] = M1.h==M2.h - isEqual[3] = (M1.x0 == M2.x0) - isEqual[5] = (M1.n == M2.n) - isEqual[6] = (M1.nc == M2.nc) - isEqual[7] = (M1.nf == M2.nf) - isEqual[8] = (M1.ne == M2.ne) - - # check fields that might be empty - if !(isempty(M1.Af)) && !(isempty(M2.Af)) - isEqual[12] = (M1.Af == M2.Af) - end - if !(isempty(M1.Ae)) && !(isempty(M2.Ae)) - isEqual[13] = (M1.Ae == M2.Ae) - end - if !(isempty(M1.V)) && !(isempty(M2.V)) - isEqual[15] = (M1.V == M2.V) - end - if !(isempty(M1.N)) && !(isempty(M2.N)) - isEqual[21] = (M1.N == M2.N) - end - - return all(isEqual) -end # function == - -include("getMatricesFEM.jl") - - diff --git a/src/OctreeBoxPolygon.jl b/src/createOcTree/OctreeBoxPolygon.jl similarity index 100% rename from src/OctreeBoxPolygon.jl rename to src/createOcTree/OctreeBoxPolygon.jl diff --git a/src/createOcTreeFromBox.jl b/src/createOcTree/createOcTreeFromBox.jl similarity index 100% rename from src/createOcTreeFromBox.jl rename to src/createOcTree/createOcTreeFromBox.jl diff --git a/src/createOcTreeFromImage.jl b/src/createOcTree/createOcTreeFromImage.jl similarity index 100% rename from src/createOcTreeFromImage.jl rename to src/createOcTree/createOcTreeFromImage.jl diff --git a/src/createOcTreeFromPoints.jl b/src/createOcTree/createOcTreeFromPoints.jl similarity index 100% rename from src/createOcTreeFromPoints.jl rename to src/createOcTree/createOcTreeFromPoints.jl diff --git a/src/createOcTreeFromTopo.jl b/src/createOcTree/createOcTreeFromTopo.jl similarity index 100% rename from src/createOcTreeFromTopo.jl rename to src/createOcTree/createOcTreeFromTopo.jl diff --git a/src/createOcTreeMesh.jl b/src/createOcTree/createOcTreeMesh.jl similarity index 100% rename from src/createOcTreeMesh.jl rename to src/createOcTree/createOcTreeMesh.jl diff --git a/src/initCoarseOcTree.jl b/src/createOcTree/initCoarseOcTree.jl similarity index 100% rename from src/initCoarseOcTree.jl rename to src/createOcTree/initCoarseOcTree.jl diff --git a/src/createOcTreeFromSrcRec.jl b/src/createOcTreeFromSrcRec.jl deleted file mode 100644 index 4d295f4..0000000 --- a/src/createOcTreeFromSrcRec.jl +++ /dev/null @@ -1,39 +0,0 @@ -export createOcTreeFromSrcRec - -function createOcTreeFromSrcRec(x0,n,src,rec,h,w,sigback) -# createOcTreeFromSrcRec(src,rec,w,sigback) - - f = 2*pi*w - # compute skin depth - sd = 500*sqrt(1/sigback/f) - hRec = sd/5 - if minimum(h) > hRec - warning("use smaller h, skin depth is",sd) - end - - # define the box - xmin = minimum([src[:,1];rec[:,1];rec[:,4]]) - xmax = maximum([src[:,1];rec[:,1];rec[:,4]]) - ymin = minimum([src[:,2];rec[:,2];rec[:,5]]) - ymax = maximum([src[:,2];rec[:,2];rec[:,5]]) - zmin = minimum([src[:,3];rec[:,3];rec[:,6]]) - zmax = maximum([src[:,3];rec[:,3];rec[:,6]]) - - S = createOcTreeFromBox( - x0[1], x0[2],x0[3], - n[1],n[2],n[3], - h[1],h[2],h[3], - xmin, xmax, - ymin, ymax, - zmin, zmax, - 4,4) - - S = regularizeOcTree(S) - Msh = getOcTreeMeshFV(S,h) - - EX,EY,EZ = getEdgeNumbering(MFor.S) - Smat = getEdgeIntegralOfPolygonalChain(Msh,src) - Rmat = getEdgeIntegralOfPolygonalChain(Msh,rec) - - return Msh, Smat, Rmat -end \ No newline at end of file diff --git a/src/createOcTreeFromTRX.jl b/src/createOcTreeFromTRX.jl deleted file mode 100644 index be17431..0000000 --- a/src/createOcTreeFromTRX.jl +++ /dev/null @@ -1,56 +0,0 @@ -export getOcTreeFromTRX - -function getOcTreeFromTRX(trx::Array{Transmitter,1},cellSize,padding,depthFine;numFineLayers=2,numCoarseLayers=1,outFile=[]) - # trx = readdata(data_locations.txt; only_loc = true) - #trx = readdata(dataFile) - - # find bounding box for survey area - xmin = Inf - xmax = -Inf - ymin = Inf - ymax = -Inf - zmin = Inf - zmax = -Inf - for trxi in trx - xmin = min(xmin, minimum(trxi.trxpts[1,:]), minimum(trxi.rcvpts[1,:])) - xmax = max(xmax, maximum(trxi.trxpts[1,:]), maximum(trxi.rcvpts[1,:])) - ymin = min(ymin, minimum(trxi.trxpts[2,:]), minimum(trxi.rcvpts[2,:])) - ymax = max(ymax, maximum(trxi.trxpts[2,:]), maximum(trxi.rcvpts[2,:])) - zmin = min(zmin, minimum(trxi.trxpts[3,:]), minimum(trxi.rcvpts[3,:])) - zmax = max(zmax, maximum(trxi.trxpts[3,:]), maximum(trxi.rcvpts[3,:])) - end - - if zmin != zmax - error("Flat survey required") - end - - hx = cellSize[1] - hy = cellSize[2] - hz = cellSize[3] - - # figure out domain size - nx = nextpow2(ceil(Integer,((xmax - xmin) + 2.0 * padding) / hx)) - ny = nextpow2(ceil(Integer,((ymax - ymin) + 2.0 * padding) / hy)) - nz = nextpow2(ceil(Integer,((zmax - zmin) + 2.0 * padding) / hz)) - x0 = 0.5 * (xmin + xmax - nx * hx) - y0 = 0.5 * (ymin + ymax - ny * hy) - z0 = 0.5 * (zmin + zmax - nz * hz) - - # add depth of fine cells - zmin -= depthFine - - # create Octree mesh from box - S = createOcTreeFromBox( - x0, y0, z0, nx, ny, nz, hx, hy, hz, - xmin, xmax, ymin, ymax, zmin, zmax, - numFineLayers, numCoarseLayers) - - M = getOcTreeMeshFV(S, [hx, hy, hz]; x0 = [x0, y0, z0]) - - if !isempty(outFile) - exportOcTreeRoman(outFile, M) - end - return M - -end - diff --git a/src/createQuadTreeFromImage.jl b/src/createQuadTreeFromImage.jl deleted file mode 100644 index 950eccd..0000000 --- a/src/createQuadTreeFromImage.jl +++ /dev/null @@ -1,49 +0,0 @@ -function createQuadTreeFromImage(A,tol,minsz) -# [S,A] = createOcTreeFromImage(A,tol); -# A - image -# tol - equal intensity color - -m1,m2 = size(A) - - -maxbsz = min((m1,m2)) -S = ones(Uint64,m1,m2) -Amin = copy(A) -Amax = copy(A) - -cnt = 1 - -for bsz = 2.^[0:round(log2(maxbsz))-1] - i,j = find(S.==bsz) - - I = find( (mod(i-1,2*bsz).==0) & (mod(j-1,2*bsz).==0) ) - - if ~isempty[I] - I00 = sub2ind(size(S),i[I] , j[I] ) - I01 = sub2ind(size(S),i[I]+bsz , j[I] ) - I10 = sub2ind(size(S),i[I] , j[I]+bsz) - I11 = sub2ind(size(S),i[I]+bsz , j[I]+bsz) - - Amin[I00] = min([ Amin[I00] , Amin[I01] , Amin[I10] , Amin[I11] ] ,[],2) - Amax[I00] = max([ Amax[I00] , Amax[I01] , Amax[I10] , Amax[I11] ] ,[],2) - - Ic = find( (Amax[I00]-Amin[I00]) .<= tol ); - - S(I00[Ic]) = 2*bsz; - S(I01[Ic]) = 0; - S(I10[Ic]) = 0; - S(I11[Ic]) = 0; - end; - -end - -S = sparse(S); -%% ------------------------------------------------------- -%% ------------------------------------------------------- -%% ------------------------------------------------------- -function ret = isPowerOfTwo(n) -if 2^floor(log2(n)) == n - ret = 1; -else - ret = 0; -end; \ No newline at end of file diff --git a/src/getCellCenteredGrid.jl b/src/getCellCenteredGrid.jl deleted file mode 100644 index 725aefd..0000000 --- a/src/getCellCenteredGrid.jl +++ /dev/null @@ -1,18 +0,0 @@ -export getCellCenteredGrid - - -function getCellCenteredGrid(M::OcTreeMesh) -# Xc = getCellCenteredGrid(M::OcTreeMesh) - -i,j,k,bsz = find3(M.S) -X = i+bsz/2 -Y = j+bsz/2 -Z = k+bsz/2 - -X = (X .- 1) * M.h[1] .+ M.x0[1] -Y = (Y .- 1) * M.h[2] .+ M.x0[2] -Z = (Z .- 1) * M.h[3] .+ M.x0[3] - -return [X Y Z] - -end diff --git a/src/getCellNumbering.jl b/src/getCellNumbering.jl deleted file mode 100644 index dde1437..0000000 --- a/src/getCellNumbering.jl +++ /dev/null @@ -1,16 +0,0 @@ -export getCellNumbering - -function getCellNumbering(S::SparseArray3D) - -sz = collect(size(S.SV)) -colptr = copy(S.SV.colptr) -rowval = copy(S.SV.rowval) - -nz = length(rowval) -nzval = collect(1:nz) - -SV = SparseMatrixCSC(sz[1],sz[2], colptr, rowval, nzval) - -CN = SparseArray3D(SV, S.sz) -return CN -end # function getCellNumbering diff --git a/src/getEdgeGrids.jl b/src/getEdgeGrids.jl deleted file mode 100644 index 6baab80..0000000 --- a/src/getEdgeGrids.jl +++ /dev/null @@ -1,25 +0,0 @@ -export getEdgeGrids - -function getEdgeGrids(M::OcTreeMesh) -# EX,EY,EZ = getEdgeGrids(M::OcTreeMesh) - -EX,EY,EZ = getEdgeSize(M.S) -i,j,k,esz = find3(EX) -EXX = (i+esz/2 .-1) * M.h[1] .+ M.x0[1] -EXY = (j .-1) * M.h[2] .+ M.x0[2] -EXZ = (k .-1) * M.h[3] .+ M.x0[3] - - -i,j,k,esz = find3(EY) -EYX = (i .-1) * M.h[1] .+ M.x0[1] -EYY = (j+esz/2 .-1) * M.h[2] .+ M.x0[2] -EYZ = (k .-1) * M.h[3] .+ M.x0[3] - -i,j,k,esz = find3(EZ) -EZX = (i .-1) * M.h[1] .+ M.x0[1] -EZY = (j .-1) * M.h[2] .+ M.x0[2] -EZZ = (k+esz/2 .-1) * M.h[3] .+ M.x0[3] - -return [EXX EXY EXZ],[EYX EYY EYZ],[EZX EZY EZZ] - -end diff --git a/src/getEdgeNumbering.jl b/src/getEdgeNumbering.jl deleted file mode 100644 index ef6b919..0000000 --- a/src/getEdgeNumbering.jl +++ /dev/null @@ -1,41 +0,0 @@ -export getEdgeNumbering - -function getEdgeNumbering(M::OcTreeMesh) - if isempty(M.NEX) || isempty(M.NEY) || isempty(M.NEZ) - M.NEX,M.NEY,M.NEZ = getEdgeNumbering(M.S) - end - return M.NEX,M.NEY,M.NEZ -end - -function getEdgeNumbering(S) -# [EX,EY,EZ] = getEdgeNumbering(S) -# Numbering of the edges of the structure - -m1,m2,m3 = S.sz -i,j,k,bsz = find3(S) - - -sizeEX = (m1,m2+1,m3+1) -ii = [ i ; i ; i ; i ;] -jj = [ j ; j+bsz ; j ; j+bsz ;] -kk = [ k ; k ; k+bsz ; k+bsz ;] -ii,jj,kk = ind2sub(sizeEX,sort(unique(sub2ind(sizeEX,ii,jj,kk)))) # make em unique -EX = sparse3(ii,jj,kk,1:length(ii), collect(sizeEX)) - -sizeEY = (m1+1,m2,m3+1) -ii = [ i ; i+bsz ; i ; i+bsz ;] -jj = [ j ; j ; j ; j ;] -kk = [ k ; k ; k+bsz ; k+bsz ;] -ii,jj,kk = ind2sub(sizeEY,sort(unique(sub2ind(sizeEY,ii,jj,kk)))) # make em unique -EY = sparse3(ii,jj,kk,1:length(ii), collect(sizeEY)) - -sizeEZ = (m1+1,m2+1,m3) -ii = [ i ; i+bsz ; i ; i+bsz ;] -jj = [ j ; j ; j+bsz ; j+bsz ;] -kk = [ k ; k ; k ; k ;] -ii,jj,kk = ind2sub(sizeEZ,sort(unique(sub2ind(sizeEZ,ii,jj,kk)))) # make em unique -EZ = sparse3(ii,jj,kk,1:length(ii), collect(sizeEZ)) - -return EX, EY, EZ - -end diff --git a/src/getFaceGrids.jl b/src/getFaceGrids.jl deleted file mode 100644 index 9e0c075..0000000 --- a/src/getFaceGrids.jl +++ /dev/null @@ -1,25 +0,0 @@ -export getFaceGrids - -function getFaceGrids(M::OcTreeMesh) -# FX,FY,FZ = getFaceGrids(M::OcTreeMesh) - -FX,FY,FZ = getFaceSize(M.S) -i,j,k,fsz = find3(FX) -FXX = (i .-1) * M.h[1] .+ M.x0[1] -FXY = (j+fsz/2 .-1) * M.h[2] .+ M.x0[2] -FXZ = (k+fsz/2 .-1) * M.h[3] .+ M.x0[3] - - -i,j,k,fsz = find3(FY) -FYX = (i+fsz/2 .-1) * M.h[1] .+ M.x0[1] -FYY = (j .-1) * M.h[2] .+ M.x0[2] -FYZ = (k+fsz/2 .-1) * M.h[3] .+ M.x0[3] - -i,j,k,fsz = find3(FZ) -FZX = (i+fsz/2 .-1) * M.h[1] .+ M.x0[1] -FZY = (j+fsz/2 .-1) * M.h[2] .+ M.x0[2] -FZZ = (k .-1) * M.h[3] .+ M.x0[3] - -return [FXX FXY FXZ],[FYX FYY FYZ],[FZX FZY FZZ] - -end diff --git a/src/getFaceNumbering.jl b/src/getFaceNumbering.jl deleted file mode 100644 index 8723f5c..0000000 --- a/src/getFaceNumbering.jl +++ /dev/null @@ -1,56 +0,0 @@ -export getFaceNumbering - -function getFaceNumbering(M::OcTreeMesh) - if isempty(M.NFX) || isempty(M.NFY) || isempty(M.NFZ) - M.NFX,M.NFY,M.NFZ = getFaceNumbering(M.S) - end - return M.NFX,M.NFY,M.NFZ -end - - - -function getFaceNumbering(S) - -m1,m2,m3 = S.sz -i,j,k,bsz = find3(S) - -## mark upper mark lower -## | | -## v v -ii = vcat( i , i+bsz ) -jj = vcat( j , j ) -kk = vcat( k , k ) - -sizeFX = (m1+1, m2, m3) -I = sort(unique(sub2ind(sizeFX, ii,jj,kk))) ## create unique sorted linear indices -ii,jj,kk = ind2sub(sizeFX, I) ## linear indices to nd indices -FX = sparse3(ii,jj,kk,1:length(ii), collect(sizeFX)) - - -## mark left mark right -## | | -## v v -ii = vcat( i , i ) -jj = vcat( j , j+bsz) -kk = vcat( k , k ) -sizeFY = (m1, m2+1, m3) -I = sort(unique(sub2ind(sizeFY, ii,jj,kk))) ## create unique sorted linear indices -ii,jj,kk = ind2sub(sizeFY, I) ## linear indices to nd indices -FY = sparse3(ii,jj,kk,1:length(ii), collect(sizeFY)) - - -## mark front mark back -## | | -## v v -ii = vcat( i , i ) -jj = vcat( j , j ) -kk = vcat( k , k+bsz ) - -sizeFZ = (m1, m2, m3+1) -I = sort(unique(sub2ind(sizeFZ, ii,jj,kk))) ## create unique sorted linear indices -ii,jj,kk = ind2sub(sizeFZ, I) ## linear indices to nd indices -FZ = sparse3(ii,jj,kk,1:length(ii), collect(sizeFZ)) - -return FX, FY, FZ - -end # function getFaceNumbering diff --git a/src/getGrids.jl b/src/getGrids.jl new file mode 100644 index 0000000..f6e71c9 --- /dev/null +++ b/src/getGrids.jl @@ -0,0 +1,75 @@ +export getCellCenteredGrid, getNodalGrid, getEdgeGrids, getFaceGrids + +function getCellCenteredGrid(M::OcTreeMesh) + # Xc = getCellCenteredGrid(M::OcTreeMesh) + + i,j,k,bsz = find3(M.S) + X = i+bsz/2 + Y = j+bsz/2 + Z = k+bsz/2 + + X = (X .- 1) * M.h[1] .+ M.x0[1] + Y = (Y .- 1) * M.h[2] .+ M.x0[2] + Z = (Z .- 1) * M.h[3] .+ M.x0[3] + + return [X Y Z] +end + +function getNodalGrid(M::OcTreeMesh) + # X = getNodalGrid(M::OcTreeMesh) + + N = getNodalNumbering(M.S) + X,Y,Z = find3(N) + + X = (X .- 1) * M.h[1] .+ M.x0[1] + Y = (Y .- 1) * M.h[2] .+ M.x0[2] + Z = (Z .- 1) * M.h[3] .+ M.x0[3] + + return [X Y Z] +end + +function getEdgeGrids(M::OcTreeMesh) + # EX,EY,EZ = getEdgeGrids(M::OcTreeMesh) + + EX,EY,EZ = getEdgeSize(M.S) + i,j,k,esz = find3(EX) + EXX = (i+esz/2 .-1) * M.h[1] .+ M.x0[1] + EXY = (j .-1) * M.h[2] .+ M.x0[2] + EXZ = (k .-1) * M.h[3] .+ M.x0[3] + + + i,j,k,esz = find3(EY) + EYX = (i .-1) * M.h[1] .+ M.x0[1] + EYY = (j+esz/2 .-1) * M.h[2] .+ M.x0[2] + EYZ = (k .-1) * M.h[3] .+ M.x0[3] + + i,j,k,esz = find3(EZ) + EZX = (i .-1) * M.h[1] .+ M.x0[1] + EZY = (j .-1) * M.h[2] .+ M.x0[2] + EZZ = (k+esz/2 .-1) * M.h[3] .+ M.x0[3] + + return [EXX EXY EXZ],[EYX EYY EYZ],[EZX EZY EZZ] +end + +function getFaceGrids(M::OcTreeMesh) + # FX,FY,FZ = getFaceGrids(M::OcTreeMesh) + + FX,FY,FZ = getFaceSize(M.S) + i,j,k,fsz = find3(FX) + FXX = (i .-1) * M.h[1] .+ M.x0[1] + FXY = (j+fsz/2 .-1) * M.h[2] .+ M.x0[2] + FXZ = (k+fsz/2 .-1) * M.h[3] .+ M.x0[3] + + + i,j,k,fsz = find3(FY) + FYX = (i+fsz/2 .-1) * M.h[1] .+ M.x0[1] + FYY = (j .-1) * M.h[2] .+ M.x0[2] + FYZ = (k+fsz/2 .-1) * M.h[3] .+ M.x0[3] + + i,j,k,fsz = find3(FZ) + FZX = (i+fsz/2 .-1) * M.h[1] .+ M.x0[1] + FZY = (j+fsz/2 .-1) * M.h[2] .+ M.x0[2] + FZZ = (k .-1) * M.h[3] .+ M.x0[3] + + return [FXX FXY FXZ],[FYX FYY FYZ],[FZX FZY FZZ] +end diff --git a/src/getMatricesFEM.jl b/src/getMatricesFEM.jl deleted file mode 100644 index 0574591..0000000 --- a/src/getMatricesFEM.jl +++ /dev/null @@ -1,243 +0,0 @@ -export getMatricesFEM, getLocalElementMatrices, getMassMatrixFEM, getDiffMassMatrixFEM - -function getMatricesFEM(M::OcTreeMeshFEM,sigma) -# K, M, Msig = getMatricesFEM(S,h,sigma) -# - -S = M.S -h = M.h - -Nx,Ny,Nz = getEdgeNumbering(M) - -ii,jj,kk,bsz = find3(S) - -nex = nnz(Nx) -ney = nnz(Ny) -nez = nnz(Nz) -nc = nnz(S) -ne = nex + ney + nez - -Ke,Me,Ge = getLocalElementMatrices(h) - -IE = zeros(nc,12); -# identify edges -IE[:,1] = Nx.SV[sub2ind(size(Nx),ii,jj,kk)] -IE[:,2] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk)] -IE[:,3] = Nx.SV[sub2ind(size(Nx),ii,jj,kk+bsz)] -IE[:,4] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk+bsz)] -IE[:,5] = Ny.SV[sub2ind(size(Ny),ii,jj,kk)] + nex -IE[:,6] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk)] + nex -IE[:,7] = Ny.SV[sub2ind(size(Ny),ii,jj,kk+bsz)] + nex -IE[:,8] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk+bsz)] + nex -IE[:,9] = Nz.SV[sub2ind(size(Nz),ii,jj,kk)] + nex + ney -IE[:,10] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj,kk)] + nex + ney -IE[:,11] = Nz.SV[sub2ind(size(Nz),ii,jj+bsz,kk)] + nex + ney -IE[:,12] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj+bsz,kk)] + nex + ney - -IE = round(Int64,IE) -ii = vec((kron(ones(Int64,1,12),IE))') -jj = vec((kron(IE,ones(Int64,1,12)))') -kc = vec(kron(bsz,Ke[:])) -km = vec(kron(bsz.^3,Me[:])) -ks = vec(kron(sigma.*(bsz.^3),Me[:])) - -K = sparse(ii,jj,vec(full(kc)),ne,ne) -M = sparse(ii,jj,vec(full(km)),ne,ne) - -Msig = sparse(ii,jj,vec(full(ks)),ne,ne) - - -return K, M, Msig - -end - -function getMassMatrixFEM(M::OcTreeMeshFEM,sigma) -# K, M, Msig = getMatricesFEM(S,h,sigma) -# - -S = M.S -h = M.h - -Nx,Ny,Nz = getEdgeNumbering(M) - -ii,jj,kk,bsz = find3(S) - -nex = nnz(Nx) -ney = nnz(Ny) -nez = nnz(Nz) -nc = nnz(S) -ne = nex + ney + nez - -Ke,Me,Ge = getLocalElementMatrices(h) - -IE = zeros(nc,12); -# identify edges -IE[:,1] = Nx.SV[sub2ind(size(Nx),ii,jj,kk)] -IE[:,2] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk)] -IE[:,3] = Nx.SV[sub2ind(size(Nx),ii,jj,kk+bsz)] -IE[:,4] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk+bsz)] -IE[:,5] = Ny.SV[sub2ind(size(Ny),ii,jj,kk)] + nex -IE[:,6] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk)] + nex -IE[:,7] = Ny.SV[sub2ind(size(Ny),ii,jj,kk+bsz)] + nex -IE[:,8] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk+bsz)] + nex -IE[:,9] = Nz.SV[sub2ind(size(Nz),ii,jj,kk)] + nex + ney -IE[:,10] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj,kk)] + nex + ney -IE[:,11] = Nz.SV[sub2ind(size(Nz),ii,jj+bsz,kk)] + nex + ney -IE[:,12] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj+bsz,kk)] + nex + ney - -IE = round(Int64,IE) -ii = vec((kron(ones(Int64,1,12),IE))') -jj = vec((kron(IE,ones(Int64,1,12)))') -ks = vec(kron(sigma.*(bsz.^3),Me[:])) - -Msig = sparse(ii,jj,vec(full(ks)),ne,ne) - -return Msig - -end - - -function getDiffMassMatrixFEM(M::OcTreeMeshFEM,u::Vector) - - - S = M.S; h = M.h - - Nx,Ny,Nz = getEdgeNumbering(M) - - ii,jj,kk,bsz = find3(S) - - nex = nnz(Nx) - ney = nnz(Ny) - nez = nnz(Nz) - nc = nnz(S) - ne = nex + ney + nez - - Ke,Me,Ge = getLocalElementMatrices(h) - - IE = zeros(nc,12) - - # identify edges - IE[:,1] = Nx.SV[sub2ind(size(Nx),ii,jj,kk)] - IE[:,2] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk)] - IE[:,3] = Nx.SV[sub2ind(size(Nx),ii,jj,kk+bsz)] - IE[:,4] = Nx.SV[sub2ind(size(Nx),ii,jj+bsz,kk+bsz)] - IE[:,5] = Ny.SV[sub2ind(size(Ny),ii,jj,kk)] + nex - IE[:,6] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk)] + nex - IE[:,7] = Ny.SV[sub2ind(size(Ny),ii,jj,kk+bsz)] + nex - IE[:,8] = Ny.SV[sub2ind(size(Ny),ii+bsz,jj,kk+bsz)] + nex - IE[:,9] = Nz.SV[sub2ind(size(Nz),ii,jj,kk)] + nex + ney - IE[:,10] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj,kk)] + nex + ney - IE[:,11] = Nz.SV[sub2ind(size(Nz),ii,jj+bsz,kk)] + nex + ney - IE[:,12] = Nz.SV[sub2ind(size(Nz),ii+bsz,jj+bsz,kk)] + nex + ney - - IE = round(Int64,IE) - U = diagm(sparse(bsz.^3))*u[IE] - - ii = vec(IE') - jj = kron([1:nnz(S);],ones(Int64,12)) - vv = vec(Me*U') - - dMsig = sparse(ii,jj,vv,ne,nc) - - return dMsig - -end - -function getLocalElementMatrices(h) -# [K,M] = getElementMatrices(a,b,c) -# Compute stiffness and mass matrix for brick element of size a x b x c. -# -# The edge degrees of freedom are enumerated lexicographically as follows: -# -# x-edges: -# +---4---+ -# /| /| -# +---3---+ | -# | | | | z -# | +--2--|-+ | y -# |/ |/ |/ -# +---1---+ +--- x -# -# y-edges: -# +-------+ -# 7| 8| -# +-------+ | -# | | | | z -# | +-----|-+ | y -# |5 |6 |/ -# +-------+ +--- x -# -# z-edges: -# +-------+ -# /| /| -# +11-----+12 -# | | | | z -# 9 +----10-+ | y -# |/ |/ |/ -# +-------+ +--- x -# - - -a = h[1]; b = h[2]; c = h[3] - -# unit cube edge mass matrix times 36 -Mei = [ - 4 2 2 1 - 2 4 1 2 - 2 1 4 2 - 1 2 2 4] - -Me = kron(speye(3),Mei) - -# unit cube face mass matrix times 6 -Mfi = [ - 2 1 - 1 2] - -Mf = kron(speye(3),Mfi) - -# topological curl operator (edge integral to face integral degrees of -# freedom) -CURL = [ - 0 0 0 0 1 0 -1 0 -1 0 1 0 - 0 0 0 0 0 1 0 -1 0 -1 0 1 - -1 0 1 0 0 0 0 0 1 -1 0 0 - 0 -1 0 1 0 0 0 0 0 0 1 -1 - 1 -1 0 0 -1 1 0 0 0 0 0 0 - 0 0 1 -1 0 0 -1 1 0 0 0 0] - -CURL = sparse(CURL) - -E = diagm(sparse([a a a a b b b b c c c c])) # edge length -F = diagm(sparse(1./[b*c b*c a*c a*c a*b a*b])) # face area -V = a*b*c # cell volume - -# scale by metric -CURL = F * (CURL * E) -Me = Me * (V / 36.0) -Mf = Mf * (V / 6.0) - -# assemble -K = CURL' * Mf * CURL -M = Me - -# Topological Gradient -GRAD = spzeros(12,8) -GRAD[1,[1;2;]] = [-1 1] -GRAD[2,[3;4;]] = [-1 1] -GRAD[3,[5;6;]] = [-1 1] -GRAD[4,[7;8;]] = [-1 1] -GRAD[5,[1;3;]] = [-1 1] -GRAD[6,[2;4;]] = [-1 1] -GRAD[7,[5;7;]] = [-1 1] -GRAD[8,[6;8;]] = [-1 1] -GRAD[9,[1;5;]] = [-1 1] -GRAD[10,[2;6;]] = [-1 1] -GRAD[11,[3;7;]] = [-1 1] -GRAD[12,[4;8;]] = [-1 1] - -GRAD = E\GRAD - -return K,M,GRAD - -end \ No newline at end of file diff --git a/src/getNodalGrid.jl b/src/getNodalGrid.jl deleted file mode 100644 index ce5830e..0000000 --- a/src/getNodalGrid.jl +++ /dev/null @@ -1,14 +0,0 @@ -export getNodalGrid - -function getNodalGrid(M::OcTreeMesh) -# X = getNodalGrid(M::OcTreeMesh) - -N = getNodalNumbering(M.S) -X,Y,Z = find3(N) - -X = (X .- 1) * M.h[1] .+ M.x0[1] -Y = (Y .- 1) * M.h[2] .+ M.x0[2] -Z = (Z .- 1) * M.h[3] .+ M.x0[3] - -return [X Y Z] -end diff --git a/src/getNodalNumbering.jl b/src/getNodalNumbering.jl deleted file mode 100644 index a476c1d..0000000 --- a/src/getNodalNumbering.jl +++ /dev/null @@ -1,32 +0,0 @@ -export getNodalNumbering - -function getNodalNumbering(M::OcTreeMesh) - return getNodalNumbering(M.S) -end - - -function getNodalNumbering(S) -# N = getEdgeNumbering(S) -# Numbering of the nodes of an OcTree structure - -m1,m2,m3 = S.sz -i,j,k,bsz = find3(S) - -nind = [i j k; - i j k+bsz; - i j+bsz k; - i j+bsz k+bsz; - i+bsz j k; - i+bsz j k+bsz; - i+bsz j+bsz k; - i+bsz j+bsz k+bsz ] - -Ni = sparse3(nind[:,1],nind[:,2],nind[:,3],nind[:,3],[m1+1,m2+1,m3+1]) - - -i,j,k = find3(Ni) -N = sparse3(i,j,k,1:length(i), [m1+1,m2+1,m3+1]); - -return N - -end diff --git a/src/getNumbering.jl b/src/getNumbering.jl new file mode 100644 index 0000000..6629cff --- /dev/null +++ b/src/getNumbering.jl @@ -0,0 +1,143 @@ + +export getCellNumbering, getEdgeNumbering, getFaceNumbering, getNodalNumbering + +function getCellNumbering(M::OcTreeMesh) + return getCellNumbering(M.S) +end + +function getCellNumbering(S::SparseArray3D) + sz = collect(size(S.SV)) + colptr = copy(S.SV.colptr) + rowval = copy(S.SV.rowval) + + nz = length(rowval) + nzval = collect(1:nz) + + SV = SparseMatrixCSC(sz[1],sz[2], colptr, rowval, nzval) + + CN = SparseArray3D(SV, S.sz) + return CN +end + +# ------------------------------------------------------------ + +function getEdgeNumbering(M::OcTreeMesh) + if isempty(M.NEX) || isempty(M.NEY) || isempty(M.NEZ) + M.NEX,M.NEY,M.NEZ = getEdgeNumbering(M.S) + end + return M.NEX,M.NEY,M.NEZ +end + +function getEdgeNumbering(S) + # [EX,EY,EZ] = getEdgeNumbering(S) + # Numbering of the edges of the structure + + m1,m2,m3 = S.sz + i,j,k,bsz = find3(S) + + + sizeEX = (m1,m2+1,m3+1) + ii = [ i ; i ; i ; i ;] + jj = [ j ; j+bsz ; j ; j+bsz ;] + kk = [ k ; k ; k+bsz ; k+bsz ;] + ii,jj,kk = ind2sub(sizeEX,sort(unique(sub2ind(sizeEX,ii,jj,kk)))) # make em unique + EX = sparse3(ii,jj,kk,1:length(ii), collect(sizeEX)) + + sizeEY = (m1+1,m2,m3+1) + ii = [ i ; i+bsz ; i ; i+bsz ;] + jj = [ j ; j ; j ; j ;] + kk = [ k ; k ; k+bsz ; k+bsz ;] + ii,jj,kk = ind2sub(sizeEY,sort(unique(sub2ind(sizeEY,ii,jj,kk)))) # make em unique + EY = sparse3(ii,jj,kk,1:length(ii), collect(sizeEY)) + + sizeEZ = (m1+1,m2+1,m3) + ii = [ i ; i+bsz ; i ; i+bsz ;] + jj = [ j ; j ; j+bsz ; j+bsz ;] + kk = [ k ; k ; k ; k ;] + ii,jj,kk = ind2sub(sizeEZ,sort(unique(sub2ind(sizeEZ,ii,jj,kk)))) # make em unique + EZ = sparse3(ii,jj,kk,1:length(ii), collect(sizeEZ)) + + return EX, EY, EZ +end + +# ------------------------------------------------------------ + +function getFaceNumbering(M::OcTreeMesh) + if isempty(M.NFX) || isempty(M.NFY) || isempty(M.NFZ) + M.NFX,M.NFY,M.NFZ = getFaceNumbering(M.S) + end + return M.NFX,M.NFY,M.NFZ +end + +function getFaceNumbering(S) + + m1,m2,m3 = S.sz + i,j,k,bsz = find3(S) + + ## mark upper mark lower + ## | | + ## v v + ii = vcat( i , i+bsz ) + jj = vcat( j , j ) + kk = vcat( k , k ) + + sizeFX = (m1+1, m2, m3) + I = sort(unique(sub2ind(sizeFX, ii,jj,kk))) ## create unique sorted linear indices + ii,jj,kk = ind2sub(sizeFX, I) ## linear indices to nd indices + FX = sparse3(ii,jj,kk,1:length(ii), collect(sizeFX)) + + ## mark left mark right + ## | | + ## v v + ii = vcat( i , i ) + jj = vcat( j , j+bsz) + kk = vcat( k , k ) + sizeFY = (m1, m2+1, m3) + I = sort(unique(sub2ind(sizeFY, ii,jj,kk))) ## create unique sorted linear indices + ii,jj,kk = ind2sub(sizeFY, I) ## linear indices to nd indices + FY = sparse3(ii,jj,kk,1:length(ii), collect(sizeFY)) + + ## mark front mark back + ## | | + ## v v + ii = vcat( i , i ) + jj = vcat( j , j ) + kk = vcat( k , k+bsz ) + + sizeFZ = (m1, m2, m3+1) + I = sort(unique(sub2ind(sizeFZ, ii,jj,kk))) ## create unique sorted linear indices + ii,jj,kk = ind2sub(sizeFZ, I) ## linear indices to nd indices + FZ = sparse3(ii,jj,kk,1:length(ii), collect(sizeFZ)) + + return FX, FY, FZ +end + +# ------------------------------------------------------------ + +function getNodalNumbering(M::OcTreeMesh) + return getNodalNumbering(M.S) +end + +function getNodalNumbering(S) + # N = getEdgeNumbering(S) + # Numbering of the nodes of an OcTree structure + + m1,m2,m3 = S.sz + i,j,k,bsz = find3(S) + + nind = [i j k; + i j k+bsz; + i j+bsz k; + i j+bsz k+bsz; + i+bsz j k; + i+bsz j k+bsz; + i+bsz j+bsz k; + i+bsz j+bsz k+bsz ] + + Ni = sparse3(nind[:,1],nind[:,2],nind[:,3],nind[:,3],[m1+1,m2+1,m3+1]) + + i,j,k = find3(Ni) + N = sparse3(i,j,k,1:length(i), [m1+1,m2+1,m3+1]); + + return N +end diff --git a/src/plot.jl b/src/plot.jl deleted file mode 100644 index fd25f0b..0000000 --- a/src/plot.jl +++ /dev/null @@ -1,68 +0,0 @@ -export plot - -if Pkg.installed("MATLAB") == Nothing() - - function missingPackage() - - warn("Install Julia package MATLAB to enable OcTree mesh plotting.") - return - - end - - plot(mesh::OcTreeMesh) = missingPackage() - plot(mesh::OcTreeMesh, u::Array{Int64,1}) = missingPackage() - plot(mesh::OcTreeMesh, u::Array{Float64,1}) = missingPackage() - -else - - using MATLAB - - function plot(mesh::OcTreeMesh; f = 1) - - u = iround(log2(nonzeros(mesh.S)) .+ 1) - - plot(mesh, u; f = f) - - return - - end - - function plot(mesh::OcTreeMesh, u::Array{Int64,1}; f = 1) - - # to avoid incompatibilities between Julia's and Matlab's sparse3 we pass the indices and values - i,j,k,bsz = find3(mesh.S) - n = mesh.n - h = mesh.h - x0 = mesh.x0 - - # set path to Matlab function plotOcTree.m which resides in the same directory like this Julia function - (dname,fname) = splitdir(@__FILE__()) - mxcall(:addpath, 0, dname) - - # plot - mxcall(:plotOcTreeMesh, 0, f, i, j, k, bsz, n, h, x0, u) - - return - - end - - function plot(mesh::OcTreeMesh, u::Array{Float64,1}; f = 1) - - # to avoid incompatibilities between Julia's and Matlab's sparse3 we pass the indices and values - i,j,k,bsz = find3(mesh.S) - n = mesh.n - h = mesh.h - x0 = mesh.x0 - - # set path to Matlab function plotOcTree.m which resides in the same directory like this Julia function - (dname,fname) = splitdir(@__FILE__()) - mxcall(:addpath, 0, dname) - - # plot - mxcall(:plotOcTreeMesh, 0, f, i, j, k, bsz, n, h, x0, u) - - return - - end - -end \ No newline at end of file diff --git a/src/plotOcTreeMesh.m b/src/plotOcTreeMesh.m deleted file mode 100644 index cfb0002..0000000 --- a/src/plotOcTreeMesh.m +++ /dev/null @@ -1,386 +0,0 @@ -function plotOcTreeMesh(f, i, j, k, bsz, n, h, x0, cdata) - -% Julia provides indices as int64 but Matlab requires double -f = double(f); -i = double(i); -j = double(j); -k = double(k); -bsz = double(bsz); -n = double(n); -cdata = double(cdata); - -nx = n(1); -ny = n(2); -nz = n(3); -hx = h(1); -hy = h(2); -hz = h(3); -z0 = x0(3); -y0 = x0(2); -x0 = x0(1); - -clim = sort(cdata(isfinite(cdata))); -clim = clim([1 end]); -c = mean(clim); -d = diff(clim); -if d < eps(single(c)) - if abs(c) < realmin('single') - clim = [-1 1]; - else - clim = c + eps(single(c)) * [-1 1]; - end -end -cmap = jet(320); -cmap = cmap(33:288, :); -clear c d - -% variables accessed by draw and callback functions -vertices = []; -faces = []; -slices = []; -iSlice = []; -nSlices = []; -edgeColor = [0.5 0.5 0.5]; -hPatch = []; -ratio = []; -tfun = []; -iPlane = 0; - -% plot -figure(f) -clf -axis equal -caxis(clim); -colormap(cmap); -box on -set(gca, ... - 'Layer', 'top'); -colorbar; -% We have to deactivate interactive plotting before changing the scroll -% callback. -plotedit off -zoom off -pan off -rotate3d off -datacursormode off -brush off -set(gcf, ... - 'KeyPressFcn', @keyFcn, ... - 'WindowScrollWheelFcn', @scrollFcn, ... - 'Toolbar', 'figure'); - -% checkbox show mesh -uicontrol( ... - 'Style', 'checkbox', ... - 'Units', 'normalized', ... - 'Position', [0.02 0.10 0.24 0.04], ... - 'Value', true, ... - 'HorizontalAlignment', 'left', ... - 'String', 'Show mesh', ... - 'TooltipString', 'Check to show mesh', ... - 'Callback', @meshFcn, ... - 'Tag', 'Mesh'); - -% checkbox aspect ratio -uicontrol( ... - 'Style', 'checkbox', ... - 'Units', 'normalized', ... - 'Position', [0.02 0.06 0.24 0.04], ... - 'Value', true, ... - 'HorizontalAlignment', 'left', ... - 'String', 'Axis equal', ... - 'TooltipString', 'Check for axis equal', ... - 'Callback', @axisFcn, ... - 'Tag', 'AxisEqual'); - -% checkbox vertical axis direction up -uicontrol( ... - 'Style', 'checkbox', ... - 'Units', 'normalized', ... - 'Position', [0.02 0.02 0.24 0.04], ... - 'Value', false, ... - 'HorizontalAlignment', 'left', ... - 'String', 'Reverse y-axis', ... - 'TooltipString', 'Check to reverse y-axis', ... - 'Callback', @yAxisFcn, ... - 'Tag', 'YDir'); - -% popup for x|y|z plane -h = uicontrol( ... - 'Style', 'popup', ... - 'Units', 'normalized', ... - 'Position', [0.86 0.02 0.12 0.04], ... - 'HorizontalAlignment', 'left', ... - 'String', {'x','y','z'}, ... - 'TooltipString', 'Select slice direction', ... - 'Value', 3, ... - 'Callback', @planeFcn, ... - 'Tag','Plane'); - -% finally, trigger planeFcn to draw -feval(get(h, 'Callback'), h); - -% nested callback functions - function draw - - cla - hPatch = patch( ... - 'Vertices', vertices, ... - 'Faces', faces(slices{iSlice}, :), ... - 'FaceVertexCData', cdata(slices{iSlice}), ... - 'FaceColor', 'flat', ... - 'EdgeColor', edgeColor, ... - 'ButtonDownFcn', @resetFcn); - title(tfun(iSlice)); - drawnow - - end - - function keyFcn(hObject, hEvent) - - if hObject ~= gcf - figure(hObject); - end - - switch hEvent.Key - case 'uparrow' - iSlice = min(iSlice + 1, nSlices); - case 'downarrow' - iSlice = max(1, iSlice - 1); - otherwise - return - end - - draw; - - end - - function scrollFcn(hObject, hEvent) - - if hObject ~= gcf - figure(hObject); - end - - inc = hEvent.VerticalScrollCount; - iSlice = min(max(1, iSlice + inc), nSlices); - draw; - - end - - function meshFcn(hObject, ~, ~) - - if get(hObject, 'Value') == get(hObject, 'Max') - edgeColor = [0.5 0.5 0.5]; - else - edgeColor = 'none'; - end - set(hPatch, 'EdgeColor', edgeColor); - - end - - function axisFcn(hObject, ~, ~) - - if get(hObject, 'Value') == get(hObject, 'Max') - set(gca, 'DataAspectRatio', [1 1 1]); - else - set(gca, 'DataAspectRatio', [ratio 1]); - end - - - end - - function resetFcn(hObject, ~) - - % hObj is handle to patch. 'SelectionType' is defined for figure that - % is two levels up from patch. - type = get(get(get(hObject, 'Parent'), 'Parent'), 'SelectionType'); - - if strcmp(type, 'open') - - iSlice = floor((nSlices + 1) / 2); - draw; - - end - - end - - function planeFcn(hObject, ~) - - jPlane = get(hObject, 'Value'); - if iPlane == jPlane - % same plane, nothing to do - return - end - - % get current axis limits - alim = axis; - switch iPlane - case 0 % initialization - xmin = x0; - xmax = x0 + nx * hx; - ymin = y0; - ymax = y0 + ny * hy; - zmin = z0; - zmax = z0 + nz * hz; - case 1 - ymin = alim(1); - ymax = alim(2); - zmin = alim(3); - zmax = alim(4); - dy = ymax - ymin; - dz = zmax - zmin; - dx = max(dy, dz); - xc = x0 + (iSlice - 0.5) * hx; - xmin = xc - 0.5 * dx; - xmax = xc + 0.5 * dx; - case 2 - xmin = alim(1); - xmax = alim(2); - zmin = alim(3); - zmax = alim(4); - dx = xmax - xmin; - dz = zmax - zmin; - dy = max(dx, dz); - yc = y0 + (iSlice - 0.5) * hy; - ymin = yc - 0.5 * dy; - ymax = yc + 0.5 * dy; - case 3 - xmin = alim(1); - xmax = alim(2); - ymin = alim(3); - ymax = alim(4); - dx = xmax - xmin; - dy = ymax - ymin; - dz = max(dx, dy); - zc = z0 + (iSlice - 0.5) * hz; - zmin = zc - 0.5 * dz; - zmax = zc + 0.5 * dz; - end - - iPlane = jPlane; - - hYDir = findobj(gcf, 'Type', 'uicontrol', 'Tag', 'YDir'); - - switch iPlane - - case 1 - - vertices = [ - j k - j+bsz k - j+bsz k+bsz - j k+bsz]; - [vertices, ~, faces] = unique(vertices, 'rows'); - faces = reshape(faces, [], 4); - - vertices(:,1) = y0 + (vertices(:,1) - 1) * hy; - vertices(:,2) = z0 + (vertices(:,2) - 1) * hz; - - xlim([ymin ymax]); - ylim([zmin zmax]); - xlabel('y in m'); - ylabel('z in m'); - tfun = @(i)(sprintf('x = %gm', x0 + (i - 0.5) * hx)); - % update checkbox for vertical axis - set(hYDir, ... - 'String', 'Reverse z-axis', ... - 'TooltipString', 'Check to reverse z-axis'); - ratio = [hy hz]; - - m = i; - nSlices = nx; - iSlice = round((0.5 * (xmin + xmax) - x0) / hx + 0.5); - iSlice = min(max(1, iSlice), nSlices); - - case 2 - - vertices = [ - i k - i+bsz k - i+bsz k+bsz - i k+bsz]; - [vertices, ~, faces] = unique(vertices, 'rows'); - faces = reshape(faces, [], 4); - - vertices(:,1) = x0 + (vertices(:,1) - 1) * hx; - vertices(:,2) = z0 + (vertices(:,2) - 1) * hz; - - xlim([xmin xmax]); - ylim([zmin zmax]); - xlabel('x in m'); - ylabel('z in m'); - tfun = @(j)(sprintf('y = %gm', y0 + (j - 0.5) * hy)); - % update checkbox for vertical axis - set(hYDir, ... - 'String', 'Reverse z-axis', ... - 'TooltipString', 'Check to reverse z-axis'); - ratio = [hx hz]; - - m = j; - nSlices = ny; - iSlice = round((0.5 * (ymin + ymax) - y0) / hy + 0.5); - iSlice = min(max(1, iSlice), nSlices); - - case 3 - - vertices = [ - i j - i+bsz j - i+bsz j+bsz - i j+bsz]; - [vertices, ~, faces] = unique(vertices, 'rows'); - faces = reshape(faces, [], 4); - - vertices(:,1) = x0 + (vertices(:,1) - 1) * hx; - vertices(:,2) = y0 + (vertices(:,2) - 1) * hy; - - xlim([xmin xmax]); - ylim([ymin ymax]); - xlabel('x in m'); - ylabel('y in m'); - tfun = @(k)(sprintf('z = %gm', z0 + (k - 0.5) * hz)); - % update checkbox for vertical axis - set(hYDir, ... - 'String', 'Reverse y-axis', ... - 'TooltipString', 'Check to reverse y-axis'); - ratio = [hx hy]; - - m = k; - nSlices = nz; - iSlice = round((0.5 * (zmin + zmax) - z0) / hz + 0.5); - iSlice = min(max(1, iSlice), nSlices); - - end - - bval = reshape(unique(bsz), 1, []); - nbsz = length(bval); - - slices = cell(nSlices, 1); - - for ibsz = 1:nbsz - - b = bval(ibsz); - M = find(bsz == b); - N = repmat(m(M), 1, b) + repmat(0:b-1, length(M), 1); - M = repmat(M, 1, b); - - s = accumarray(N(:), M(:), [nSlices, 1], @(x)({x})); - slices = cellfun(@(x,y)([x;y]), slices, s, 'UniformOutput', false); - - end - - draw; - - end - - function yAxisFcn(hObject, ~, ~) - - if get(hObject, 'Value') == get(hObject, 'Max') - set(gca, 'YDir', 'reverse'); - else - set(gca, 'YDir', 'normal'); - end - - end - -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a3aeb11..e51e4a2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ + using JOcTree using Base.Test -# write your own tests here -@test 1 == 1 +println("==== test input & output ====") +include("testIO.jl") diff --git a/test/testIO.jl b/test/testIO.jl new file mode 100644 index 0000000..67fea94 --- /dev/null +++ b/test/testIO.jl @@ -0,0 +1,19 @@ + +S = initializeOctree([256, 256, 256]) +for i in 1:4 + refineInd = unique(rand(1:nnz(S),5)) + S = splitCells(S, refineInd) +end +meshOut = getOcTreeMeshFV(S, rand(1.:100.0,3); x0=rand(1.:100.0,3)) +modelOut = rand(1:0.1:100, meshOut.nc) +exportUBCOcTreeMesh("mesh.msh", meshOut) +exportUBCOcTreeModel("model.mod", meshOut, modelOut) + +meshIn = importUBCOcTreeMesh("mesh.msh") +modelIn = importUBCOcTreeModel("model.mod", meshIn); + +@test meshIn==meshOut +@test modelIn==modelOut + +rm("mesh.msh") +rm("model.mod")