pointsource_pde: comparison src/laser

-:c4f923017349
+:a0d61ac3719e
+import basix
 import numpy as np
-from dolfinx import fem, geometry
+from dolfinx import fem
-from dolfinx_access import cell_FunctionSpace_f64
 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:
 """
 A linear operator class for laser sampling of a concentration (fem.Function)
 V,
 domain,
 nx,
 ny,
 beams,
-num_segments=100,
 domain_size=[-2, 2, -2, 2],
 ):
 self.V = V
 self.domain = domain
 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)
 A = get_mass_matrix(V)
 # SVD operator norm
 _u, s, _v = svds(self.R, k=1)
 self._opnorm = s[0]  # / sigma
 # print(f"Laser sampling opnorm = {self._opnorm:.6f}")
+# ------------------------------------------------------------------
+# 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.
+# ------------------------------------------------------------------
+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)]
+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)
+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
+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
 def matrix(self, beams):
-"""Combined: generate beams + segment + cell collision + dof weighting"""
+"""Combined: exact beam/cell intersection + exact per-cell line integration."""
-dofmap = self.V.dofmap
-tdim = self.domain.topology.dim
-# 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))
+# R = np.zeros((self.M, n_dofs))
+rows, cols, data = [], [], []
-xmin, xmax, ymin, ymax = self.domain_size
+for cell_idx in range(len(self.domain.geometry.dofmap)):
-for beam_idx, (start_pt, end_pt) in enumerate(beams):
+pts = self._cell_vertices_2d(cell_idx)
-# Segment into midpoints
+for beam_idx, (start_pt, end_pt) in enumerate(beams):
-vec = end_pt - start_pt
+clip = self._clip_segment_to_cell(pts, start_pt, end_pt)
-total_len = np.linalg.norm(vec)
+if clip is not None:
-seg_len = total_len / self.num_segments
+t0, t1 = clip
+if t1 - t0 > 1e-10:
-midpoints = []
+beam_v = end_pt - start_pt
-for i in range(self.num_segments):
+x0 = start_pt + t0 * beam_v
-t = (i + 0.5) / self.num_segments
+x1 = start_pt + t1 * beam_v
-midpoint = start_pt + t * vec
+midpt = start_pt + (0.5 * (t1 + t0)) * beam_v
-midpoints.append(midpoint)
+"""
-# midpoints = np.array(midpoints, dtype=np.float64)  # Force float64
+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]
-# Ensure shape (N, 3), C-contiguous, read-only for DOLFINx
+equals the integral. Uses the 3-point Simpson rule (endpoints + midpoint),
-# assert midpoints.shape[1] == 2, (
+which is exact since phi_h has degree <= 2 along the line.
-#     "Expected 2D points, got shape[1]=2 but need padding to 3D"
+"""
-# )
+x_geom = self.domain.geometry.x
-# midpoints = np.pad(
+verts = self.domain.geometry.dofmap[cell_idx]
-#     midpoints, ((0, 0), (0, 1)), mode="constant"
+v0, v1, v2 = (
-# )  # (N, 3) with z=0
+x_geom[verts[0], :2],
-# midpoints = np.ascontiguousarray(midpoints)  # C-order
+x_geom[verts[1], :2],
-# midpoints.setflags(write=False)  # Read-only requirement
+x_geom[verts[2], :2],
+)
-# Find colliding cells
+J = np.column_stack([v1 - v0, v2 - v0])
-# cell_candidates = geometry.compute_collisions_points(bb_tree, midpoints)
+"""Map a physical point p to reference coordinates on an affine triangle."""
-# colliding_cells = geometry.compute_colliding_cells(
+pts_ref = np.array(
-#    self.domain, cell_candidates, midpoints
+[
-# )
+np.linalg.solve(J, x0 - v0),
+np.linalg.solve(J, midpt - v0),
-# Cell→DOFs→weight accumulation
+np.linalg.solve(J, x1 - v0),
-# 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
+tab = self._element.tabulate(0, pts_ref)[
-for point in midpoints:
+0
-pad_point = np.array([point[0], point[1], 0.0])
+]  # (npoints, ndof_local, value_size)
-cell_idx = cell_FunctionSpace_f64(self.V._cpp_object, pad_point)
+tab = tab[
-if cell_idx >= 0:  # Valid cell index
+..., 0
-# Get global DOFs
+]  # scalar element -> drop trailing value_size=1 axis
-cell_dofs_local = dofmap.cell_dofs(cell_idx)
-# Convert local dofs to numpy int32 array (DOLFINx requirement)
+length = (t1 - t0) * np.linalg.norm(beam_v)
-# local_dofs_array = np.asarray(cell_dofs_local, dtype=np.int32)
+weights = np.array([1.0, 4.0, 1.0]) * length / 6.0
-# local_dofs_array.setflags(write=False)  # Read-only requirement
+local_dofs = self.V.dofmap.cell_dofs(cell_idx)
+# global_dofs = self.V.dofmap.index_map.local_to_global(local_dofs)
-# Batch convert all local indices to global in one call
+# R[beam_idx, local_dofs] += weights @ tab
-cell_dofs_global = dofmap.index_map.local_to_global(cell_dofs_local)
+rows.append(np.full(len(local_dofs), beam_idx, dtype=np.int64))
-cell_dofs = np.unique(cell_dofs_global)
+cols.append(np.asarray(local_dofs, dtype=np.int64))
+data.append(weights @ tab)
-weight = seg_len / len(cell_dofs)
-R[beam_idx, cell_dofs] += weight  # ADD for overlaps!
+# return R
-return R
+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 coo_matrix((data, (rows, cols)), shape=(self.M, n_dofs)).tocsr()
 def apply(self, x):
 """
 This does not work with MPI
 """
 # We have  ⟨x,R^*v⟩_V = ⟨Ax.array,[R^*v].array⟩_{ℝ^m}
 # But Rx = R₀ x.array, so we need
 # ⟨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
+self.solver.solve(tmp.x.petsc_vec, res.x.petsc_vec)
-# TODO: is convection_diffusion actualy based on norms in ℝ^n?
+res.x.scatter_forward()
-# res = fem.Function(self.V)
+return res
-# res.x.array[:] = 0.0
-# 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?")
 return self._opnorm

comparison: src/laser_sampling.py

src/laser_sampling.py

Mercurial > repos > pointsource_pde / file comparison

comparison: src/laser_sampling.py

src/laser_sampling.py