--- a/src/laser_sampling.py Mon Jun 29 00:20:42 2026 +0300
+++ b/src/laser_sampling.py Sun Oct 04 13:02:26 2026 -0500
@@ -1,12 +1,15 @@
+import basix
import numpy as np
-from dolfinx import fem, geometry
-from dolfinx_access import cell_FunctionSpace_f64
+from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc
from pointsource_pde.dolfinx_extras import get_mass_matrix
+from scipy.sparse import coo_matrix
from scipy.sparse.linalg import svds
from slepc4py import SLEPc
+from dolfinx_access import cell_FunctionSpace_f64
+
class LaserSampling:
"""
@@ -20,7 +23,6 @@
nx,
ny,
beams,
- num_segments=100,
domain_size=[-2, 2, -2, 2],
):
self.V = V
@@ -28,7 +30,10 @@
self.nx, self.ny = nx, ny
self.M = len(beams)
self.domain_size = domain_size
- self.num_segments = num_segments
+
+ # Basix element matching V, used to tabulate phi_h exactly at points
+ # self._element = self._create_basix_element()
+ self._element = self.V.element.basix_element
# Build R matrix (single function does everything)
self.R = self.matrix(beams)
@@ -60,71 +65,133 @@
self._opnorm = s[0] # / sigma
# print(f"Laser sampling opnorm = {self._opnorm:.6f}")
- def matrix(self, beams):
- """Combined: generate beams + segment + cell collision + dof weighting"""
- dofmap = self.V.dofmap
- tdim = self.domain.topology.dim
+ # ------------------------------------------------------------------
+ # Exact per-cell line integration
+ #
+ # A beam is split into the ordered sub-segments where it actually
+ # crosses each element (via Cyrus-Beck clipping against the cell's
+ # polygon), and each sub-segment is integrated exactly using the
+ # 3-point Simpson rule (exact for phi_h of degree <= 2, since phi_h
+ # pulled back along a straight line stays degree <= 2 in the
+ # parameter t). This replaces the previous approach of sampling many
+ # midpoints and splitting weight equally across a cell's dofs.
+ # ------------------------------------------------------------------
- # Geometry trees
-
- bb_tree = geometry.bb_tree(self.domain, tdim)
-
- n_dofs = self.V.dofmap.index_map.size_global
- R = np.zeros((self.M, n_dofs))
+ def _cell_vertices_2d(self, cell_idx):
+ """Physical (x, y) vertices of a cell, ordered CCW around its centroid."""
+ x_geom = self.domain.geometry.x
+ verts = self.domain.geometry.dofmap[cell_idx]
+ pts = x_geom[verts, :2]
+ centroid = pts.mean(axis=0)
+ angles = np.arctan2(pts[:, 1] - centroid[1], pts[:, 0] - centroid[0])
+ return pts[np.argsort(angles)]
- xmin, xmax, ymin, ymax = self.domain_size
+ def _clip_segment_to_cell(self, pts, x, y, eps=1e-12):
+ """
+ Cyrus-Beck clip of the segment x->y against the convex polygon `cell_idx`.
+ Returns (t0, t1) in [0, 1] such that x + t*(y - x), t in [t0, t1], is the
+ part of the segment inside the cell -- or None if there's no overlap.
+ Works for any convex cell (triangle or quad); assumes straight-sided
+ (affine, P1) geometry.
+ """
+ n_verts = len(pts)
+ d = y - x
+ t_lo, t_hi = 0.0, 1.0
+
+ for i in range(n_verts):
+ v0 = pts[i]
+ v1 = pts[(i + 1) % n_verts]
+ edge = v1 - v0
+ normal = np.array([edge[1], -edge[0]]) # outward normal (CCW polygon)
- for beam_idx, (start_pt, end_pt) in enumerate(beams):
- # Segment into midpoints
- vec = end_pt - start_pt
- total_len = np.linalg.norm(vec)
- seg_len = total_len / self.num_segments
+ denom = normal @ d
+ numer = normal @ (v0 - x)
+
+ if abs(denom) < eps:
+ if numer < -eps: # segment parallel to edge, and entirely outside it
+ return None
+ continue # parallel and inside: this edge imposes no constraint
- midpoints = []
- for i in range(self.num_segments):
- t = (i + 0.5) / self.num_segments
- midpoint = start_pt + t * vec
- midpoints.append(midpoint)
- # midpoints = np.array(midpoints, dtype=np.float64) # Force float64
+ t_edge = numer / denom
+ if denom > 0:
+ t_hi = min(t_hi, t_edge)
+ else:
+ t_lo = max(t_lo, t_edge)
+
+ if t_lo > t_hi + eps:
+ return None
+
+ if t_lo > t_hi:
+ return None
+ return t_lo, t_hi
- # Ensure shape (N, 3), C-contiguous, read-only for DOLFINx
- # assert midpoints.shape[1] == 2, (
- # "Expected 2D points, got shape[1]=2 but need padding to 3D"
- # )
- # midpoints = np.pad(
- # midpoints, ((0, 0), (0, 1)), mode="constant"
- # ) # (N, 3) with z=0
- # midpoints = np.ascontiguousarray(midpoints) # C-order
- # midpoints.setflags(write=False) # Read-only requirement
-
- # Find colliding cells
- # cell_candidates = geometry.compute_collisions_points(bb_tree, midpoints)
- # colliding_cells = geometry.compute_colliding_cells(
- # self.domain, cell_candidates, midpoints
- # )
+ def matrix(self, beams):
+ """Combined: exact beam/cell intersection + exact per-cell line integration."""
+ n_dofs = self.V.dofmap.index_map.size_global
+ # R = np.zeros((self.M, n_dofs))
+ rows, cols, data = [], [], []
- # Cell→DOFs→weight accumulation
- # cell_lists = list(colliding_cells.array)
- # for seg_i, cell_idx_raw in enumerate(cell_lists):
- # cell_idx = int(cell_idx_raw) # Convert np.int32 → Python int
- for point in midpoints:
- pad_point = np.array([point[0], point[1], 0.0])
- cell_idx = cell_FunctionSpace_f64(self.V._cpp_object, pad_point)
- if cell_idx >= 0: # Valid cell index
- # Get global DOFs
- cell_dofs_local = dofmap.cell_dofs(cell_idx)
- # Convert local dofs to numpy int32 array (DOLFINx requirement)
- # local_dofs_array = np.asarray(cell_dofs_local, dtype=np.int32)
- # local_dofs_array.setflags(write=False) # Read-only requirement
+ for cell_idx in range(len(self.domain.geometry.dofmap)):
+ pts = self._cell_vertices_2d(cell_idx)
+ for beam_idx, (start_pt, end_pt) in enumerate(beams):
+ clip = self._clip_segment_to_cell(pts, start_pt, end_pt)
+ if clip is not None:
+ t0, t1 = clip
+ if t1 - t0 > 1e-10:
+ beam_v = end_pt - start_pt
+ x0 = start_pt + t0 * beam_v
+ x1 = start_pt + t1 * beam_v
+ midpt = start_pt + (0.5 * (t1 + t0)) * beam_v
+ """
+ Exact int_{x0}^{x1} phi_h ds for phi_h in V restricted to `cell_idx`, as a
+ (global_dofs, row_values) pair such that row_values @ phi_h.x.array[global_dofs]
+ equals the integral. Uses the 3-point Simpson rule (endpoints + midpoint),
+ which is exact since phi_h has degree <= 2 along the line.
+ """
+ x_geom = self.domain.geometry.x
+ verts = self.domain.geometry.dofmap[cell_idx]
+ v0, v1, v2 = (
+ x_geom[verts[0], :2],
+ x_geom[verts[1], :2],
+ x_geom[verts[2], :2],
+ )
+ J = np.column_stack([v1 - v0, v2 - v0])
+ """Map a physical point p to reference coordinates on an affine triangle."""
+ pts_ref = np.array(
+ [
+ np.linalg.solve(J, x0 - v0),
+ np.linalg.solve(J, midpt - v0),
+ np.linalg.solve(J, x1 - v0),
+ ]
+ )
+ tab = self._element.tabulate(0, pts_ref)[
+ 0
+ ] # (npoints, ndof_local, value_size)
+ tab = tab[
+ ..., 0
+ ] # scalar element -> drop trailing value_size=1 axis
- # Batch convert all local indices to global in one call
- cell_dofs_global = dofmap.index_map.local_to_global(cell_dofs_local)
- cell_dofs = np.unique(cell_dofs_global)
+ length = (t1 - t0) * np.linalg.norm(beam_v)
+ weights = np.array([1.0, 4.0, 1.0]) * length / 6.0
+ local_dofs = self.V.dofmap.cell_dofs(cell_idx)
+ # global_dofs = self.V.dofmap.index_map.local_to_global(local_dofs)
+ # R[beam_idx, local_dofs] += weights @ tab
+ rows.append(np.full(len(local_dofs), beam_idx, dtype=np.int64))
+ cols.append(np.asarray(local_dofs, dtype=np.int64))
+ data.append(weights @ tab)
+
+ # return R
- weight = seg_len / len(cell_dofs)
- R[beam_idx, cell_dofs] += weight # ADD for overlaps!
+ if rows:
+ rows = np.concatenate(rows)
+ cols = np.concatenate(cols)
+ data = np.concatenate(data)
+ else:
+ rows = np.array([], dtype=np.int64)
+ cols = np.array([], dtype=np.int64)
+ data = np.array([], dtype=np.float64)
- return R
+ return coo_matrix((data, (rows, cols)), shape=(self.M, n_dofs)).tocsr()
def apply(self, x):
"""
@@ -143,15 +210,12 @@
# ⟨x.array,R₀^* v⟩_{ℝ^n} = ⟨Ax.array,[R^*v].array⟩_{ℝ^m}
# Taking [R^*v].array=A^{-1}R₀^* v, the conjugate works correctly.
tmp = fem.Function(self.V)
- np.matmul(self.R.T, j, out=tmp.x.array)
+ res = fem.Function(self.V)
+ tmp.x.array[:] = self.R.T @ j
tmp.x.scatter_forward()
- return tmp
- # TODO: is convection_diffusion actualy based on norms in ℝ^n?
- # res = fem.Function(self.V)
- # res.x.array[:] = 0.0
- # self.solver.solve(tmp.x.petsc_vec, res.x.petsc_vec)
- # res.x.scatter_forward()
- # return res
+ self.solver.solve(tmp.x.petsc_vec, res.x.petsc_vec)
+ res.x.scatter_forward()
+ return res
def opnorm(self, *args):
# raise NotImplementedError("Need mesh weighting?")