Merge branch 'main' into compathelper/new_version/2024-05-25-00-12-48…

…-771-01042506695
jlchan · Jun 28, 2024 · 2c45c94 · 2c45c94
2 parents d343195 + 5f0ff52
commit 2c45c94
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 25 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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))