fix surface vectors

Author: ChrLackner

ngsolve_webgpu/shaders/surface_vectors.wgsl

@@ -17,8 +17,6 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
     let p = loadTriangle(trigId).p;
 
     let gridsize = u_gridsize;
-    let pmin = vec3f(-1, -1, -1);
-    let rad = 1.0;
 
     var dir: u32 =0;
     var dir1: u32 =0;
@@ -43,8 +41,7 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
 
     for (var k: u32 = 0; k < 3; k++)
       {
-        p2d[k] = vec2f((p[k][dir1] - pmin[dir1]) / (2*rad),
-                       (p[k][dir2] - pmin[dir2]) / (2*rad));
+        p2d[k] = vec2f(p[k][dir1], p[k][dir2]);
       }
 
     var min2d = min(min(p2d[0], p2d[1]), p2d[2]);
@@ -57,20 +54,20 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
     let mdet = determinant(m);
 
     let minv = 1.0/mdet * mat2x2f( m[1][1], -m[0][1], -m[1][0], m[0][0] );
-    
-    for (var s = 0.0; s <= 1.; s += 1.0 * gridsize) {
-      if (s >= min2d.x && s <= max2d.x) 
-      {
-        for (var t = 0.; t <= 1.; t += 1.0 * gridsize) {
-          if (t >= min2d.y && t <= max2d.y)
-            {
+
+    let s_start = ceil(min2d.x / gridsize) * gridsize;
+    let t_start = ceil(min2d.y / gridsize) * gridsize;
+
+    for (var s = s_start; s <= max2d.x; s += gridsize) {
+        for (var t = t_start; t <= max2d.y; t += gridsize) {
               let lam = minv * (vec2f(s, t) - p2d[0]);
 
               if (lam.x >= 0 && lam.y >= 0 && lam.x+lam.y <= 1)
                 {
-                  var cp = lam.x*p[0] + lam.y * p[1] + (1.0 - lam.x - lam.y) * p[2];
+                  var cp = p[0] + lam.x * (p[1] - p[0]) + lam.y * (p[2] - p[0]);
+                  let eval_lam = vec2f(1.0 - lam.x - lam.y, lam.x);
 
-                  let v_ri = evalTrigVec3ReIm(&u_function_values_2d, trigId, lam, 0u);
+                  let v_ri = evalTrigVec3ReIm(&u_function_values_2d, trigId, eval_lam, 0u);
                   let val = sqrt(dot(v_ri.re, v_ri.re) + dot(v_ri.im, v_ri.im));
 
 #ifdef SCALE_BY_VALUE
@@ -82,7 +79,7 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
 #endif SCALE_BY_VALUE
                   let index = atomicAdd(&count_vectors, 1);
                   if (mesh.is_curved == 1u) {
-                    cp = evalTrigVec3(&mesh.data, trigId, lam, mesh.offset_curvature_2d);
+                    cp = evalTrigVec3(&mesh.data, trigId, eval_lam, mesh.offset_curvature_2d);
                   }
                   cp += 0.5 * gridsize * normalize(n);
 
@@ -99,9 +96,7 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
                   directions_imag[index*3+2] = dir_im[2];
 #endif IS_COMPLEX
                 }
-  }
-    }
-      }
+        }
     }
   }
 }

ngsolve_webgpu/shaders/surface_vectors_count.wgsl

@@ -10,8 +10,6 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
     let p = loadTriangle(trigId).p;
 
     let gridsize = u_gridsize;
-    let pmin = vec3f(-1, -1, -1);
-    let rad = 1.0;
 
     var dir: u32 =0;
     var dir1: u32 =0;
@@ -36,8 +34,7 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
 
     for (var k: u32 = 0; k < 3; k++)
       {
-        p2d[k] = vec2f((p[k][dir1] - pmin[dir1]) / (2*rad),
-                       (p[k][dir2] - pmin[dir2]) / (2*rad));
+        p2d[k] = vec2f(p[k][dir1], p[k][dir2]);
       }
 
     var min2d = min(min(p2d[0], p2d[1]), p2d[2]);
@@ -50,22 +47,19 @@ fn compute_surface_vectors(@builtin(global_invocation_id) id: vec3<u32>) {
     let mdet = determinant(m);
 
     let minv = 1.0/mdet * mat2x2f( m[1][1], -m[0][1], -m[1][0], m[0][0] );
-    
-    for (var s = 0.0; s <= 1.; s += 1.0 * gridsize) {
-      if (s >= min2d.x && s <= max2d.x) 
-      {
-        for (var t = 0.; t <= 1.; t += 1.0 * gridsize) {
-          if (t >= min2d.y && t <= max2d.y)
-            {
+
+    let s_start = ceil(min2d.x / gridsize) * gridsize;
+    let t_start = ceil(min2d.y / gridsize) * gridsize;
+
+    for (var s = s_start; s <= max2d.x; s += gridsize) {
+        for (var t = t_start; t <= max2d.y; t += gridsize) {
               let lam = minv * (vec2f(s, t) - p2d[0]);
 
               if (lam.x >= 0 && lam.y >= 0 && lam.x+lam.y <= 1)
                 {
                   atomicAdd(&count_vectors, 1);
                 }
-  }
-    }
-      }
+        }
     }
   }
 }

ngsolve_webgpu/vectors.py

@@ -76,7 +76,11 @@ def __init__(
         self.scale_by_value = scale_by_value
         mesh = function_data.mesh_data
         bbox = mesh.get_bounding_box()
-        self.box_size = np.linalg.norm(np.array(bbox[1]) - np.array(bbox[0]))
+        # Use the longest side of the bounding box (not its diagonal) so that
+        # thin slabs and cubes of similar lateral extent get a similar sampling
+        # density. The diagonal would couple the spacing to the aspect ratio.
+        sides = np.array(bbox[1]) - np.array(bbox[0])
+        self.box_size = float(np.max(sides))
         self.set_grid_size(grid_size)
         self.__buffers = {}
         self.symmetry = symmetry
@@ -177,6 +181,7 @@ def compute_vectors(self):
             label="n_vectors",
             usage=BufferUsage.STORAGE | BufferUsage.COPY_DST | BufferUsage.COPY_SRC,
             reuse=self.u_nvectors,
+            use_cache=False,
         )
         self.n_vectors = 1
         self.allocate_buffers()
@@ -396,6 +401,7 @@ def __init__(self, function_data: FunctionData, grid_size: float = 20, clipping:
         self.compute_entry_point = "compute_surface_vectors"
         self.u_ntrigs = None
         super().__init__(function_data=function_data, grid_size=grid_size, clipping=clipping, colormap=colormap, symmetry=symmetry, vector_symmetry=vector_symmetry, scale_by_value=scale_by_value)
+        self.gpu_objects.clipping = self.clipping
 
     def update(self, options):
         self.n_search_els = self.function_data.mesh_data.ngs_mesh.GetNE(ngs.BND)
@@ -463,6 +469,7 @@ def compute_vectors(self):
             label="n_vectors",
             usage=BufferUsage.STORAGE | BufferUsage.COPY_DST | BufferUsage.COPY_SRC,
             reuse=self.u_nvectors,
+            use_cache=False,
         )
         self.n_vectors = 1
         self.allocate_buffers()