From fd9dc7ca98d0ae9be51163faa0250336df959f49 Mon Sep 17 00:00:00 2001
From: jpekkila <johannes.pekkila@aalto.fi>
Date: Tue, 15 Oct 2019 17:37:57 +0300
Subject: [PATCH] Added periodic boundconds to utils

---
 src/utils/memory.c | 65 ++++++++++++++++++++++++++++++++++++++++++++++
 src/utils/memory.h |  2 ++
 2 files changed, 67 insertions(+)

diff --git a/src/utils/memory.c b/src/utils/memory.c
index 669bfe8..d74d5a2 100644
--- a/src/utils/memory.c
+++ b/src/utils/memory.c
@@ -74,6 +74,71 @@ acMeshRandomize(AcMesh* mesh)
     return AC_SUCCESS;
 }
 
+AcResult
+acMeshApplyPeriodicBounds(AcMesh* mesh)
+{
+    const AcMeshInfo info = mesh->info;
+    for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
+        const int3 start = (int3){0, 0, 0};
+        const int3 end   = (int3){info.int_params[AC_mx], info.int_params[AC_my],
+                                info.int_params[AC_mz]};
+
+        const int nx = info.int_params[AC_nx];
+        const int ny = info.int_params[AC_ny];
+        const int nz = info.int_params[AC_nz];
+
+        const int nx_min = info.int_params[AC_nx_min];
+        const int ny_min = info.int_params[AC_ny_min];
+        const int nz_min = info.int_params[AC_nz_min];
+
+        // The old kxt was inclusive, but our mx_max is exclusive
+        const int nx_max = info.int_params[AC_nx_max];
+        const int ny_max = info.int_params[AC_ny_max];
+        const int nz_max = info.int_params[AC_nz_max];
+
+        for (int k_dst = start.z; k_dst < end.z; ++k_dst) {
+            for (int j_dst = start.y; j_dst < end.y; ++j_dst) {
+                for (int i_dst = start.x; i_dst < end.x; ++i_dst) {
+
+                    // If destination index is inside the computational domain, return since
+                    // the boundary conditions are only applied to the ghost zones
+                    if (i_dst >= nx_min && i_dst < nx_max && j_dst >= ny_min && j_dst < ny_max &&
+                        k_dst >= nz_min && k_dst < nz_max)
+                        continue;
+
+                    // Find the source index
+                    // Map to nx, ny, nz coordinates
+                    int i_src = i_dst - nx_min;
+                    int j_src = j_dst - ny_min;
+                    int k_src = k_dst - nz_min;
+
+                    // Translate (s.t. the index is always positive)
+                    i_src += nx;
+                    j_src += ny;
+                    k_src += nz;
+
+                    // Wrap
+                    i_src %= nx;
+                    j_src %= ny;
+                    k_src %= nz;
+
+                    // Map to mx, my, mz coordinates
+                    i_src += nx_min;
+                    j_src += ny_min;
+                    k_src += nz_min;
+
+                    const size_t src_idx = acVertexBufferIdx(i_src, j_src, k_src, info);
+                    const size_t dst_idx = acVertexBufferIdx(i_dst, j_dst, k_dst, info);
+                    ERRCHK(src_idx < acVertexBufferSize(info));
+                    ERRCHK(dst_idx < acVertexBufferSize(info));
+                    mesh->vertex_buffer[w][dst_idx] = mesh->vertex_buffer[w][src_idx];
+                }
+            }
+        }
+    }
+    return AC_SUCCESS;
+}
+
 AcResult
 acMeshClear(AcMesh* mesh)
 {
diff --git a/src/utils/memory.h b/src/utils/memory.h
index f955476..6c3409d 100644
--- a/src/utils/memory.h
+++ b/src/utils/memory.h
@@ -39,6 +39,8 @@ AcResult acMeshSet(const AcReal value, AcMesh* mesh);
 
 AcResult acMeshRandomize(AcMesh* mesh);
 
+AcResult acMeshApplyPeriodicBounds(AcMesh* mesh);
+
 AcResult acMeshClear(AcMesh* mesh);
 
 #ifdef __cplusplus