Analytic SolCx: expose exact stress tensor (SolCx_stress)

src/underworld3/function/AnalyticSolCx.c

@@ -43,3 +43,17 @@ double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const dou
     (void) n; (void) z;
     return ( x < x_c ) ? eta_A : eta_B;
 }
+
+tensor2 SolCx_stress( double eta_A, double eta_B, double x_c, int n, const double x, const double z )
+{
+    double pos[2];
+    double total_stress[3];   /* [0]=sigma_xx, [1]=sigma_zz, [2]=sigma_xz (full Cauchy) */
+    tensor2 out;
+    pos[0] = x;
+    pos[1] = z;
+    _Velic_solCx( pos, eta_A, eta_B, x_c, n, NULL, NULL, total_stress, NULL );
+    out.xx = total_stress[0];
+    out.zz = total_stress[1];
+    out.xz = total_stress[2];
+    return out;
+}

src/underworld3/function/AnalyticSolCx.h

@@ -24,8 +24,9 @@
    the result depends solely on the arguments, so redundant calls (e.g. the
    x and z components evaluated separately) can be collapsed by the compiler. */
 
-vec2   SolCx_velocity(  double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
-double SolCx_pressure(  double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
-double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
+vec2    SolCx_velocity(  double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
+double  SolCx_pressure(  double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
+double  SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
+tensor2 SolCx_stress(    double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const));
 
 #endif

src/underworld3/function/analytic.pyx

@@ -12,15 +12,20 @@ cdef extern from "AnalyticSolNL.h" nogil:
     ctypedef struct vec2:
         double x
         double z
+    ctypedef struct tensor2:
+        double xx
+        double zz
+        double xz
     vec2   SolNL_velocity(  double eta0, unsigned n, double r, double x, double y )
     vec2   SolNL_bodyforce( double eta0, unsigned n, double r, double x, double y )
     double SolNL_viscosity( double eta0, unsigned n, double r, double x, double y )
 
 cdef extern from "AnalyticSolCx.h" nogil:
-    # vec2 already declared above (shared typedef via AnalyticSolNL.h)
-    vec2   SolCx_velocity(  double eta_A, double eta_B, double x_c, int n, double x, double z )
-    double SolCx_pressure(  double eta_A, double eta_B, double x_c, int n, double x, double z )
-    double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, double x, double z )
+    # vec2 / tensor2 already declared above (shared typedefs via AnalyticSolNL.h)
+    vec2    SolCx_velocity(  double eta_A, double eta_B, double x_c, int n, double x, double z )
+    double  SolCx_pressure(  double eta_A, double eta_B, double x_c, int n, double x, double z )
+    double  SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, double x, double z )
+    tensor2 SolCx_stress(    double eta_A, double eta_B, double x_c, int n, double x, double z )
 
 class sympy_function_printable(sympy.Function):
     """
@@ -131,6 +136,25 @@ class AnalyticSolCx_viscosity(AnalyticSolCx_base):
         from sympy import sympify
         return sympify(SolCx_viscosity( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ))
 
+# Exact total (Cauchy) stress components sigma_ij from the Velic kernel's
+# total_stress array: xx, zz (= yy in the x-z plane), xz. Dynamic topography on
+# the top boundary (n = y_hat) is -n.sigma.n = -sigma_zz.
+class AnalyticSolCx_stress_xx(AnalyticSolCx_base):
+    _printstr = "SolCx_stress({}).xx"
+    def _eval_evalf(self,prec):
+        from sympy import sympify
+        return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).xx)
+class AnalyticSolCx_stress_yy(AnalyticSolCx_base):
+    _printstr = "SolCx_stress({}).zz"
+    def _eval_evalf(self,prec):
+        from sympy import sympify
+        return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).zz)
+class AnalyticSolCx_stress_xy(AnalyticSolCx_base):
+    _printstr = "SolCx_stress({}).xz"
+    def _eval_evalf(self,prec):
+        from sympy import sympify
+        return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).xz)
+
 
 class SolCx:
     r"""Velic *SolCx* analytic Stokes solution — the canonical free-slip,
@@ -175,6 +199,10 @@ class SolCx:
         self.fn_velocity  = _sp.Matrix([AnalyticSolCx_velocity_x(*p, x, y),
                                         AnalyticSolCx_velocity_y(*p, x, y)])
         self.fn_pressure  = AnalyticSolCx_pressure(*p, x, y)
+        # exact total (Cauchy) stress components (sigma_xx, sigma_yy, sigma_xy)
+        self.fn_stress_xx = AnalyticSolCx_stress_xx(*p, x, y)
+        self.fn_stress_yy = AnalyticSolCx_stress_yy(*p, x, y)
+        self.fn_stress_xy = AnalyticSolCx_stress_xy(*p, x, y)
         # tie-break at x==x_c matches the compiled SolCx_viscosity (eta_B for x>=x_c)
         self.fn_viscosity = _sp.Piecewise((self.eta_A, x < self.x_c), (self.eta_B, True))
         self.fn_bodyforce = _sp.Matrix([_sp.Integer(0),
@@ -195,4 +223,26 @@ class SolCx:
         """Relative L2 error of a velocity MeshVariable against the analytic."""
         import numpy as _np
         ua = self.evaluate_velocity(velocity_var.coords)
-        return float(_np.linalg.norm(velocity_var.data - ua) / _np.linalg.norm(ua))
+        return float(_np.linalg.norm(velocity_var.data - ua) / _np.linalg.norm(ua))
+
+    def evaluate_stress(self, coords):
+        """Exact total (Cauchy) stress at ``coords`` (N×2 array), via the
+        compiled kernel. Returns an (N, 3) array of (sigma_xx, sigma_yy,
+        sigma_xy)."""
+        import numpy as _np
+        cdef Py_ssize_t i
+        cdef double xi, yi
+        cdef tensor2 s
+        out = _np.empty((len(coords), 3), dtype=_np.float64)
+        for i in range(len(coords)):
+            xi = float(coords[i, 0]); yi = float(coords[i, 1])
+            s = SolCx_stress(self.eta_A, self.eta_B, self.x_c, self.n, xi, yi)
+            out[i, 0] = s.xx
+            out[i, 1] = s.zz
+            out[i, 2] = s.xz
+        return out
+
+    def topography_top(self, coords):
+        """Exact dynamic topography -n.sigma.n on the top boundary (n = y_hat),
+        i.e. -sigma_yy, at ``coords`` (N×2 array). Returns a length-N array."""
+        return -self.evaluate_stress(coords)[:, 1]

tests/test_1015_analytic_solcx.py

@@ -56,3 +56,32 @@ def test_solcx_stokes_converges_to_analytic():
     e32 = _solve_solcx(32)
     assert e32 < e16, f"no convergence: e16={e16:.2e} e32={e32:.2e}"
     assert e32 < 1.0e-3, f"poor match to analytic: e32={e32:.2e}"
+
+
+def test_solcx_stress_binding():
+    """The exact total-stress kernel binding returns finite, non-trivial values."""
+    pts = [(0.25, 1.0), (0.7, 0.2), (0.9, 0.6), (0.15, 0.4)]
+    sxx = [float(A.AnalyticSolCx_stress_xx(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts]
+    syy = [float(A.AnalyticSolCx_stress_yy(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts]
+    sxy = [float(A.AnalyticSolCx_stress_xy(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts]
+    for s in sxx + syy + sxy:
+        assert np.isfinite(s)
+    # the stress field is non-trivial (the dynamic-topography signal on the top is -syy)
+    assert max(abs(s) for s in syy) > 1.0e-3
+    assert max(abs(s) for s in sxx) > 1.0e-3
+
+
+def test_solcx_stress_pressure_consistency():
+    """Convention-independent check: in 2D sigma = -p I + tau with tau traceless,
+    so the mean normal stress (sigma_xx + sigma_yy)/2 must equal -p at every
+    point. This validates the new stress wrapper against the kernel's own
+    pressure without assuming any UW3 sign convention."""
+    pts = [(0.2, 0.3), (0.7, 0.8), (0.5, 0.5), (0.9, 0.1), (0.35, 0.65)]
+    for xi, yi in pts:
+        sxx = float(A.AnalyticSolCx_stress_xx(1.0, 1.0e6, 0.5, 1, xi, yi).evalf())
+        syy = float(A.AnalyticSolCx_stress_yy(1.0, 1.0e6, 0.5, 1, xi, yi).evalf())
+        p = float(A.AnalyticSolCx_pressure(1.0, 1.0e6, 0.5, 1, xi, yi).evalf())
+        mean_normal = 0.5 * (sxx + syy)
+        assert abs(mean_normal + p) < 1.0e-10, (
+            f"(sxx+syy)/2 != -p at ({xi},{yi}): {mean_normal:.3e} vs {-p:.3e}"
+        )