diff --git a/Project.toml b/Project.toml index 9711f0fc..8f858656 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StartUpDG" uuid = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f" authors = ["Jesse Chan", "Yimin Lin"] -version = "1.0.1" +version = "1.0.2" [deps] ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" diff --git a/src/MeshData.jl b/src/MeshData.jl index 2b1fdc3e..0e368151 100644 --- a/src/MeshData.jl +++ b/src/MeshData.jl @@ -332,7 +332,7 @@ function MeshData(VX, VY, VZ, EToV, rd::RefElemData{3}; is_periodic=(false, fals x, y, z = (x -> V1 * x[transpose(EToV)]).((VX, VY, VZ)) #Compute connectivity maps: uP = exterior value used in DG numerical fluxes - (; r, s, t, Vf ) = rd + (; r, s, t, Vf) = rd xf, yf, zf = (x -> Vf * x).((x, y, z)) mapM, mapP, mapB = build_node_maps(rd, FToF, (xf, yf, zf)) mapM = reshape(mapM, :, K) diff --git a/src/connectivity_functions.jl b/src/connectivity_functions.jl index 0a9f18bb..aa845354 100644 --- a/src/connectivity_functions.jl +++ b/src/connectivity_functions.jl @@ -166,7 +166,7 @@ end """ make_periodic(md::MeshData{Dim}, is_periodic...) where {Dim} - make_periodic(md::MeshData{Dim}, is_periodic = ntuple(x->true,Dim)) where {Dim} + make_periodic(md::MeshData{Dim}, is_periodic = ntuple(x -> true, Dim)) where {Dim} make_periodic(md::MeshData, is_periodic = true) Returns new MeshData such that the node maps `mapP` and face maps `FToF` are now periodic. @@ -230,7 +230,6 @@ function build_periodic_boundary_maps!(xf, yf, is_periodic_x, is_periodic_y, ymin, ymax = extrema(yc) LX, LY = map((x -> x[2] - x[1]) ∘ extrema, (xf, yf)) - #NODETOL = 100 * max(eps.((LX, LY))...) NODETOL = tol * max(LX, LY) if abs(abs(xmax - xmin) - LX) > NODETOL && is_periodic_x error("periodicity requested in x, but LX = $LX while abs(xmax-xmin) = $(abs(xmax-xmin))") @@ -288,7 +287,6 @@ function build_periodic_boundary_maps!(xf, yf, zf, Flist = 1:length(FToF) Bfaces = findall(vec(FToF) .== Flist) - # xc, yc, zc = compute_boundary_centroids(xf, yf, zf) xb, yb, zb = xf[mapB], yf[mapB], zf[mapB] Nfp = length(xf) ÷ NfacesTotal Nbfaces = length(xb) ÷ Nfp @@ -322,16 +320,16 @@ function build_periodic_boundary_maps!(xf, yf, zf, yfaces = map(x -> x[1], findall(@. (@. abs(yc - ymax) < NODETOL * LY) | (@. abs(yc - ymin) < NODETOL * LY))) zfaces = map(x -> x[1], findall(@. (@. abs(zc - zmax) < NODETOL * LZ) | (@. abs(zc - zmin) < NODETOL * LZ))) - D = zeros(eltype(xb), size(xb,1), size(xb,1)) - ids = zeros(Int, size(xb, 1)) + p = zeros(Int, size(xb, 1)) if is_periodic_x # find matches in x faces for i in xfaces, j in xfaces - if i!=j - if abs(yc[i] - yc[j]) < NODETOL * LY && abs(zc[i] - zc[j]) < NODETOL * LZ && abs(abs(xc[i] - xc[j]) - LX) < NODETOL * LX - # create distance matrix - @. D = abs(yb[:,i] - yb[:,j]') + abs(zb[:,i] - zb[:,j]') - map!(x->x[1], ids, findall(@. D < NODETOL * LY)) - @. mapPB[:,i] = mapMB[ids,j] + if i!=j + if abs(yc[i] - yc[j]) < NODETOL * LY && + abs(zc[i] - zc[j]) < NODETOL * LZ && + abs(abs(xc[i] - xc[j]) - LX) < NODETOL * LX + + match_coordinate_vectors!(p, (yb[:,i], zb[:,i]), (yb[:,j], zb[:,j])) + @. mapPB[:,i] = mapMB[p,j] FToF[Bfaces[i]] = Bfaces[j] end @@ -343,10 +341,12 @@ function build_periodic_boundary_maps!(xf, yf, zf, if is_periodic_y for i in yfaces, j = yfaces if i!=j - if abs(xc[i] - xc[j]) < NODETOL * LX && abs(zc[i] - zc[j]) < NODETOL * LZ && abs(abs(yc[i] - yc[j]) - LY) < NODETOL * LY - @. D = abs(xb[:,i] - xb[:,j]') + abs(zb[:,i] - zb[:,j]') - map!(x->x[1], ids, findall(@. D < NODETOL * LX)) - @. mapPB[:,i] = mapMB[ids,j] + if abs(xc[i] - xc[j]) < NODETOL * LX && + abs(zc[i] - zc[j]) < NODETOL * LZ && + abs(abs(yc[i] - yc[j]) - LY) < NODETOL * LY + + match_coordinate_vectors!(p, (xb[:,i], zb[:,i]), (xb[:,j], zb[:,j])) + @. mapPB[:,i] = mapMB[p,j] FToF[Bfaces[i]] = Bfaces[j] end @@ -358,10 +358,12 @@ function build_periodic_boundary_maps!(xf, yf, zf, if is_periodic_z for i in zfaces, j in zfaces if i!=j - if abs(xc[i] - xc[j]) < NODETOL * LX && abs(yc[i] - yc[j]) < NODETOL * LY && abs(abs(zc[i] - zc[j]) - LZ) < NODETOL * LZ - @. D = abs(xb[:,i] - xb[:,j]') + abs(yb[:,i] - yb[:,j]') - map!(x->x[1], ids, findall(@. D < NODETOL * LX)) - @. mapPB[:,i] = mapMB[ids,j] + if abs(xc[i] - xc[j]) < NODETOL * LX && + abs(yc[i] - yc[j]) < NODETOL * LY && + abs(abs(zc[i] - zc[j]) - LZ) < NODETOL * LZ + + match_coordinate_vectors!(p, (xb[:,i], yb[:,i]), (xb[:,j], yb[:,j])) + @. mapPB[:,i] = mapMB[p,j] FToF[Bfaces[i]] = Bfaces[j] end diff --git a/src/ref_elem_utils.jl b/src/ref_elem_utils.jl index 1fa75801..1049868d 100644 --- a/src/ref_elem_utils.jl +++ b/src/ref_elem_utils.jl @@ -30,13 +30,13 @@ end function map_face_nodes(::Tet, face_nodes...) r, s = face_nodes e = ones(size(r)) - rf = [r; r; -e; r] - sf = [-e; s; r; s] - tf = [s; -(e + r + s); s; -e] + rf = [r; -(e + r + s); -e; r] + sf = [-e; r; r; s] + tf = [s; s; s; -e] return rf, sf, tf end -function init_face_data(elem::Tri; quad_rule_face = gauss_quad(0,0,N)) +function init_face_data(elem::Tri; quad_rule_face = gauss_quad(0, 0, N)) r1D, w1D = quad_rule_face e = ones(size(r1D)) z = zeros(size(r1D))