Add C implementations of shell-patch derivative kernels (WithShell support)

New files provide C equivalents of Fortran diff_new_sh.f90 and kodiss_sh.f90:
- fderivs_sh_c.C: first derivatives in shell (rho, sigma, R) coords
- fdderivs_sh_c.C: second derivatives in shell coords
- fderivs_shc_c.C: shell first derivs + chain rule to Cartesian
- fdderivs_shc_c.C: shell second derivs + chain rule to Cartesian
- kodiss_sh_c.C: Kreiss-Oliger dissipation on shell patches

Also add symmetry_stbd() C implementation and shell fh indexing to share_func.h.
Controlled by USE_CXX_SHELL_KERNELS switch (default 0, experimental).

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

AMSS_NCKU_source/share_func.h

@@ -289,4 +289,84 @@ static inline void symmetry_bd(int ord,
 
     symmetry_bd_impl(ord, ord - 1, extc, func, funcc, SoA);
 }
+
+/*
+ * symmetry_stbd — shell-patch (staggered boundary) ghost fill.
+ *
+ * Fortran: funcc(-ord+1:extc1+ord, -ord+1:extc2+ord, extc3)
+ * Only 2 SoA values (x/y). No z symmetry fill.
+ * Ghost on BOTH positive and negative sides of x and y.
+ * Reflection uses i+2 (skips boundary) instead of i+1.
+ * nx = extc1 + 2*ord, ny = extc2 + 2*ord
+ */
+static inline void symmetry_stbd(int ord,
+                                 const int extc[3],
+                                 const double *func,
+                                 double *funcc,
+                                 const double SoA[2])
+{
+    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
+    const int nx = extc1 + 2 * ord;
+    const int ny = extc2 + 2 * ord;
+    const int sh = ord - 1;
+    const size_t snx = (size_t)nx;
+    const size_t splane = snx * (size_t)ny;
+
+    /* 1) Copy interior: funcc(1:extc1, 1:extc2, 1:extc3) = func */
+    for (int k0 = 0; k0 < extc3; ++k0) {
+        const double *src = func + (size_t)k0 * (size_t)extc2 * (size_t)extc1;
+        const size_t kbase = (size_t)k0 * splane;
+        for (int j0 = 0; j0 < extc2; ++j0) {
+            double *dst = funcc + kbase + (size_t)(sh + j0) * snx + (size_t)sh;
+            const double *s = src + (size_t)j0 * (size_t)extc1;
+            for (int i0 = 0; i0 < extc1; ++i0) dst[i0] = s[i0];
+        }
+    }
+
+    /* 2) x-direction ghost fill */
+    const double s1 = SoA[0];
+    for (int k0 = 0; k0 < extc3; ++k0) {
+        const size_t kbase = (size_t)k0 * splane;
+        for (int j0 = 0; j0 < extc2; ++j0) {
+            const size_t off = kbase + (size_t)(sh + j0) * snx;
+            /* left side: funcc(-i) = funcc(i+2) * s1 */
+            for (int i = 0; i < ord; ++i) {
+                funcc[off + (size_t)(sh - i)]       = funcc[off + (size_t)(sh + i + 2)] * s1;
+                /* right side: funcc(extc1+1+i) = funcc(extc1-1-i) * s1 */
+                funcc[off + (size_t)(sh + extc1 + 1 + i)] = funcc[off + (size_t)(sh + extc1 - 1 - i)] * s1;
+            }
+        }
+    }
+
+    /* 3) y-direction ghost fill */
+    const double s2 = SoA[1];
+    for (int i = 0; i < nx; ++i) {
+        for (int k0 = 0; k0 < extc3; ++k0) {
+            const size_t kbase = (size_t)k0 * splane;
+            /* bottom: funcc(:,-i,:) = funcc(:,i+2,:) * s2 */
+            for (int jj = 0; jj < ord; ++jj) {
+                funcc[kbase + (size_t)(sh - jj) * snx + (size_t)i] =
+                    funcc[kbase + (size_t)(sh + jj + 2) * snx + (size_t)i] * s2;
+                /* top: funcc(:,extc2+1+jj,:) = funcc(:,extc2-1-jj,:) * s2 */
+                funcc[kbase + (size_t)(sh + extc2 + 1 + jj) * snx + (size_t)i] =
+                    funcc[kbase + (size_t)(sh + extc2 - 1 - jj) * snx + (size_t)i] * s2;
+            }
+        }
+    }
+}
+
+/*
+ * Indexing for shell fh buffer: Fortran fh(-ord+1:extc1+ord, -ord+1:extc2+ord, extc3)
+ * C 0-based: ii = iF + ord - 1
+ * nx = extc1 + 2*ord, ny = extc2 + 2*ord
+ */
+static inline size_t idx_fh_stbd(int iF, int jF, int kF, int ord, const int extc[3]) {
+    const int sh = ord - 1;
+    const int nx = extc[0] + 2 * ord;
+    const int ny = extc[1] + 2 * ord;
+    const int ii = iF + sh;
+    const int jj = jF + sh;
+    const int kk = kF - 1;  // Fortran 1-based kF → C 0-based
+    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
+}
 #endif