fix surface vectors on curved meshes

Author: ChrLackner

ngsolve_webgpu/shaders/surface_vectors.wgsl

@@ -14,8 +14,14 @@
 fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
   let n_trigs = mesh.num_trigs;
   for (var trigId = id.x; trigId<n_trigs; trigId+=256*1024) {
-    let p = loadTriangle(trigId).p;
+    var p = loadTriangle(trigId).p;
 
+    if (mesh.is_curved == 1u) {
+      p[0] = evalTrigVec3(&mesh.data, trigId, vec2f(1.0, 0.0), mesh.offset_curvature_2d);
+      p[1] = evalTrigVec3(&mesh.data, trigId, vec2f(0.0, 1.0), mesh.offset_curvature_2d);
+      p[2] = evalTrigVec3(&mesh.data, trigId, vec2f(0.0, 0.0), mesh.offset_curvature_2d);
+    }
+
     let gridsize = u_gridsize;
 
     var dir: u32 =0;
@@ -78,9 +84,6 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
                   let dir_im = v_ri.im / max(val, 1e-10);
 #endif SCALE_BY_VALUE
                   let index = atomicAdd(&count_vectors, 1);
-                  if (mesh.is_curved == 1u) {
-                    cp = evalTrigVec3(&mesh.data, trigId, eval_lam, mesh.offset_curvature_2d);
-                  }
                   cp += 0.5 * gridsize * normalize(n);
 
                   positions[index*3+0] = cp[0];

ngsolve_webgpu/shaders/surface_vectors_count.wgsl

@@ -3,11 +3,48 @@
 @group(0) @binding(21) var<storage, read_write> count_vectors: atomic<u32>;
 @group(0) @binding(31) var<uniform> u_gridsize: f32;
 
+// Inline curved position evaluation using mesh curvature data (de Casteljau).
+// Avoids importing eval/trig which would declare extra storage bindings.
+fn evalCurvedPos(id: u32, lam: vec2f) -> vec3f {
+    let offset_ = mesh.offset_curvature_2d;
+    let order = i32(mesh.data[offset_ + 1]);
+    let ncomp = u32(mesh.data[offset_ + 0]);
+    let ndof = u32((order + 1) * (order + 2) / 2);
+    let stride = ncomp;
+    let offset = offset_ + ndof * id * stride + 3u;
+    let b = vec3f(lam.x, lam.y, 1.0 - lam.x - lam.y);
+    let dy = order + 1;
+
+    var w: array<vec3f, 15>;
+    for (var i: u32 = 0u; i < min(ndof, 15u); i++) {
+        w[i] = vec3f(mesh.data[offset + i * stride],
+                     mesh.data[offset + i * stride + 1],
+                     mesh.data[offset + i * stride + 2]);
+    }
+
+    for (var n = order; n > 0; n--) {
+        var i0 = 0;
+        for (var iy = 0; iy < n; iy++) {
+            for (var ix = 0; ix < n - iy; ix++) {
+                w[i0 + ix] = mat3x3<f32>(w[i0 + ix], w[i0 + ix + 1], w[i0 + ix + dy - iy]) * b;
+            }
+            i0 += dy - iy;
+        }
+    }
+    return w[0];
+}
+
 @compute @workgroup_size(256)
 fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
   let n_trigs = mesh.num_trigs;
   for (var trigId = id.x; trigId<n_trigs; trigId+=256*1024) {
-    let p = loadTriangle(trigId).p;
+    var p = loadTriangle(trigId).p;
+
+    if (mesh.is_curved == 1u) {
+      p[0] = evalCurvedPos(trigId, vec2f(1.0, 0.0));
+      p[1] = evalCurvedPos(trigId, vec2f(0.0, 1.0));
+      p[2] = evalCurvedPos(trigId, vec2f(0.0, 0.0));
+    }
 
     let gridsize = u_gridsize;