From 138fe2372890938bef8d91ee59a971dde2d9fe77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 27 Jul 2026 15:34:15 +0000 Subject: [PATCH 01/17] Add H1 boundary integral and unit tests for unit cube domain. --- .gitignore | 2 + src/struphy/feec/boundary_integrals.py | 224 ++++++++++++++++++ src/struphy/feec/mass_kernels.py | 55 +++++ .../feec/tests/test_boundary_integrals.py | 92 +++++++ 4 files changed, 373 insertions(+) create mode 100644 src/struphy/feec/boundary_integrals.py create mode 100644 src/struphy/feec/tests/test_boundary_integrals.py diff --git a/.gitignore b/.gitignore index c0e89a792..c602eab65 100644 --- a/.gitignore +++ b/.gitignore @@ -103,3 +103,5 @@ lib64 pyvenv.cfg *profile_output*.txt *kernels.txt + +examples/TwoFluidQuasiNeutralToy \ No newline at end of file diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py new file mode 100644 index 000000000..e88f4c43d --- /dev/null +++ b/src/struphy/feec/boundary_integrals.py @@ -0,0 +1,224 @@ +import cunumpy as xp +from typing import Callable + +from feectools.linalg.stencil import StencilVector +from feectools.linalg.block import BlockVector + +from struphy.feec import mass_kernels +from struphy.feec.mass import WeightedMassOperators +from struphy.feec.psydac_derham import Derham, SplineFunction +from struphy.geometry.base import Domain + +class BoundaryIntegralOperator: + """ + Assembles the boundary integral vector for H1 basis functions. + + Computes the six surface integrals + + I_i' = int_{partial Omega_i'} psi_h Tr(alpha) sqrt(g) |DF^-T n_hat_i| dS + + and adds them together into a single StencilVector v such that + + I = psi^T v + + for any discrete test function psi_h in V^0_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain = mass_ops.domain + + # H1 space info + self._space = self._derham.fem_spaces["0"] + self._space_key = "0" + + # 3D quadrature grid info for H1 space + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + + # for each of the 6 faces, extract surface quadrature grid and geometric weights + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + + # take quadrature points in the two surface directions + surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + + # build 2D meshgrid over surface + self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) + + # compute geometric weights + fixed_val = 0.0 if face_idx < 3 else 1.0 + + surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain.jacobian_det(*e_1d)) # metric + + DFinv = self._domain.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., :, normal_dir] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) # jacobian + + surface_geom_weights = sqrt_g * norm_DFinv_n + surface_geom_weights = xp.squeeze(surface_geom_weights) + self._surface_geom_weights.append(surface_geom_weights) + + + # extract surface spans, weights, bases for 2D quadrature + self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) # global index of the last nonzero spline + self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) # quadrature weights + self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) # spline values + + def _assemble_face( + self, + face_idx: int, + fun_weights: xp.ndarray, + dofs: StencilVector, + ): + """ + Assembles the contribution of a single face to the boundary integral vector. + + Parameters + ---------- + face_idx : int + Index of the face (0 to 5). + + fun_weights : xp.ndarray + Function alpha evaluated at the surface quadrature points, + already multiplied by the surface Jacobian. + + dofs : StencilVector + Output vector to accumulate into. + """ + + boundary_index = 0 if face_idx < 3 else -1 + + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + pads = fem_space.coeff_space.pads + + mass_kernels.surface_kernel_3d_vec( + *self._surface_spans[face_idx], + *fem_space.degree, + *starts, + *pads, + *self._surface_wts[face_idx], + *self._surface_bases[face_idx], + boundary_index, + fun_weights, + dofs._data, + ) + + + def assemble_callable( + self, + fun: Callable, + dofs: StencilVector = None, + clear: bool = True, + ) -> StencilVector: + """ + Assembles the boundary integral vector for a callable function alpha. + + Parameters + ---------- + fun : Callable + The function alpha(eta1, eta2, eta3) in logical coordinates. + + dofs : StencilVector, optional + Output vector. If None, a new zero vector is created. + + clear : bool, optional + Whether to zero the output vector before assembly. + + Returns + ------- + dofs : StencilVector + The assembled boundary integral vector v. + """ + if dofs is None: + dofs = self._space.coeff_space.zeros() + + if clear: + dofs._data[:] = 0.0 + + for face_idx in range(6): + normal_dir = face_idx % 3 + # fix the normal coordinate to 0.0 or 1.0 + fixed_val = 0.0 if face_idx < 3 else 1.0 + + surface_mesh = self._surface_quad_grid_meshes[face_idx] + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + + e = [None, None, None] + e[surf_dirs[0]] = surface_mesh[0] + e[surf_dirs[1]] = surface_mesh[1] + e[normal_dir] = xp.full_like(surface_mesh[0], fixed_val) + e1, e2, e3 = e + + fun_weights = fun(e1, e2, e3) + + # multiply by surface Jacobian + fun_weights = fun_weights * self._surface_geom_weights[face_idx] + fun_weights = xp.squeeze(fun_weights) + + self._assemble_face(face_idx, fun_weights, dofs) + + dofs.exchange_assembly_data() + dofs.update_ghost_regions() + + return dofs + + def __call__( + self, + fun: Callable | SplineFunction, + dofs: StencilVector = None, + clear: bool = True, + ) -> StencilVector: + """ + Assembles the boundary integral vector for a callable or SplineFunction alpha. + + Parameters + ---------- + fun : Callable | SplineFunction + The function alpha, either a callable or a SplineFunction. + + dofs : StencilVector, optional + Output vector. If None, a new zero vector is created. + + clear : bool, optional + Whether to zero the output vector before assembly. + + Returns + ------- + dofs : StencilVector + The assembled boundary integral vector v. + """ + if callable(fun): + return self.assemble_callable(fun, dofs=dofs, clear=clear) + else: + raise ValueError( + f"Expected callable, got {type(fun)} instead." + ) \ No newline at end of file diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 7b4f09720..f4e728b7d 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -766,3 +766,58 @@ def kernel_3d_diag( # No padding on StencilDiagonalMatrix data[i_local1, i_local2, i_local3] += value + + +def surface_kernel_3d_vec( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + boundary_index: int, + mat_fun: "float[:,:]", + data: "float[:,:,:]", +): + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + i_local0 = boundary_index - starts0 + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = ( + w1[iel1, q1] + * w2[iel2, q2] + * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + ) + + value += ( + wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] + ) + + data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value \ No newline at end of file diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py new file mode 100644 index 000000000..2d5094ae3 --- /dev/null +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -0,0 +1,92 @@ +import logging +from typing import Callable + +import cunumpy as xp +import pytest +from feectools.ddm.mpi import mpi as MPI + +from struphy import domains +from struphy.feec.boundary_integrals import BoundaryIntegralOperator +from struphy.feec.mass import WeightedMassOperators +from struphy.feec.psydac_derham import Derham +from struphy.io.options import DerhamOptions +from struphy.topology.grids import TensorProductGrid + +logger = logging.getLogger("struphy") + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_integral_callable(num_elements, degree, bcs): + """ + Tests the boundary integral operator for a callable function alpha on the + unit cube (Cuboid domain, identity mapping). + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: xp.ones_like(e1) + exact = 6.0 + + bnd_op = BoundaryIntegralOperator(mass_ops) + v = bnd_op.assemble_callable(alpha) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-10 + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_integral_callable_nonconstant(num_elements, degree, bcs): + """ + Tests the boundary integral operator for a non-constant callable alpha + on the unit cube (Cuboid domain, identity mapping). + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: e1 + e2 + e3 + exact = 9.0 + + bnd_op = BoundaryIntegralOperator(mass_ops) + v = bnd_op.assemble_callable(alpha) + + pads = v.space.pads + numerical = xp.sum(v._data[pads[0]:-pads[0], pads[1]:-pads[1], pads[2]:-pads[2]]) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-10 + + +if __name__ == "__main__": + from struphy import set_logging_level + set_logging_level(logging.INFO) + + test_boundary_integral_callable( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_integral_callable_nonconstant( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) \ No newline at end of file From 013cdb03b3a1705e72adccbf342b38d3fd92b6b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 27 Jul 2026 16:12:55 +0000 Subject: [PATCH 02/17] Add BC checks, general cuboid and hollow cylinder test cases. --- src/struphy/feec/boundary_integrals.py | 31 ++++++++- .../feec/tests/test_boundary_integrals.py | 69 +++++++++++++++++++ 2 files changed, 99 insertions(+), 1 deletion(-) diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py index e88f4c43d..b62ffc16e 100644 --- a/src/struphy/feec/boundary_integrals.py +++ b/src/struphy/feec/boundary_integrals.py @@ -55,7 +55,27 @@ def __init__( self._surface_wts = [] self._surface_bases = [] + self._active_faces = [] for face_idx in range(6): + normal_dir = face_idx % 3 + bc = self._derham.bcs[normal_dir] + + if bc is None: + self._active_faces.append(False) + elif face_idx < 3: + self._active_faces.append(bc[0] == "free") + else: + self._active_faces.append(bc[1] == "free") + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + normal_dir = face_idx % 3 surf_dirs = [d for d in range(3) if d != normal_dir] @@ -164,6 +184,9 @@ def assemble_callable( dofs._data[:] = 0.0 for face_idx in range(6): + if not self._active_faces[face_idx]: + continue + normal_dir = face_idx % 3 # fix the normal coordinate to 0.0 or 1.0 fixed_val = 0.0 if face_idx < 3 else 1.0 @@ -186,9 +209,15 @@ def assemble_callable( self._assemble_face(face_idx, fun_weights, dofs) + tmp = self._space.coeff_space.zeros() + self._assemble_face(face_idx, fun_weights, tmp) + tmp.exchange_assembly_data() + tmp.update_ghost_regions() + dofs.exchange_assembly_data() dofs.update_ghost_regions() + return dofs def __call__( @@ -221,4 +250,4 @@ def __call__( else: raise ValueError( f"Expected callable, got {type(fun)} instead." - ) \ No newline at end of file + ) diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 2d5094ae3..189431bd4 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -76,6 +76,65 @@ def test_boundary_integral_callable_nonconstant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-10 +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_integral_callable_cuboid_nontrivial(num_elements, degree, bcs): + """ + Tests the boundary integral operator for a non-constant callable alpha + on a non-cubic cuboid [0,1] x [0,2] x [0,3]. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=2.0, l3=0.0, r3=3.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: e1 + e2 + e3 + exact = 33.0 + + bnd_op = BoundaryIntegralOperator(mass_ops) + v = bnd_op.assemble_callable(alpha) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-10 + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): + """ + Tests the boundary integral operator for alpha = 1 on a HollowCylinder. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.HollowCylinder(a1=0.2, a2=1.0, Lz=4.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: xp.ones_like(e1) + exact = 11.52 * xp.pi + + bnd_op = BoundaryIntegralOperator(mass_ops) + v = bnd_op.assemble_callable(alpha) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + if __name__ == "__main__": from struphy import set_logging_level set_logging_level(logging.INFO) @@ -89,4 +148,14 @@ def test_boundary_integral_callable_nonconstant(num_elements, degree, bcs): [8, 8, 8], [2, 2, 2], (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_integral_callable_cuboid_nontrivial( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_integral_callable_hollow_cylinder( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), ) \ No newline at end of file From 0d85de3999c6ad62c60b8dc691915bf57ad5b33d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 28 Jul 2026 13:10:13 +0000 Subject: [PATCH 03/17] Fix bug in BoundaryIntegralOperator. --- src/struphy/feec/boundary_integrals.py | 3 +-- src/struphy/feec/tests/test_boundary_integrals.py | 10 +++++----- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py index b62ffc16e..ca879bdc5 100644 --- a/src/struphy/feec/boundary_integrals.py +++ b/src/struphy/feec/boundary_integrals.py @@ -98,7 +98,7 @@ def __init__( sqrt_g = xp.abs(self._domain.jacobian_det(*e_1d)) # metric DFinv = self._domain.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = DFinv[..., :, normal_dir] + DFinv_n = DFinv[..., normal_dir, :] norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) # jacobian surface_geom_weights = sqrt_g * norm_DFinv_n @@ -216,7 +216,6 @@ def assemble_callable( dofs.exchange_assembly_data() dofs.update_ghost_regions() - return dofs diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 189431bd4..928f87e03 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -82,7 +82,7 @@ def test_boundary_integral_callable_nonconstant(num_elements, degree, bcs): def test_boundary_integral_callable_cuboid_nontrivial(num_elements, degree, bcs): """ Tests the boundary integral operator for a non-constant callable alpha - on a non-cubic cuboid [0,1] x [0,2] x [0,3]. + on a non-unit cuboid [0,2]^3. """ comm = MPI.COMM_WORLD @@ -90,11 +90,11 @@ def test_boundary_integral_callable_cuboid_nontrivial(num_elements, degree, bcs) derham_opts = DerhamOptions(degree=degree, bcs=bcs) derham = Derham(grid, derham_opts, comm=comm) - domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=2.0, l3=0.0, r3=3.0) + domain = domains.Cuboid(l1=0.0, r1=2.0, l2=0.0, r2=2.0, l3=0.0, r3=2.0) mass_ops = WeightedMassOperators(derham, domain) alpha = lambda e1, e2, e3: e1 + e2 + e3 - exact = 33.0 + exact = 36.0 bnd_op = BoundaryIntegralOperator(mass_ops) v = bnd_op.assemble_callable(alpha) @@ -108,7 +108,7 @@ def test_boundary_integral_callable_cuboid_nontrivial(num_elements, degree, bcs) @pytest.mark.parametrize("num_elements", [[8, 8, 8]]) @pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): """ Tests the boundary integral operator for alpha = 1 on a HollowCylinder. @@ -157,5 +157,5 @@ def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): test_boundary_integral_callable_hollow_cylinder( [8, 8, 8], [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), + (("free", "free"), None, ("free", "free")), ) \ No newline at end of file From d932c175cb6a4417f2111a87083a16810ef25927 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 28 Jul 2026 15:33:37 +0000 Subject: [PATCH 04/17] Add H1 boundary mass matrix and unit tests. --- src/struphy/feec/boundary_integrals.py | 281 +++++++++++++++++- src/struphy/feec/mass_kernels.py | 81 ++++- .../feec/tests/test_boundary_integrals.py | 160 ++++++++++ 3 files changed, 518 insertions(+), 4 deletions(-) diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py index ca879bdc5..21c82e0b2 100644 --- a/src/struphy/feec/boundary_integrals.py +++ b/src/struphy/feec/boundary_integrals.py @@ -1,19 +1,23 @@ -import cunumpy as xp +import logging from typing import Callable -from feectools.linalg.stencil import StencilVector +import cunumpy as xp +from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL from feectools.linalg.block import BlockVector +from feectools.linalg.stencil import StencilMatrix, StencilVector from struphy.feec import mass_kernels +from struphy.feec.linear_operators import LinOpWithTransp from struphy.feec.mass import WeightedMassOperators from struphy.feec.psydac_derham import Derham, SplineFunction from struphy.geometry.base import Domain +from struphy.utils.pyccel import Pyccelkernel class BoundaryIntegralOperator: """ Assembles the boundary integral vector for H1 basis functions. - Computes the six surface integrals + Computes the surface integrals I_i' = int_{partial Omega_i'} psi_h Tr(alpha) sqrt(g) |DF^-T n_hat_i| dS @@ -250,3 +254,274 @@ def __call__( raise ValueError( f"Expected callable, got {type(fun)} instead." ) + + +class BoundaryMassOperator(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H1 basis functions. + + Computes the six surface integrals + + S_i'_{ijk,lmn} = int_{partial Omega_i'} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n_hat_i| dS + + and adds them together into a single StencilMatrix S such that + + I = psi^T S alpha + + for any discrete test function psi_h and spline function alpha_h in V^0_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + + # H1 space info + self._space = self._derham.fem_spaces["0"] + self._space_key = "0" + + # 3D quadrature grid info for H1 space + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + + # boundary and extraction operators + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + # initialize StencilMatrix + fem_space = self._tensor_fem_spaces[0] + self._mat = StencilMatrix( + fem_space.coeff_space, + fem_space.coeff_space, + backend=PSYDAC_BACKEND_GPYCCEL, + precompiled=True, + ) + + # build composite operator B * E * M * E^T * B^T + self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T + self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T + + # set domain and codomain + self._domain = self._M0.domain + self._codomain = self._M0.codomain + self._dtype = fem_space.coeff_space.dtype + + # allocate temporaries + self._temp_WB = self._W_boundary_op.domain.zeros() + self._temp_WE = self._W_extraction_op.domain.zeros() + self._temp_VB = self._V_boundary_op.domain.zeros() + self._temp_VE = self._V_extraction_op.domain.zeros() + self._temp_mat = self._mat.domain.zeros() + + # determine which faces to integrate over based on bcs + self._active_faces = [] + for face_idx in range(6): + normal_dir = face_idx % 3 + bc = self._derham.bcs[normal_dir] + + if bc is None: + self._active_faces.append(False) + elif face_idx < 3: + self._active_faces.append(bc[0] == "free") + else: + self._active_faces.append(bc[1] == "free") + + # for each active face, extract surface quadrature grid and geometric weights + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + fixed_val = 0.0 if face_idx < 3 else 1.0 + + surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) + + surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) + DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., normal_dir, :] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) + + surface_geom_weights = xp.squeeze(sqrt_g * norm_DFinv_n) + self._surface_geom_weights.append(surface_geom_weights) + + self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) + self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) + self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) + + # load assembly kernel + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat) + + self.assemble() + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face( + self, + face_idx: int, + mat: StencilMatrix, + ): + """ + Assembles the contribution of a single face to the boundary mass matrix. + + Parameters + ---------- + face_idx : int + Index of the face (0 to 5). + + mat : StencilMatrix + Output matrix to accumulate into. + """ + boundary_index = 0 if face_idx < 3 else -1 # TODO apparently not working + + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + pads = fem_space.coeff_space.pads + + self._assembly_kernel( + *self._surface_spans[face_idx], + *fem_space.degree, + *fem_space.degree, + *starts, + *pads, + *self._surface_wts[face_idx], + *self._surface_bases[face_idx], + *self._surface_bases[face_idx], + boundary_index, + self._surface_geom_weights[face_idx], + mat._data, + ) + + def assemble( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + if clear: + self._mat._data[:] = 0.0 + + for face_idx in range(6): + if not self._active_faces[face_idx]: + continue + self._assemble_face(face_idx, self._mat) + + self._mat.exchange_assembly_data() + self._mat.update_ghost_regions() + + def dot(self, v, out=None, apply_bc=True): + """ + Applies the boundary mass matrix to a vector. + + Parameters + ---------- + v : StencilVector + Input vector (spline coefficients of alpha_h). + + out : StencilVector, optional + Output vector. If None, a new zero vector is created. + + apply_bc : bool + Whether to apply boundary operators. + + Returns + ------- + out : StencilVector + The result S * v. + """ + if out is None: + out = self.codomain.zeros() + + if apply_bc: + self._V_boundary_op_T.dot(v, out=self._temp_VB) + self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) + self._W_boundary_op.dot(self._temp_WB, out=out) + else: + self._V_extraction_op_T.dot(v, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=out) + + return out + + def transpose(self, conjugate=False): + """ + Returns self since the boundary mass matrix is symmetric. + """ + return self + + def __call__( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + self.assemble(clear=clear) + return self + + + def toarray(self): + return self._M0.toarray() + + + def tosparse(self): + return self._M0.tosparse() \ No newline at end of file diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index f4e728b7d..a8e31fe49 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -820,4 +820,83 @@ def surface_kernel_3d_vec( wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] ) - data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value \ No newline at end of file + data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value + + +def surface_kernel_3d_mat( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + qi0: int, + qi1: int, + qi2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + boundary_index: int, + mat_fun: "float[:,:]", + data: "float[:,:,:,:,:,:]", +): + ne1 = spans1.size + ne2 = spans2.size + + nq1 = shape(w1)[1] + nq2 = shape(w2)[1] + + i_local0 = boundary_index - starts0 + + for iel1 in range(ne1): + for iel2 in range(ne2): + for il1 in range(pi1 + 1): + for il2 in range(pi2 + 1): + i_global1 = spans1[iel1] - pi1 + il1 + i_global2 = spans2[iel2] - pi2 + il2 + + i_local1 = i_global1 - starts1 + i_local2 = i_global2 - starts2 + + for jl1 in range(qi1 + 1): + for jl2 in range(qi2 + 1): + j_global1 = spans1[iel1] - qi1 + jl1 + j_global2 = spans2[iel2] - qi2 + jl2 + + j_local1 = j_global1 - i_global1 + j_local2 = j_global2 - i_global2 + + value = 0.0 + + for q1 in range(nq1): + for q2 in range(nq2): + wvol = ( + w1[iel1, q1] + * w2[iel2, q2] + * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] + ) + + value += ( + wvol + * bi1[iel1, il1, 0, q1] + * bi2[iel2, il2, 0, q2] + * bj1[iel1, jl1, 0, q1] + * bj2[iel2, jl2, 0, q2] + ) + + data[ + pads0 + i_local0, + pads1 + i_local1, + pads2 + i_local2, + pads0, + pads1 + j_local1, + pads2 + j_local2, + ] += value diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 928f87e03..e96a29dda 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -12,6 +12,9 @@ from struphy.io.options import DerhamOptions from struphy.topology.grids import TensorProductGrid +from struphy.feec.boundary_integrals import BoundaryIntegralOperator, BoundaryMassOperator +from struphy.feec.mass import L2Projector, WeightedMassOperators + logger = logging.getLogger("struphy") @@ -135,6 +138,142 @@ def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-3 +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = 1 on the unit cube. + S * alpha should give the same result as the callable version. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: xp.ones_like(e1) + exact = 6.0 + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_mass = BoundaryMassOperator(mass_ops) + v = bnd_mass.dot(alpha_h) + + logger.info(f"bnd_mass._mat._data max = {xp.max(xp.abs(bnd_mass._mat._data))}") + logger.info(f"alpha_h max = {xp.max(xp.abs(alpha_h.toarray()))}") + logger.info(f"v max = {xp.max(xp.abs(v.toarray()))}") + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 on the unit cube. + S * alpha should give the same result as the callable version. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: e1 + e2 + e3 + exact = 9.0 + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_mass = BoundaryMassOperator(mass_ops) + v = bnd_mass.dot(alpha_h) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 + on a non-unit cuboid [0,2]^3. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=2.0, l2=0.0, r2=2.0, l3=0.0, r3=2.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: e1 + e2 + e3 + exact = 36.0 + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_mass = BoundaryMassOperator(mass_ops) + v = bnd_mass.dot(alpha_h) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-3 + + +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) +def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = 1 on a HollowCylinder. + """ + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.HollowCylinder(a1=0.2, a2=1.0, Lz=4.0) + mass_ops = WeightedMassOperators(derham, domain) + + alpha = lambda e1, e2, e3: xp.ones_like(e1) + exact = 11.52 * xp.pi + + P = L2Projector("H1", mass_ops) + alpha_h = P(alpha) + + bnd_mass = BoundaryMassOperator(mass_ops) + v = bnd_mass.dot(alpha_h) + + numerical = xp.sum(v.toarray()) + + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-2 + + + if __name__ == "__main__": from struphy import set_logging_level set_logging_level(logging.INFO) @@ -158,4 +297,25 @@ def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): [8, 8, 8], [2, 2, 2], (("free", "free"), None, ("free", "free")), + ) + + test_boundary_mass_unit_cube_constant( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_mass_unit_cube_nonconstant( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_mass_cuboid_nontrivial( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + ) + test_boundary_mass_hollow_cylinder( + [8, 8, 8], + [2, 2, 2], + (("free", "free"), None, ("free", "free")), ) \ No newline at end of file From df5a25e5188d55e52105d7a7db6ba76fa4ba119b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Wed, 29 Jul 2026 11:32:38 +0000 Subject: [PATCH 05/17] Fix indexing in _assemble_face(). --- src/struphy/feec/boundary_integrals.py | 16 ++++++++++++++-- .../feec/tests/test_boundary_integrals.py | 4 ---- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py index 21c82e0b2..64c82a0c4 100644 --- a/src/struphy/feec/boundary_integrals.py +++ b/src/struphy/feec/boundary_integrals.py @@ -137,7 +137,13 @@ def _assemble_face( Output vector to accumulate into. """ - boundary_index = 0 if face_idx < 3 else -1 + normal_dir = face_idx % 3 + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + ends = [int(end) for end in fem_space.coeff_space.ends] + pads = fem_space.coeff_space.pads + + boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] fem_space = self._tensor_fem_spaces[0] starts = [int(start) for start in fem_space.coeff_space.starts] @@ -418,7 +424,13 @@ def _assemble_face( mat : StencilMatrix Output matrix to accumulate into. """ - boundary_index = 0 if face_idx < 3 else -1 # TODO apparently not working + normal_dir = face_idx % 3 + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + ends = [int(end) for end in fem_space.coeff_space.ends] + pads = fem_space.coeff_space.pads + + boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] fem_space = self._tensor_fem_spaces[0] starts = [int(start) for start in fem_space.coeff_space.starts] diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index e96a29dda..e6733c336 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -164,10 +164,6 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): bnd_mass = BoundaryMassOperator(mass_ops) v = bnd_mass.dot(alpha_h) - logger.info(f"bnd_mass._mat._data max = {xp.max(xp.abs(bnd_mass._mat._data))}") - logger.info(f"alpha_h max = {xp.max(xp.abs(alpha_h.toarray()))}") - logger.info(f"v max = {xp.max(xp.abs(v.toarray()))}") - numerical = xp.sum(v.toarray()) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") From 6b689425d2b0091b8270ecb1f8cec9d515667382 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Wed, 29 Jul 2026 13:40:39 +0000 Subject: [PATCH 06/17] Add factory class for boundary mass matrices and H(div), H(curl) drafts. --- src/struphy/feec/boundary_mass.py | 639 ++++++++++++++++++ .../feec/tests/test_boundary_integrals.py | 167 +---- 2 files changed, 649 insertions(+), 157 deletions(-) create mode 100644 src/struphy/feec/boundary_mass.py diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py new file mode 100644 index 000000000..4ab7d0a28 --- /dev/null +++ b/src/struphy/feec/boundary_mass.py @@ -0,0 +1,639 @@ +import logging +from typing import Callable + +import cunumpy as xp +from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL +from feectools.linalg.block import BlockVector +from feectools.linalg.stencil import StencilMatrix, StencilVector + +from struphy.feec import mass_kernels +from struphy.feec.linear_operators import LinOpWithTransp +from struphy.feec.mass import WeightedMassOperators +from struphy.feec.psydac_derham import Derham, SplineFunction +from struphy.geometry.base import Domain +from struphy.utils.pyccel import Pyccelkernel + + +class BoundaryOperators: + """ + Collection of boundary integral operators and boundary mass operators + for the H1, H(curl) and H(div) spaces. + + Analogous to WeightedMassOperators but for surface integrals. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool] | None = None, + ): + + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain = mass_ops.domain + + # shared surface setup for all spaces + # active faces based on bcs + if active_faces is not None: + # use provided active faces directly + self._active_faces = active_faces + else: + # default: integrate on free faces based on bcs + self._active_faces = [] + for face_idx in range(6): + normal_dir = face_idx % 3 + bc = self._derham.bcs[normal_dir] + if bc is None: + self._active_faces.append(False) + elif face_idx < 3: + self._active_faces.append(bc[0] == "free") + else: + self._active_faces.append(bc[1] == "free") + + # TODO: shared surface quad grids, geom weights, spans, wts, bases + # for each space (H1, Hcurl, Hdiv) — these differ because the + # quadrature grids are different for each space + + ################################################## + # H1 boundary operators (scalar, normal trace) # + ################################################## + + @property + def S0(self) -> "BoundaryMassOperatorH1": + """ + Boundary mass matrix for H1: + + S0_{ijk,lmn} = int_{partial Omega} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n| dS + """ + if not hasattr(self, "_S0"): + self._S0 = BoundaryMassOperatorH1(self._mass_ops, self._active_faces) + return self._S0 + + ################################################## + # H(curl) boundary operators (tangential trace) # + ################################################## + + @property + def S1(self) -> "BoundaryMassOperatorHCurl": + """ + Boundary mass matrix for H(curl): + + S1_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^1_{mu,ijk} x n) . Lambda^1_{nu,lmn} sqrt(g) |DF^-T n| dS + + Encodes the bilinear form for the tangential trace u x n against H(curl) test functions. + """ + if not hasattr(self, "_S1"): + self._S1 = BoundaryMassOperatorHCurl(self._mass_ops, self._active_faces) + return self._S1 + + ################################################## + # H(div) boundary operators (normal trace) # + ################################################## + + @property + def S2(self) -> "BoundaryMassOperatorHDiv": + """ + Boundary mass matrix for H(div): + + S2_{(mu,ijk),(nu,lmn)} = int_{partial Omega} Lambda^2_{mu,ijk} . n Lambda^2_{nu,lmn} . n sqrt(g) |DF^-T n| dS + + Encodes the bilinear form for the normal trace u . n against H(div) test functions. + """ + if not hasattr(self, "_S2"): + self._S2 = BoundaryMassOperatorHDiv(self._mass_ops, self._active_faces) + return self._S2 + + +class BoundaryMassOperatorH1(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H1 basis functions. + + Computes the six surface integrals + + S_i'_{ijk,lmn} = int_{partial Omega_i'} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n_hat_i| dS + + and adds them together into a single StencilMatrix S such that + + I = psi^T S alpha + + for any discrete test function psi_h and spline function alpha_h in V^0_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool] + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + + # H1 space info + self._space = self._derham.fem_spaces["0"] + self._space_key = "0" + + # 3D quadrature grid info for H1 space + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + + # boundary and extraction operators + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + # initialize StencilMatrix + fem_space = self._tensor_fem_spaces[0] + self._mat = StencilMatrix( + fem_space.coeff_space, + fem_space.coeff_space, + backend=PSYDAC_BACKEND_GPYCCEL, + precompiled=True, + ) + + # build composite operator B * E * M * E^T * B^T + self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T + self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T + + # set domain and codomain + self._domain = self._M0.domain + self._codomain = self._M0.codomain + self._dtype = fem_space.coeff_space.dtype + + # allocate temporaries + self._temp_WB = self._W_boundary_op.domain.zeros() + self._temp_WE = self._W_extraction_op.domain.zeros() + self._temp_VB = self._V_boundary_op.domain.zeros() + self._temp_VE = self._V_extraction_op.domain.zeros() + self._temp_mat = self._mat.domain.zeros() + + # determine which faces to integrate over based on bcs + self._active_faces = active_faces + + # for each active face, extract surface quadrature grid and geometric weights + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + fixed_val = 0.0 if face_idx < 3 else 1.0 + + surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) + + surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) + DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., normal_dir, :] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) + + surface_geom_weights = xp.squeeze(sqrt_g * norm_DFinv_n) + self._surface_geom_weights.append(surface_geom_weights) + + self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) + self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) + self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) + + # load assembly kernel + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat) + + self.assemble() + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face( + self, + face_idx: int, + mat: StencilMatrix, + ): + """ + Assembles the contribution of a single face to the boundary mass matrix. + + Parameters + ---------- + face_idx : int + Index of the face (0 to 5). + + mat : StencilMatrix + Output matrix to accumulate into. + """ + normal_dir = face_idx % 3 + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + ends = [int(end) for end in fem_space.coeff_space.ends] + pads = fem_space.coeff_space.pads + + boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] + + fem_space = self._tensor_fem_spaces[0] + starts = [int(start) for start in fem_space.coeff_space.starts] + pads = fem_space.coeff_space.pads + + self._assembly_kernel( + *self._surface_spans[face_idx], + *fem_space.degree, + *fem_space.degree, + *starts, + *pads, + *self._surface_wts[face_idx], + *self._surface_bases[face_idx], + *self._surface_bases[face_idx], + boundary_index, + self._surface_geom_weights[face_idx], + mat._data, + ) + + def assemble( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + if clear: + self._mat._data[:] = 0.0 + + for face_idx in range(6): + if not self._active_faces[face_idx]: + continue + self._assemble_face(face_idx, self._mat) + + self._mat.exchange_assembly_data() + self._mat.update_ghost_regions() + + def dot(self, v, out=None, apply_bc=True): + """ + Applies the boundary mass matrix to a vector. + + Parameters + ---------- + v : StencilVector + Input vector (spline coefficients of alpha_h). + + out : StencilVector, optional + Output vector. If None, a new zero vector is created. + + apply_bc : bool + Whether to apply boundary operators. + + Returns + ------- + out : StencilVector + The result S * v. + """ + if out is None: + out = self.codomain.zeros() + + if apply_bc: + self._V_boundary_op_T.dot(v, out=self._temp_VB) + self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) + self._W_boundary_op.dot(self._temp_WB, out=out) + else: + self._V_extraction_op_T.dot(v, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=out) + + return out + + def transpose(self, conjugate=False): + """ + Returns self since the boundary mass matrix is symmetric. + """ + return self + + def __call__( + self, + clear: bool = True, + ): + """ + Assembles the boundary mass matrix. + + Parameters + ---------- + clear : bool, optional + Whether to zero the matrix before assembly. + """ + self.assemble(clear=clear) + return self + + + def toarray(self): + return self._M0.toarray() + + + def tosparse(self): + return self._M0.tosparse() + + +class BoundaryMassOperatorHCurl(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H(curl) basis functions. + + Computes the surface integrals + + S1_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^1_{mu,ijk} x n) . Lambda^1_{nu,lmn} sqrt(g) |DF^-T n| dS + + such that I = u^T S1 alpha for any discrete H(curl) functions u_h and alpha_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + active_faces : list[bool] + Which of the six faces to integrate over. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool], + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + self._active_faces = active_faces + + self._space_key = "1" + self._space = self._derham.fem_spaces[self._space_key] + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + # TODO: initialize BlockLinearOperator (3x3 blocks) + self._mat = None + self._M = None + self._M0 = None + self._domain = None + self._codomain = None + self._dtype = None + + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + fixed_val = 0.0 if face_idx < 3 else 1.0 + + # TODO: use correct component quadrature points for H(curl) + surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) + DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., normal_dir, :] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) + self._surface_geom_weights.append(xp.squeeze(sqrt_g * norm_DFinv_n)) + + # TODO: surface meshes, spans, wts, bases per component + self._surface_quad_grid_meshes.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + + # TODO: load assembly kernel + self._assembly_kernel = None + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face(self, face_idx: int, mat): + # TODO + pass + + def assemble(self, clear: bool = True): + # TODO + pass + + def dot(self, v, out=None, apply_bc=True): + # TODO + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def toarray(self): + # TODO + raise NotImplementedError + + def tosparse(self): + # TODO + raise NotImplementedError + + +class BoundaryMassOperatorHDiv(LinOpWithTransp): + """ + Assembles the boundary mass matrix for H(div) basis functions. + + Computes the surface integrals + + S2_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^2_{mu,ijk} . n) (Lambda^2_{nu,lmn} . n) sqrt(g) |DF^-T n| dS + + such that I = u^T S2 alpha for any discrete H(div) functions u_h and alpha_h. + + Parameters + ---------- + mass_ops : WeightedMassOperators + Mass operators object, contains geometry and derham. + active_faces : list[bool] + Which of the six faces to integrate over. + """ + + def __init__( + self, + mass_ops: WeightedMassOperators, + active_faces: list[bool], + ): + self._mass_ops = mass_ops + self._derham = mass_ops.derham + self._domain_obj = mass_ops.domain + self._active_faces = active_faces + + self._space_key = "2" + self._space = self._derham.fem_spaces[self._space_key] + self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts + self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans + self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts + self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases + self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + + self._V_extraction_op = self._derham.extraction_ops[self._space_key] + self._W_extraction_op = self._derham.extraction_ops[self._space_key] + self._V_boundary_op = self._derham.boundary_ops[self._space_key] + self._W_boundary_op = self._derham.boundary_ops[self._space_key] + + self._V_extraction_op_T = self._V_extraction_op.T + self._W_extraction_op_T = self._W_extraction_op.T + self._V_boundary_op_T = self._V_boundary_op.T + self._W_boundary_op_T = self._W_boundary_op.T + + # TODO: initialize BlockLinearOperator (3x3 blocks, only diagonal nonzero) + self._mat = None + self._M = None + self._M0 = None + self._domain = None + self._codomain = None + self._dtype = None + + self._surface_quad_grid_meshes = [] + self._surface_geom_weights = [] + self._surface_spans = [] + self._surface_wts = [] + self._surface_bases = [] + + for face_idx in range(6): + if not self._active_faces[face_idx]: + self._surface_quad_grid_meshes.append(None) + self._surface_geom_weights.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + continue + + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] + fixed_val = 0.0 if face_idx < 3 else 1.0 + + # TODO: use correct component quadrature points for H(div) + surf_pts_1d = [self._quad_grid_pts[normal_dir][d].flatten() for d in surf_dirs] + e_1d = [None, None, None] + e_1d[surf_dirs[0]] = surf_pts_1d[0] + e_1d[surf_dirs[1]] = surf_pts_1d[1] + e_1d[normal_dir] = xp.array([fixed_val]) + + sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) + DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) + DFinv_n = DFinv[..., normal_dir, :] + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) + self._surface_geom_weights.append(xp.squeeze(sqrt_g * norm_DFinv_n)) + + # TODO: surface meshes, spans, wts, bases for normal_dir component only + self._surface_quad_grid_meshes.append(None) + self._surface_spans.append(None) + self._surface_wts.append(None) + self._surface_bases.append(None) + + # TODO: load assembly kernel + self._assembly_kernel = None + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._dtype + + def _assemble_face(self, face_idx: int, mat): + # TODO + pass + + def assemble(self, clear: bool = True): + # TODO + pass + + def dot(self, v, out=None, apply_bc=True): + # TODO + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def toarray(self): + # TODO + raise NotImplementedError + + def tosparse(self): + # TODO + raise NotImplementedError \ No newline at end of file diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index e6733c336..568d71901 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -6,145 +6,21 @@ from feectools.ddm.mpi import mpi as MPI from struphy import domains -from struphy.feec.boundary_integrals import BoundaryIntegralOperator -from struphy.feec.mass import WeightedMassOperators +from struphy.feec.boundary_mass import BoundaryOperators +from struphy.feec.mass import L2Projector, WeightedMassOperators from struphy.feec.psydac_derham import Derham from struphy.io.options import DerhamOptions from struphy.topology.grids import TensorProductGrid -from struphy.feec.boundary_integrals import BoundaryIntegralOperator, BoundaryMassOperator -from struphy.feec.mass import L2Projector, WeightedMassOperators - logger = logging.getLogger("struphy") -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -def test_boundary_integral_callable(num_elements, degree, bcs): - """ - Tests the boundary integral operator for a callable function alpha on the - unit cube (Cuboid domain, identity mapping). - """ - comm = MPI.COMM_WORLD - - grid = TensorProductGrid(num_elements=num_elements) - derham_opts = DerhamOptions(degree=degree, bcs=bcs) - derham = Derham(grid, derham_opts, comm=comm) - - domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) - mass_ops = WeightedMassOperators(derham, domain) - - alpha = lambda e1, e2, e3: xp.ones_like(e1) - exact = 6.0 - - bnd_op = BoundaryIntegralOperator(mass_ops) - v = bnd_op.assemble_callable(alpha) - - numerical = xp.sum(v.toarray()) - - logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - - assert xp.abs(numerical - exact) < 1e-10 - - -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -def test_boundary_integral_callable_nonconstant(num_elements, degree, bcs): - """ - Tests the boundary integral operator for a non-constant callable alpha - on the unit cube (Cuboid domain, identity mapping). - """ - comm = MPI.COMM_WORLD - - grid = TensorProductGrid(num_elements=num_elements) - derham_opts = DerhamOptions(degree=degree, bcs=bcs) - derham = Derham(grid, derham_opts, comm=comm) - - domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) - mass_ops = WeightedMassOperators(derham, domain) - - alpha = lambda e1, e2, e3: e1 + e2 + e3 - exact = 9.0 - - bnd_op = BoundaryIntegralOperator(mass_ops) - v = bnd_op.assemble_callable(alpha) - - pads = v.space.pads - numerical = xp.sum(v._data[pads[0]:-pads[0], pads[1]:-pads[1], pads[2]:-pads[2]]) - - logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - - assert xp.abs(numerical - exact) < 1e-10 - - -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -def test_boundary_integral_callable_cuboid_nontrivial(num_elements, degree, bcs): - """ - Tests the boundary integral operator for a non-constant callable alpha - on a non-unit cuboid [0,2]^3. - """ - comm = MPI.COMM_WORLD - - grid = TensorProductGrid(num_elements=num_elements) - derham_opts = DerhamOptions(degree=degree, bcs=bcs) - derham = Derham(grid, derham_opts, comm=comm) - - domain = domains.Cuboid(l1=0.0, r1=2.0, l2=0.0, r2=2.0, l3=0.0, r3=2.0) - mass_ops = WeightedMassOperators(derham, domain) - - alpha = lambda e1, e2, e3: e1 + e2 + e3 - exact = 36.0 - - bnd_op = BoundaryIntegralOperator(mass_ops) - v = bnd_op.assemble_callable(alpha) - - numerical = xp.sum(v.toarray()) - - logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - - assert xp.abs(numerical - exact) < 1e-10 - - -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) -def test_boundary_integral_callable_hollow_cylinder(num_elements, degree, bcs): - """ - Tests the boundary integral operator for alpha = 1 on a HollowCylinder. - """ - comm = MPI.COMM_WORLD - - grid = TensorProductGrid(num_elements=num_elements) - derham_opts = DerhamOptions(degree=degree, bcs=bcs) - derham = Derham(grid, derham_opts, comm=comm) - - domain = domains.HollowCylinder(a1=0.2, a2=1.0, Lz=4.0) - mass_ops = WeightedMassOperators(derham, domain) - - alpha = lambda e1, e2, e3: xp.ones_like(e1) - exact = 11.52 * xp.pi - - bnd_op = BoundaryIntegralOperator(mass_ops) - v = bnd_op.assemble_callable(alpha) - - numerical = xp.sum(v.toarray()) - - logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - - assert xp.abs(numerical - exact) < 1e-3 - - @pytest.mark.parametrize("num_elements", [[8, 8, 8]]) @pytest.mark.parametrize("degree", [[2, 2, 2]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = 1 on the unit cube. - S * alpha should give the same result as the callable version. """ comm = MPI.COMM_WORLD @@ -161,8 +37,8 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_mass = BoundaryMassOperator(mass_ops) - v = bnd_mass.dot(alpha_h) + bnd_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -177,7 +53,6 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 on the unit cube. - S * alpha should give the same result as the callable version. """ comm = MPI.COMM_WORLD @@ -194,8 +69,8 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_mass = BoundaryMassOperator(mass_ops) - v = bnd_mass.dot(alpha_h) + bnd_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -227,8 +102,8 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_mass = BoundaryMassOperator(mass_ops) - v = bnd_mass.dot(alpha_h) + bnd_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -259,8 +134,8 @@ def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_mass = BoundaryMassOperator(mass_ops) - v = bnd_mass.dot(alpha_h) + bnd_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -269,32 +144,10 @@ def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-2 - if __name__ == "__main__": from struphy import set_logging_level set_logging_level(logging.INFO) - test_boundary_integral_callable( - [8, 8, 8], - [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), - ) - test_boundary_integral_callable_nonconstant( - [8, 8, 8], - [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), - ) - test_boundary_integral_callable_cuboid_nontrivial( - [8, 8, 8], - [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), - ) - test_boundary_integral_callable_hollow_cylinder( - [8, 8, 8], - [2, 2, 2], - (("free", "free"), None, ("free", "free")), - ) - test_boundary_mass_unit_cube_constant( [8, 8, 8], [2, 2, 2], From 822b2079ea076f85f6ee20b255ac39ed837403ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Wed, 29 Jul 2026 14:06:27 +0000 Subject: [PATCH 07/17] Add draft surface kernels for H(div) and H(curl). --- src/struphy/feec/mass_kernels.py | 58 +++++++++++++++++++++++++++++++- 1 file changed, 57 insertions(+), 1 deletion(-) diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index a8e31fe49..796f6096d 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -823,7 +823,7 @@ def surface_kernel_3d_vec( data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value -def surface_kernel_3d_mat( +def surface_kernel_3d_mat_h1( spans1: "int[:]", spans2: "int[:]", pi0: int, @@ -900,3 +900,59 @@ def surface_kernel_3d_mat( pads1 + j_local1, pads2 + j_local2, ] += value + + +def surface_kernel_3d_mat_hdiv( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + qi0: int, + qi1: int, + qi2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + boundary_index: int, + mat_fun: "float[:,:]", + data: "float[:,:,:,:,:,:]", +): + pass + + +def surface_kernel_3d_mat_hcurl( + spans1: "int[:]", + spans2: "int[:]", + pi0: int, + pi1: int, + pi2: int, + qi0: int, + qi1: int, + qi2: int, + starts0: int, + starts1: int, + starts2: int, + pads0: int, + pads1: int, + pads2: int, + w1: "float[:,:]", + w2: "float[:,:]", + bi1: "float[:,:,:,:]", + bi2: "float[:,:,:,:]", + bj1: "float[:,:,:,:]", + bj2: "float[:,:,:,:]", + boundary_index: int, + n_cross_weight: "float[:,:]", + data: "float[:,:,:,:,:,:]", +): + pass \ No newline at end of file From 642295fecc1e7ee2085db5f5fc10ff954b206639 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 3 Aug 2026 11:08:48 +0000 Subject: [PATCH 08/17] Claude refactor for surface_kernel_3d_mat_h1. --- src/struphy/feec/boundary_integrals.py | 539 ------------------ src/struphy/feec/boundary_mass.py | 5 +- src/struphy/feec/mass_kernels.py | 109 ++-- .../feec/tests/test_boundary_integrals.py | 52 +- 4 files changed, 98 insertions(+), 607 deletions(-) delete mode 100644 src/struphy/feec/boundary_integrals.py diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py deleted file mode 100644 index 64c82a0c4..000000000 --- a/src/struphy/feec/boundary_integrals.py +++ /dev/null @@ -1,539 +0,0 @@ -import logging -from typing import Callable - -import cunumpy as xp -from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL -from feectools.linalg.block import BlockVector -from feectools.linalg.stencil import StencilMatrix, StencilVector - -from struphy.feec import mass_kernels -from struphy.feec.linear_operators import LinOpWithTransp -from struphy.feec.mass import WeightedMassOperators -from struphy.feec.psydac_derham import Derham, SplineFunction -from struphy.geometry.base import Domain -from struphy.utils.pyccel import Pyccelkernel - -class BoundaryIntegralOperator: - """ - Assembles the boundary integral vector for H1 basis functions. - - Computes the surface integrals - - I_i' = int_{partial Omega_i'} psi_h Tr(alpha) sqrt(g) |DF^-T n_hat_i| dS - - and adds them together into a single StencilVector v such that - - I = psi^T v - - for any discrete test function psi_h in V^0_h. - - Parameters - ---------- - mass_ops : WeightedMassOperators - Mass operators object, contains geometry and derham. - """ - - def __init__( - self, - mass_ops: WeightedMassOperators, - ): - self._mass_ops = mass_ops - self._derham = mass_ops.derham - self._domain = mass_ops.domain - - # H1 space info - self._space = self._derham.fem_spaces["0"] - self._space_key = "0" - - # 3D quadrature grid info for H1 space - self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts - self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans - self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts - self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases - self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces - - # for each of the 6 faces, extract surface quadrature grid and geometric weights - self._surface_quad_grid_meshes = [] - self._surface_geom_weights = [] - self._surface_spans = [] - self._surface_wts = [] - self._surface_bases = [] - - self._active_faces = [] - for face_idx in range(6): - normal_dir = face_idx % 3 - bc = self._derham.bcs[normal_dir] - - if bc is None: - self._active_faces.append(False) - elif face_idx < 3: - self._active_faces.append(bc[0] == "free") - else: - self._active_faces.append(bc[1] == "free") - - for face_idx in range(6): - if not self._active_faces[face_idx]: - self._surface_quad_grid_meshes.append(None) - self._surface_geom_weights.append(None) - self._surface_spans.append(None) - self._surface_wts.append(None) - self._surface_bases.append(None) - continue - - normal_dir = face_idx % 3 - surf_dirs = [d for d in range(3) if d != normal_dir] - - # take quadrature points in the two surface directions - surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] - - # build 2D meshgrid over surface - self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) - - # compute geometric weights - fixed_val = 0.0 if face_idx < 3 else 1.0 - - surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] - - e_1d = [None, None, None] - e_1d[surf_dirs[0]] = surf_pts_1d[0] - e_1d[surf_dirs[1]] = surf_pts_1d[1] - e_1d[normal_dir] = xp.array([fixed_val]) - - sqrt_g = xp.abs(self._domain.jacobian_det(*e_1d)) # metric - - DFinv = self._domain.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = DFinv[..., normal_dir, :] - norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) # jacobian - - surface_geom_weights = sqrt_g * norm_DFinv_n - surface_geom_weights = xp.squeeze(surface_geom_weights) - self._surface_geom_weights.append(surface_geom_weights) - - - # extract surface spans, weights, bases for 2D quadrature - self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) # global index of the last nonzero spline - self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) # quadrature weights - self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) # spline values - - def _assemble_face( - self, - face_idx: int, - fun_weights: xp.ndarray, - dofs: StencilVector, - ): - """ - Assembles the contribution of a single face to the boundary integral vector. - - Parameters - ---------- - face_idx : int - Index of the face (0 to 5). - - fun_weights : xp.ndarray - Function alpha evaluated at the surface quadrature points, - already multiplied by the surface Jacobian. - - dofs : StencilVector - Output vector to accumulate into. - """ - - normal_dir = face_idx % 3 - fem_space = self._tensor_fem_spaces[0] - starts = [int(start) for start in fem_space.coeff_space.starts] - ends = [int(end) for end in fem_space.coeff_space.ends] - pads = fem_space.coeff_space.pads - - boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] - - fem_space = self._tensor_fem_spaces[0] - starts = [int(start) for start in fem_space.coeff_space.starts] - pads = fem_space.coeff_space.pads - - mass_kernels.surface_kernel_3d_vec( - *self._surface_spans[face_idx], - *fem_space.degree, - *starts, - *pads, - *self._surface_wts[face_idx], - *self._surface_bases[face_idx], - boundary_index, - fun_weights, - dofs._data, - ) - - - def assemble_callable( - self, - fun: Callable, - dofs: StencilVector = None, - clear: bool = True, - ) -> StencilVector: - """ - Assembles the boundary integral vector for a callable function alpha. - - Parameters - ---------- - fun : Callable - The function alpha(eta1, eta2, eta3) in logical coordinates. - - dofs : StencilVector, optional - Output vector. If None, a new zero vector is created. - - clear : bool, optional - Whether to zero the output vector before assembly. - - Returns - ------- - dofs : StencilVector - The assembled boundary integral vector v. - """ - if dofs is None: - dofs = self._space.coeff_space.zeros() - - if clear: - dofs._data[:] = 0.0 - - for face_idx in range(6): - if not self._active_faces[face_idx]: - continue - - normal_dir = face_idx % 3 - # fix the normal coordinate to 0.0 or 1.0 - fixed_val = 0.0 if face_idx < 3 else 1.0 - - surface_mesh = self._surface_quad_grid_meshes[face_idx] - normal_dir = face_idx % 3 - surf_dirs = [d for d in range(3) if d != normal_dir] - - e = [None, None, None] - e[surf_dirs[0]] = surface_mesh[0] - e[surf_dirs[1]] = surface_mesh[1] - e[normal_dir] = xp.full_like(surface_mesh[0], fixed_val) - e1, e2, e3 = e - - fun_weights = fun(e1, e2, e3) - - # multiply by surface Jacobian - fun_weights = fun_weights * self._surface_geom_weights[face_idx] - fun_weights = xp.squeeze(fun_weights) - - self._assemble_face(face_idx, fun_weights, dofs) - - tmp = self._space.coeff_space.zeros() - self._assemble_face(face_idx, fun_weights, tmp) - tmp.exchange_assembly_data() - tmp.update_ghost_regions() - - dofs.exchange_assembly_data() - dofs.update_ghost_regions() - - return dofs - - def __call__( - self, - fun: Callable | SplineFunction, - dofs: StencilVector = None, - clear: bool = True, - ) -> StencilVector: - """ - Assembles the boundary integral vector for a callable or SplineFunction alpha. - - Parameters - ---------- - fun : Callable | SplineFunction - The function alpha, either a callable or a SplineFunction. - - dofs : StencilVector, optional - Output vector. If None, a new zero vector is created. - - clear : bool, optional - Whether to zero the output vector before assembly. - - Returns - ------- - dofs : StencilVector - The assembled boundary integral vector v. - """ - if callable(fun): - return self.assemble_callable(fun, dofs=dofs, clear=clear) - else: - raise ValueError( - f"Expected callable, got {type(fun)} instead." - ) - - -class BoundaryMassOperator(LinOpWithTransp): - """ - Assembles the boundary mass matrix for H1 basis functions. - - Computes the six surface integrals - - S_i'_{ijk,lmn} = int_{partial Omega_i'} Lambda^0_{ijk} Lambda^0_{lmn} sqrt(g) |DF^-T n_hat_i| dS - - and adds them together into a single StencilMatrix S such that - - I = psi^T S alpha - - for any discrete test function psi_h and spline function alpha_h in V^0_h. - - Parameters - ---------- - mass_ops : WeightedMassOperators - Mass operators object, contains geometry and derham. - """ - - def __init__( - self, - mass_ops: WeightedMassOperators, - ): - self._mass_ops = mass_ops - self._derham = mass_ops.derham - self._domain_obj = mass_ops.domain - - # H1 space info - self._space = self._derham.fem_spaces["0"] - self._space_key = "0" - - # 3D quadrature grid info for H1 space - self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts - self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans - self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts - self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases - self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces - - # boundary and extraction operators - self._V_extraction_op = self._derham.extraction_ops[self._space_key] - self._W_extraction_op = self._derham.extraction_ops[self._space_key] - self._V_boundary_op = self._derham.boundary_ops[self._space_key] - self._W_boundary_op = self._derham.boundary_ops[self._space_key] - - self._V_extraction_op_T = self._V_extraction_op.T - self._W_extraction_op_T = self._W_extraction_op.T - self._V_boundary_op_T = self._V_boundary_op.T - self._W_boundary_op_T = self._W_boundary_op.T - - # initialize StencilMatrix - fem_space = self._tensor_fem_spaces[0] - self._mat = StencilMatrix( - fem_space.coeff_space, - fem_space.coeff_space, - backend=PSYDAC_BACKEND_GPYCCEL, - precompiled=True, - ) - - # build composite operator B * E * M * E^T * B^T - self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T - self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T - - # set domain and codomain - self._domain = self._M0.domain - self._codomain = self._M0.codomain - self._dtype = fem_space.coeff_space.dtype - - # allocate temporaries - self._temp_WB = self._W_boundary_op.domain.zeros() - self._temp_WE = self._W_extraction_op.domain.zeros() - self._temp_VB = self._V_boundary_op.domain.zeros() - self._temp_VE = self._V_extraction_op.domain.zeros() - self._temp_mat = self._mat.domain.zeros() - - # determine which faces to integrate over based on bcs - self._active_faces = [] - for face_idx in range(6): - normal_dir = face_idx % 3 - bc = self._derham.bcs[normal_dir] - - if bc is None: - self._active_faces.append(False) - elif face_idx < 3: - self._active_faces.append(bc[0] == "free") - else: - self._active_faces.append(bc[1] == "free") - - # for each active face, extract surface quadrature grid and geometric weights - self._surface_quad_grid_meshes = [] - self._surface_geom_weights = [] - self._surface_spans = [] - self._surface_wts = [] - self._surface_bases = [] - - for face_idx in range(6): - if not self._active_faces[face_idx]: - self._surface_quad_grid_meshes.append(None) - self._surface_geom_weights.append(None) - self._surface_spans.append(None) - self._surface_wts.append(None) - self._surface_bases.append(None) - continue - - normal_dir = face_idx % 3 - surf_dirs = [d for d in range(3) if d != normal_dir] - fixed_val = 0.0 if face_idx < 3 else 1.0 - - surf_pts = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] - self._surface_quad_grid_meshes.append(xp.meshgrid(*surf_pts, indexing="ij")) - - surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] - e_1d = [None, None, None] - e_1d[surf_dirs[0]] = surf_pts_1d[0] - e_1d[surf_dirs[1]] = surf_pts_1d[1] - e_1d[normal_dir] = xp.array([fixed_val]) - - sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) - DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = DFinv[..., normal_dir, :] - norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) - - surface_geom_weights = xp.squeeze(sqrt_g * norm_DFinv_n) - self._surface_geom_weights.append(surface_geom_weights) - - self._surface_spans.append([self._spans_l[0][d] for d in surf_dirs]) - self._surface_wts.append([self._wts_l[0][d] for d in surf_dirs]) - self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) - - # load assembly kernel - self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat) - - self.assemble() - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._codomain - - @property - def dtype(self): - return self._dtype - - def _assemble_face( - self, - face_idx: int, - mat: StencilMatrix, - ): - """ - Assembles the contribution of a single face to the boundary mass matrix. - - Parameters - ---------- - face_idx : int - Index of the face (0 to 5). - - mat : StencilMatrix - Output matrix to accumulate into. - """ - normal_dir = face_idx % 3 - fem_space = self._tensor_fem_spaces[0] - starts = [int(start) for start in fem_space.coeff_space.starts] - ends = [int(end) for end in fem_space.coeff_space.ends] - pads = fem_space.coeff_space.pads - - boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] - - fem_space = self._tensor_fem_spaces[0] - starts = [int(start) for start in fem_space.coeff_space.starts] - pads = fem_space.coeff_space.pads - - self._assembly_kernel( - *self._surface_spans[face_idx], - *fem_space.degree, - *fem_space.degree, - *starts, - *pads, - *self._surface_wts[face_idx], - *self._surface_bases[face_idx], - *self._surface_bases[face_idx], - boundary_index, - self._surface_geom_weights[face_idx], - mat._data, - ) - - def assemble( - self, - clear: bool = True, - ): - """ - Assembles the boundary mass matrix. - - Parameters - ---------- - clear : bool, optional - Whether to zero the matrix before assembly. - """ - if clear: - self._mat._data[:] = 0.0 - - for face_idx in range(6): - if not self._active_faces[face_idx]: - continue - self._assemble_face(face_idx, self._mat) - - self._mat.exchange_assembly_data() - self._mat.update_ghost_regions() - - def dot(self, v, out=None, apply_bc=True): - """ - Applies the boundary mass matrix to a vector. - - Parameters - ---------- - v : StencilVector - Input vector (spline coefficients of alpha_h). - - out : StencilVector, optional - Output vector. If None, a new zero vector is created. - - apply_bc : bool - Whether to apply boundary operators. - - Returns - ------- - out : StencilVector - The result S * v. - """ - if out is None: - out = self.codomain.zeros() - - if apply_bc: - self._V_boundary_op_T.dot(v, out=self._temp_VB) - self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) - self._mat.dot(self._temp_mat, out=self._temp_WE) - self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) - self._W_boundary_op.dot(self._temp_WB, out=out) - else: - self._V_extraction_op_T.dot(v, out=self._temp_mat) - self._mat.dot(self._temp_mat, out=self._temp_WE) - self._W_extraction_op.dot(self._temp_WE, out=out) - - return out - - def transpose(self, conjugate=False): - """ - Returns self since the boundary mass matrix is symmetric. - """ - return self - - def __call__( - self, - clear: bool = True, - ): - """ - Assembles the boundary mass matrix. - - Parameters - ---------- - clear : bool, optional - Whether to zero the matrix before assembly. - """ - self.assemble(clear=clear) - return self - - - def toarray(self): - return self._M0.toarray() - - - def tosparse(self): - return self._M0.tosparse() \ No newline at end of file diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 4ab7d0a28..8e01e55b7 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -14,7 +14,7 @@ from struphy.utils.pyccel import Pyccelkernel -class BoundaryOperators: +class BoundaryIntegralOperators: """ Collection of boundary integral operators and boundary mass operators for the H1, H(curl) and H(div) spaces. @@ -230,7 +230,7 @@ def __init__( self._surface_bases.append([self._bases_l[0][d] for d in surf_dirs]) # load assembly kernel - self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat) + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat_h1) self.assemble() @@ -284,6 +284,7 @@ def _assemble_face( *self._surface_bases[face_idx], *self._surface_bases[face_idx], boundary_index, + normal_dir, self._surface_geom_weights[face_idx], mat._data, ) diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 796f6096d..7d54d5574 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -21,7 +21,7 @@ def kernel_1d_mat( data: "float[:,:]", ): """ - Performs the integration of Lambda_i * mat_fun(eta1) * Lambda_l for the basis functions (i, l) available on the calling process. + Performs the integration of Lambda_i * mat_fun(eta1) * Lambda_j for the basis functions (i, j) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). """ @@ -845,6 +845,7 @@ def surface_kernel_3d_mat_h1( bj1: "float[:,:,:,:]", bj2: "float[:,:,:,:]", boundary_index: int, + normal_dir: int, mat_fun: "float[:,:]", data: "float[:,:,:,:,:,:]", ): @@ -854,22 +855,43 @@ def surface_kernel_3d_mat_h1( nq1 = shape(w1)[1] nq2 = shape(w2)[1] - i_local0 = boundary_index - starts0 + starts = [starts0, starts1, starts2] + pads = [pads0, pads1, pads2] + pi = [pi0, pi1, pi2] + qi = [qi0, qi1, qi2] + + surf_dirs = [d for d in range(3) if d != normal_dir] + + pi_s1 = pi[surf_dirs[0]] + pi_s2 = pi[surf_dirs[1]] + + qi_s1 = qi[surf_dirs[0]] + qi_s2 = qi[surf_dirs[1]] + + starts_n = starts[normal_dir] + starts_s1 = starts[surf_dirs[0]] + starts_s2 = starts[surf_dirs[1]] + + pads_n = pads[normal_dir] + pads_s1 = pads[surf_dirs[0]] + pads_s2 = pads[surf_dirs[1]] + + i_local_n = boundary_index - starts_n for iel1 in range(ne1): for iel2 in range(ne2): - for il1 in range(pi1 + 1): - for il2 in range(pi2 + 1): - i_global1 = spans1[iel1] - pi1 + il1 - i_global2 = spans2[iel2] - pi2 + il2 + for il1 in range(pi_s1 + 1): + for il2 in range(pi_s2 + 1): + i_global1 = spans1[iel1] - pi_s1 + il1 + i_global2 = spans2[iel2] - pi_s2 + il2 - i_local1 = i_global1 - starts1 - i_local2 = i_global2 - starts2 + i_local1 = i_global1 - starts_s1 + i_local2 = i_global2 - starts_s2 - for jl1 in range(qi1 + 1): - for jl2 in range(qi2 + 1): - j_global1 = spans1[iel1] - qi1 + jl1 - j_global2 = spans2[iel2] - qi2 + jl2 + for jl1 in range(qi_s1 + 1): + for jl2 in range(qi_s2 + 1): + j_global1 = spans1[iel1] - qi_s1 + jl1 + j_global2 = spans2[iel2] - qi_s2 + jl2 j_local1 = j_global1 - i_global1 j_local2 = j_global2 - i_global2 @@ -892,42 +914,33 @@ def surface_kernel_3d_mat_h1( * bj2[iel2, jl2, 0, q2] ) - data[ - pads0 + i_local0, - pads1 + i_local1, - pads2 + i_local2, - pads0, - pads1 + j_local1, - pads2 + j_local2, - ] += value - - -def surface_kernel_3d_mat_hdiv( - spans1: "int[:]", - spans2: "int[:]", - pi0: int, - pi1: int, - pi2: int, - qi0: int, - qi1: int, - qi2: int, - starts0: int, - starts1: int, - starts2: int, - pads0: int, - pads1: int, - pads2: int, - w1: "float[:,:]", - w2: "float[:,:]", - bi1: "float[:,:,:,:]", - bi2: "float[:,:,:,:]", - bj1: "float[:,:,:,:]", - bj2: "float[:,:,:,:]", - boundary_index: int, - mat_fun: "float[:,:]", - data: "float[:,:,:,:,:,:]", -): - pass + if normal_dir == 0: + data[ + pads_n + i_local_n, + pads_s1 + i_local1, + pads_s2 + i_local2, + pads_n, + pads_s1 + j_local1, + pads_s2 + j_local2, + ] += value + elif normal_dir == 1: + data[ + pads_s1 + i_local1, + pads_n + i_local_n, + pads_s2 + i_local2, + pads_s1 + j_local1, + pads_n, + pads_s2 + j_local2, + ] += value + else: + data[ + pads_s1 + i_local1, + pads_s2 + i_local2, + pads_n + i_local_n, + pads_s1 + j_local1, + pads_s2 + j_local2, + pads_n, + ] += value def surface_kernel_3d_mat_hcurl( diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 568d71901..e7b09b498 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -6,7 +6,7 @@ from feectools.ddm.mpi import mpi as MPI from struphy import domains -from struphy.feec.boundary_mass import BoundaryOperators +from struphy.feec.boundary_mass import BoundaryIntegralOperators from struphy.feec.mass import L2Projector, WeightedMassOperators from struphy.feec.psydac_derham import Derham from struphy.io.options import DerhamOptions @@ -15,9 +15,14 @@ logger = logging.getLogger("struphy") -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize("num_elements", [[8, 9, 10]],) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize("bcs", [(None, None, None), + (("free", "free"), None, None), + (None, ("free", "free"), None), + (None, None, ("free", "free")), + (("free", "free"), ("free", "free"), None), + (("free", "free"), ("free", "free"), ("free", "free")),]) def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = 1 on the unit cube. @@ -31,13 +36,24 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) mass_ops = WeightedMassOperators(derham, domain) - alpha = lambda e1, e2, e3: xp.ones_like(e1) - exact = 6.0 + face_value = 2.1 + alpha = lambda e1, e2, e3: xp.ones_like(e1) * face_value + + num_faces = 0 + for face_tuple in bcs: + if face_tuple is None: + continue + if face_tuple[0] == "free": + num_faces += 1 + if face_tuple[1] == "free": + num_faces += 1 + + exact = num_faces * face_value P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_ops = BoundaryOperators(mass_ops) + bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -47,8 +63,8 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-3 -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): """ @@ -69,7 +85,7 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_ops = BoundaryOperators(mass_ops) + bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -79,8 +95,8 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-3 -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): """ @@ -102,7 +118,7 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_ops = BoundaryOperators(mass_ops) + bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -112,8 +128,8 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-3 -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2], [3, 3, 3]]) +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): """ @@ -134,7 +150,7 @@ def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): P = L2Projector("H1", mass_ops) alpha_h = P(alpha) - bnd_ops = BoundaryOperators(mass_ops) + bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) numerical = xp.sum(v.toarray()) @@ -149,8 +165,8 @@ def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): set_logging_level(logging.INFO) test_boundary_mass_unit_cube_constant( - [8, 8, 8], - [2, 2, 2], + [8, 9, 10], + [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), ) test_boundary_mass_unit_cube_nonconstant( From e3fad1ffe6d6347e800ada566f298f79eb0d1c4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 3 Aug 2026 12:24:57 +0000 Subject: [PATCH 09/17] Change to more complex unit tests. --- .../feec/tests/test_boundary_integrals.py | 64 +++++++++++++------ 1 file changed, 43 insertions(+), 21 deletions(-) diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index e7b09b498..3e8660ff5 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -65,7 +65,9 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): @pytest.mark.parametrize("num_elements", [[8, 9, 10]]) @pytest.mark.parametrize("degree", [[1, 2, 3]]) -@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize("bcs", [(("dirichlet", "free"), ("free", "free"), ("free", "free")), + (("free", "dirichlet"), ("free", "free"), ("free", "free")), + (("dirichlet", "dirichlet"), ("free", "free"), ("free", "free"))]) def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 on the unit cube. @@ -79,11 +81,19 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) mass_ops = WeightedMassOperators(derham, domain) - alpha = lambda e1, e2, e3: e1 + e2 + e3 - exact = 9.0 + if bcs[0] == ("dirichlet", "free"): + alpha = lambda e1, e2, e3: e1 + 0 * e2 + 0 * e3 + exact = 3.0 + elif bcs[0] == ("free", "dirichlet"): + alpha = lambda e1, e2, e3: 1.0 - e1 + 0 * e2 + 0 * e3 + exact = 3.0 + else: + assert bcs[0] == ("dirichlet", "dirichlet") + alpha = lambda e1, e2, e3: e1 * (1.0 - e1) + 0 * e2 + 0 * e3 + exact = 2.0 / 3.0 P = L2Projector("H1", mass_ops) - alpha_h = P(alpha) + alpha_h = P(alpha, apply_bc=True) bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) @@ -92,7 +102,7 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - assert xp.abs(numerical - exact) < 1e-3 + assert xp.abs(numerical - exact) < 2e-2 @pytest.mark.parametrize("num_elements", [[8, 9, 10]]) @@ -101,7 +111,7 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 - on a non-unit cuboid [0,2]^3. + on a non-unit cuboid [-1,1] x [-1,3] x [0,3]. """ comm = MPI.COMM_WORLD @@ -109,11 +119,11 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): derham_opts = DerhamOptions(degree=degree, bcs=bcs) derham = Derham(grid, derham_opts, comm=comm) - domain = domains.Cuboid(l1=0.0, r1=2.0, l2=0.0, r2=2.0, l3=0.0, r3=2.0) + domain = domains.Cuboid(l1=-1.0, r1=1.0, l2=-1.0, r2=3.0, l3=0.0, r3=3.0) mass_ops = WeightedMassOperators(derham, domain) alpha = lambda e1, e2, e3: e1 + e2 + e3 - exact = 36.0 + exact = 78.0 P = L2Projector("H1", mass_ops) alpha_h = P(alpha) @@ -131,21 +141,31 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): @pytest.mark.parametrize("num_elements", [[8, 9, 10]]) @pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), None, ("free", "free"))]) -def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): +def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): """ - Tests the boundary mass operator for alpha = 1 on a HollowCylinder. + Tests the boundary mass operator for alpha = exp(eta3) on a HollowCylinder. """ + import math comm = MPI.COMM_WORLD grid = TensorProductGrid(num_elements=num_elements) derham_opts = DerhamOptions(degree=degree, bcs=bcs) derham = Derham(grid, derham_opts, comm=comm) - domain = domains.HollowCylinder(a1=0.2, a2=1.0, Lz=4.0) + a1 = 0.2 + a2 = 1.0 + Lz = 4.0 + + domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) mass_ops = WeightedMassOperators(derham, domain) - alpha = lambda e1, e2, e3: xp.ones_like(e1) - exact = 11.52 * xp.pi + alpha = lambda e1, e2, e3: xp.exp(e3) + e = math.e + exact = xp.pi * ( + 2 * a1 * Lz * (e - 1) + + 2 * a2 * Lz * (e - 1) + + (a2**2 - a1**2) * (1 + e) + ) P = L2Projector("H1", mass_ops) alpha_h = P(alpha) @@ -169,18 +189,20 @@ def test_boundary_mass_hollow_cylinder(num_elements, degree, bcs): [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), ) + test_boundary_mass_unit_cube_nonconstant( - [8, 8, 8], - [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), + [8, 9, 10], + [1, 2, 3], + (("dirichlet", "free"), ("free", "free"), ("free", "free")), ) + test_boundary_mass_cuboid_nontrivial( - [8, 8, 8], - [2, 2, 2], + [8, 9, 10], + [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), ) - test_boundary_mass_hollow_cylinder( - [8, 8, 8], - [2, 2, 2], + test_boundary_mass_hollow_cylinder_nonconstant( + [8, 9, 10], + [1, 2, 3], (("free", "free"), None, ("free", "free")), ) \ No newline at end of file From 7943d0e12ceb407e8546a2950df5b8a31bb762bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 3 Aug 2026 15:50:35 +0000 Subject: [PATCH 10/17] Add initial H(curl) boundary mass operator. --- src/struphy/feec/boundary_mass.py | 348 +++++++++--------- .../feec/tests/test_boundary_integrals.py | 46 ++- 2 files changed, 216 insertions(+), 178 deletions(-) diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 8e01e55b7..9eb83d11c 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -13,6 +13,8 @@ from struphy.geometry.base import Domain from struphy.utils.pyccel import Pyccelkernel +from feectools.linalg.block import BlockLinearOperator, BlockVector + class BoundaryIntegralOperators: """ @@ -91,23 +93,6 @@ def S1(self) -> "BoundaryMassOperatorHCurl": self._S1 = BoundaryMassOperatorHCurl(self._mass_ops, self._active_faces) return self._S1 - ################################################## - # H(div) boundary operators (normal trace) # - ################################################## - - @property - def S2(self) -> "BoundaryMassOperatorHDiv": - """ - Boundary mass matrix for H(div): - - S2_{(mu,ijk),(nu,lmn)} = int_{partial Omega} Lambda^2_{mu,ijk} . n Lambda^2_{nu,lmn} . n sqrt(g) |DF^-T n| dS - - Encodes the bilinear form for the normal trace u . n against H(div) test functions. - """ - if not hasattr(self, "_S2"): - self._S2 = BoundaryMassOperatorHDiv(self._mass_ops, self._active_faces) - return self._S2 - class BoundaryMassOperatorH1(LinOpWithTransp): """ @@ -384,9 +369,13 @@ class BoundaryMassOperatorHCurl(LinOpWithTransp): Computes the surface integrals - S1_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^1_{mu,ijk} x n) . Lambda^1_{nu,lmn} sqrt(g) |DF^-T n| dS + W^{mu,nu}_{ijk,lmn} = int_{partial Omega} hat_Lambda^1_{mu,ijk} hat_R_n^{mu,nu} hat_Lambda^1_{nu,lmn} dS - such that I = u^T S1 alpha for any discrete H(curl) functions u_h and alpha_h. + where hat_R_n is the pullback of [n]_x to logical coordinates. + + The result is a 3x3 BlockLinearOperator where diagonal blocks are zero + (skew-symmetry of [n]_x) and off-diagonal blocks are assembled via + surface_kernel_3d_mat_h1. Parameters ---------- @@ -424,24 +413,53 @@ def __init__( self._V_boundary_op_T = self._V_boundary_op.T self._W_boundary_op_T = self._W_boundary_op.T - # TODO: initialize BlockLinearOperator (3x3 blocks) - self._mat = None - self._M = None - self._M0 = None - self._domain = None - self._codomain = None - self._dtype = None + V = self._space + W = self._space + + blocks = [ + [ + StencilMatrix( + Vs.coeff_space, + Ws.coeff_space, + backend=PSYDAC_BACKEND_GPYCCEL, + precompiled=True, + ) + if i != j + else None + for j, Vs in enumerate(V.spaces) + ] + for i, Ws in enumerate(W.spaces) + ] + + self._mat = BlockLinearOperator( + V.coeff_space, + W.coeff_space, + blocks=blocks, + ) - self._surface_quad_grid_meshes = [] - self._surface_geom_weights = [] + self._M = self._W_extraction_op @ self._mat @ self._V_extraction_op_T + self._M0 = self._W_boundary_op @ self._M @ self._V_boundary_op_T + + self._domain = self._M0.domain + self._codomain = self._M0.codomain + self._dtype = self._tensor_fem_spaces[0].coeff_space.dtype + + # allocate temporaries + self._temp_VB = self._V_boundary_op.domain.zeros() + self._temp_VE = self._V_extraction_op.domain.zeros() + self._temp_WB = self._W_boundary_op.domain.zeros() + self._temp_WE = self._W_extraction_op.domain.zeros() + self._temp_mat = self._mat.domain.zeros() + + # for each active face extract surface data + self._surface_R_n = [] self._surface_spans = [] self._surface_wts = [] self._surface_bases = [] for face_idx in range(6): if not self._active_faces[face_idx]: - self._surface_quad_grid_meshes.append(None) - self._surface_geom_weights.append(None) + self._surface_R_n.append(None) self._surface_spans.append(None) self._surface_wts.append(None) self._surface_bases.append(None) @@ -451,27 +469,50 @@ def __init__( surf_dirs = [d for d in range(3) if d != normal_dir] fixed_val = 0.0 if face_idx < 3 else 1.0 - # TODO: use correct component quadrature points for H(curl) - surf_pts_1d = [self._quad_grid_pts[0][d].flatten() for d in surf_dirs] + # compute R_n using surf_dirs[0] quadrature points + surf_pts_1d = [self._quad_grid_pts[surf_dirs[0]][d].flatten() for d in surf_dirs] e_1d = [None, None, None] e_1d[surf_dirs[0]] = surf_pts_1d[0] e_1d[surf_dirs[1]] = surf_pts_1d[1] e_1d[normal_dir] = xp.array([fixed_val]) - sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = DFinv[..., normal_dir, :] - norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) - self._surface_geom_weights.append(xp.squeeze(sqrt_g * norm_DFinv_n)) + DFinv_n = xp.squeeze(DFinv[..., normal_dir, :]) + norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1, keepdims=True)) + n_phys = DFinv_n / norm_DFinv_n + if face_idx >= 3: + n_phys = -n_phys + + n0 = n_phys[..., 0] + n1 = n_phys[..., 1] + n2 = n_phys[..., 2] + + R_n = xp.zeros((*n0.shape, 3, 3)) + R_n[..., 0, 1] = -n2 + R_n[..., 0, 2] = n1 + R_n[..., 1, 0] = n2 + R_n[..., 1, 2] = -n0 + R_n[..., 2, 0] = -n1 + R_n[..., 2, 1] = n0 + self._surface_R_n.append(R_n) + + # only extract for tangential components mu != normal_dir + surface_spans_per_mu = [None, None, None] + surface_wts_per_mu = [None, None, None] + surface_bases_per_mu = [None, None, None] + + for mu in surf_dirs: + surface_spans_per_mu[mu] = [self._spans_l[mu][d] for d in surf_dirs] + surface_wts_per_mu[mu] = [self._wts_l[mu][d] for d in surf_dirs] + surface_bases_per_mu[mu] = [self._bases_l[mu][d] for d in surf_dirs] + + self._surface_spans.append(surface_spans_per_mu) + self._surface_wts.append(surface_wts_per_mu) + self._surface_bases.append(surface_bases_per_mu) - # TODO: surface meshes, spans, wts, bases per component - self._surface_quad_grid_meshes.append(None) - self._surface_spans.append(None) - self._surface_wts.append(None) - self._surface_bases.append(None) + self._assembly_kernel = Pyccelkernel(mass_kernels.surface_kernel_3d_mat_h1) - # TODO: load assembly kernel - self._assembly_kernel = None + self.assemble() @property def domain(self): @@ -486,155 +527,110 @@ def dtype(self): return self._dtype def _assemble_face(self, face_idx: int, mat): - # TODO - pass - - def assemble(self, clear: bool = True): - # TODO - pass - - def dot(self, v, out=None, apply_bc=True): - # TODO - raise NotImplementedError - - def transpose(self, conjugate=False): - return self - - def toarray(self): - # TODO - raise NotImplementedError - - def tosparse(self): - # TODO - raise NotImplementedError - - -class BoundaryMassOperatorHDiv(LinOpWithTransp): - """ - Assembles the boundary mass matrix for H(div) basis functions. - - Computes the surface integrals - - S2_{(mu,ijk),(nu,lmn)} = int_{partial Omega} (Lambda^2_{mu,ijk} . n) (Lambda^2_{nu,lmn} . n) sqrt(g) |DF^-T n| dS - - such that I = u^T S2 alpha for any discrete H(div) functions u_h and alpha_h. - - Parameters - ---------- - mass_ops : WeightedMassOperators - Mass operators object, contains geometry and derham. - active_faces : list[bool] - Which of the six faces to integrate over. - """ - - def __init__( - self, - mass_ops: WeightedMassOperators, - active_faces: list[bool], - ): - self._mass_ops = mass_ops - self._derham = mass_ops.derham - self._domain_obj = mass_ops.domain - self._active_faces = active_faces - - self._space_key = "2" - self._space = self._derham.fem_spaces[self._space_key] - self._quad_grid_pts = self._derham.spline_attributes[self._space_key].quad_grid_pts - self._spans_l = self._derham.spline_attributes[self._space_key].quad_grid_spans - self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts - self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases - self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + """ + Assembles the contribution of a single face to the H(curl) boundary mass matrix. + Calls surface_kernel_3d_mat_h1 for each nonzero (mu, nu) block. + """ + normal_dir = face_idx % 3 + surf_dirs = [d for d in range(3) if d != normal_dir] - self._V_extraction_op = self._derham.extraction_ops[self._space_key] - self._W_extraction_op = self._derham.extraction_ops[self._space_key] - self._V_boundary_op = self._derham.boundary_ops[self._space_key] - self._W_boundary_op = self._derham.boundary_ops[self._space_key] + R_n = self._surface_R_n[face_idx] # shape (nq1, nq2, 3, 3) - self._V_extraction_op_T = self._V_extraction_op.T - self._W_extraction_op_T = self._W_extraction_op.T - self._V_boundary_op_T = self._V_boundary_op.T - self._W_boundary_op_T = self._W_boundary_op.T + for mu, nu, n_comp in [(0, 1, 2), (0, 2, 1), (1, 2, 0)]: + if mu == normal_dir or nu == normal_dir: + continue - # TODO: initialize BlockLinearOperator (3x3 blocks, only diagonal nonzero) - self._mat = None - self._M = None - self._M0 = None - self._domain = None - self._codomain = None - self._dtype = None + fem_space_mu = self._tensor_fem_spaces[mu] + fem_space_nu = self._tensor_fem_spaces[nu] + + starts_mu = [int(s) for s in fem_space_mu.coeff_space.starts] + ends_mu = [int(e) for e in fem_space_mu.coeff_space.ends] + pads_mu = fem_space_mu.coeff_space.pads + + starts_nu = [int(s) for s in fem_space_nu.coeff_space.starts] + ends_nu = [int(e) for e in fem_space_nu.coeff_space.ends] + pads_nu = fem_space_nu.coeff_space.pads + + boundary_index_mu = starts_mu[normal_dir] if face_idx < 3 else ends_mu[normal_dir] + boundary_index_nu = starts_nu[normal_dir] if face_idx < 3 else ends_nu[normal_dir] + + mat_fun_mu_nu = R_n[..., mu, nu] # shape (nq1, nq2) + mat_fun_nu_mu = R_n[..., nu, mu] # = -mat_fun_mu_nu + + # assemble (mu, nu) block: row=mu, col=nu + self._assembly_kernel( + *self._surface_spans[face_idx][mu], + *fem_space_mu.degree, + *fem_space_nu.degree, + *starts_mu, + *pads_mu, + *self._surface_wts[face_idx][mu], + *self._surface_bases[face_idx][mu], + *self._surface_bases[face_idx][nu], + boundary_index_mu, + normal_dir, + mat_fun_mu_nu, + mat.blocks[mu][nu]._data, + ) + + # assemble (nu, mu) block: row=nu, col=mu + self._assembly_kernel( + *self._surface_spans[face_idx][nu], + *fem_space_nu.degree, + *fem_space_mu.degree, + *starts_nu, + *pads_nu, + *self._surface_wts[face_idx][nu], + *self._surface_bases[face_idx][nu], + *self._surface_bases[face_idx][mu], + boundary_index_nu, + normal_dir, + mat_fun_nu_mu, + mat.blocks[nu][mu]._data, + ) - self._surface_quad_grid_meshes = [] - self._surface_geom_weights = [] - self._surface_spans = [] - self._surface_wts = [] - self._surface_bases = [] + def assemble(self, clear: bool = True): + """Assembles the H(curl) boundary mass matrix.""" + if clear: + for mu in range(3): + for nu in range(3): + if mu != nu: + self._mat.blocks[mu][nu]._data[:] = 0.0 for face_idx in range(6): if not self._active_faces[face_idx]: - self._surface_quad_grid_meshes.append(None) - self._surface_geom_weights.append(None) - self._surface_spans.append(None) - self._surface_wts.append(None) - self._surface_bases.append(None) continue + self._assemble_face(face_idx, self._mat) - normal_dir = face_idx % 3 - surf_dirs = [d for d in range(3) if d != normal_dir] - fixed_val = 0.0 if face_idx < 3 else 1.0 - - # TODO: use correct component quadrature points for H(div) - surf_pts_1d = [self._quad_grid_pts[normal_dir][d].flatten() for d in surf_dirs] - e_1d = [None, None, None] - e_1d[surf_dirs[0]] = surf_pts_1d[0] - e_1d[surf_dirs[1]] = surf_pts_1d[1] - e_1d[normal_dir] = xp.array([fixed_val]) - - sqrt_g = xp.abs(self._domain_obj.jacobian_det(*e_1d)) - DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = DFinv[..., normal_dir, :] - norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1)) - self._surface_geom_weights.append(xp.squeeze(sqrt_g * norm_DFinv_n)) - - # TODO: surface meshes, spans, wts, bases for normal_dir component only - self._surface_quad_grid_meshes.append(None) - self._surface_spans.append(None) - self._surface_wts.append(None) - self._surface_bases.append(None) - - # TODO: load assembly kernel - self._assembly_kernel = None - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._codomain - - @property - def dtype(self): - return self._dtype + for mu in range(3): + for nu in range(3): + if mu != nu: + self._mat.blocks[mu][nu].exchange_assembly_data() + self._mat.blocks[mu][nu].update_ghost_regions() - def _assemble_face(self, face_idx: int, mat): - # TODO - pass + def dot(self, v, out=None, apply_bc=True): + """Applies the H(curl) boundary mass matrix to a BlockVector.""" + if out is None: + out = self.codomain.zeros() - def assemble(self, clear: bool = True): - # TODO - pass + if apply_bc: + self._V_boundary_op_T.dot(v, out=self._temp_VB) + self._V_extraction_op_T.dot(self._temp_VB, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=self._temp_WB) + self._W_boundary_op.dot(self._temp_WB, out=out) + else: + self._V_extraction_op_T.dot(v, out=self._temp_mat) + self._mat.dot(self._temp_mat, out=self._temp_WE) + self._W_extraction_op.dot(self._temp_WE, out=out) - def dot(self, v, out=None, apply_bc=True): - # TODO - raise NotImplementedError + return out def transpose(self, conjugate=False): return self def toarray(self): - # TODO - raise NotImplementedError + return self._M0.toarray() def tosparse(self): - # TODO - raise NotImplementedError \ No newline at end of file + return self._M0.tosparse() \ No newline at end of file diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 3e8660ff5..3cc2eafd6 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -180,6 +180,41 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-2 +@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize("active_faces, exact", [ + ([True, False, False, False, False, False], 1.0), + ([False, False, False, True, False, False], -1.0), +]) +def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, exact): + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + u0 = lambda e1, e2, e3: xp.zeros_like(e1) + u1 = lambda e1, e2, e3: xp.ones_like(e1) + u2 = lambda e1, e2, e3: xp.zeros_like(e1) + + v0 = lambda e1, e2, e3: xp.zeros_like(e1) + v1 = lambda e1, e2, e3: xp.zeros_like(e1) + v2 = lambda e1, e2, e3: xp.ones_like(e1) + + P = L2Projector("Hcurl", mass_ops) + u_h = P([u0, u1, u2]) + v_h = P([v0, v1, v2]) + + bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) + numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) + + print(f"active_faces={active_faces}: numerical = {numerical}, exact = {exact}") + + if __name__ == "__main__": from struphy import set_logging_level set_logging_level(logging.INFO) @@ -189,7 +224,7 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), ) - + test_boundary_mass_unit_cube_nonconstant( [8, 9, 10], [1, 2, 3], @@ -205,4 +240,11 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): [8, 9, 10], [1, 2, 3], (("free", "free"), None, ("free", "free")), - ) \ No newline at end of file + ) + + test_boundary_mass_hcurl_per_face( + [16, 16, 16], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + [True, False, False, False, False, False], + 1.0) \ No newline at end of file From 3d1bc6c59728d15ec14486b265831963afc44750 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 3 Aug 2026 16:17:04 +0000 Subject: [PATCH 11/17] Fix H(curl) unit test. --- .../feec/tests/test_boundary_integrals.py | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 3cc2eafd6..39af836cc 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -180,8 +180,8 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-2 -@pytest.mark.parametrize("num_elements", [[8, 8, 8]]) -@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) +@pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) @pytest.mark.parametrize("active_faces, exact", [ ([True, False, False, False, False, False], 1.0), @@ -212,7 +212,10 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, e bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) - print(f"active_faces={active_faces}: numerical = {numerical}, exact = {exact}") + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-1 + if __name__ == "__main__": @@ -243,8 +246,9 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, e ) test_boundary_mass_hcurl_per_face( - [16, 16, 16], - [2, 2, 2], - (("free", "free"), ("free", "free"), ("free", "free")), - [True, False, False, False, False, False], - 1.0) \ No newline at end of file + [16, 16, 16], + [2, 2, 2], + (("free", "free"), ("free", "free"), ("free", "free")), + [True, False, False, False, False, False], + 1.0 + ) \ No newline at end of file From f05cc38f9451d0cde8620d6bcd3db6032367fced Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 4 Aug 2026 12:17:56 +0000 Subject: [PATCH 12/17] Remove metric coefficients from BoundaryMassOperatorHCurl and add general cuboid unit test. --- src/struphy/feec/boundary_mass.py | 152 ++++++++---------- .../feec/tests/test_boundary_integrals.py | 100 ++++++++++-- 2 files changed, 154 insertions(+), 98 deletions(-) diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 9eb83d11c..5cf2cf652 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -467,36 +467,30 @@ def __init__( normal_dir = face_idx % 3 surf_dirs = [d for d in range(3) if d != normal_dir] - fixed_val = 0.0 if face_idx < 3 else 1.0 - # compute R_n using surf_dirs[0] quadrature points - surf_pts_1d = [self._quad_grid_pts[surf_dirs[0]][d].flatten() for d in surf_dirs] - e_1d = [None, None, None] - e_1d[surf_dirs[0]] = surf_pts_1d[0] - e_1d[surf_dirs[1]] = surf_pts_1d[1] - e_1d[normal_dir] = xp.array([fixed_val]) + sign = 1.0 if face_idx < 3 else -1.0 + n_hat = xp.zeros(3) + n_hat[normal_dir] = sign + + # constant skew-symmetric cross-product matrix R_n such that R_n v = n_hat x v + R_n_const = xp.zeros((3, 3)) + R_n_const[0, 1] = -n_hat[2] + R_n_const[0, 2] = n_hat[1] + R_n_const[1, 0] = n_hat[2] + R_n_const[1, 2] = -n_hat[0] + R_n_const[2, 0] = -n_hat[1] + R_n_const[2, 1] = n_hat[0] + + # store R_n per component mu on its own quadrature grid shape + surface_R_n_per_mu = [None, None, None] + for mu in surf_dirs: + nq1 = self._spans_l[mu][surf_dirs[0]].size * self._wts_l[mu][surf_dirs[0]].shape[1] + nq2 = self._spans_l[mu][surf_dirs[1]].size * self._wts_l[mu][surf_dirs[1]].shape[1] + R_n_mu = xp.zeros((nq1, nq2, 3, 3)) + R_n_mu[..., :, :] = R_n_const + surface_R_n_per_mu[mu] = R_n_mu - DFinv = self._domain_obj.jacobian_inv(*e_1d, change_out_order=True) - DFinv_n = xp.squeeze(DFinv[..., normal_dir, :]) - norm_DFinv_n = xp.sqrt(xp.sum(DFinv_n**2, axis=-1, keepdims=True)) - n_phys = DFinv_n / norm_DFinv_n - if face_idx >= 3: - n_phys = -n_phys - - n0 = n_phys[..., 0] - n1 = n_phys[..., 1] - n2 = n_phys[..., 2] - - R_n = xp.zeros((*n0.shape, 3, 3)) - R_n[..., 0, 1] = -n2 - R_n[..., 0, 2] = n1 - R_n[..., 1, 0] = n2 - R_n[..., 1, 2] = -n0 - R_n[..., 2, 0] = -n1 - R_n[..., 2, 1] = n0 - self._surface_R_n.append(R_n) - - # only extract for tangential components mu != normal_dir + self._surface_R_n.append(surface_R_n_per_mu) surface_spans_per_mu = [None, None, None] surface_wts_per_mu = [None, None, None] surface_bases_per_mu = [None, None, None] @@ -527,67 +521,57 @@ def dtype(self): return self._dtype def _assemble_face(self, face_idx: int, mat): - """ - Assembles the contribution of a single face to the H(curl) boundary mass matrix. - Calls surface_kernel_3d_mat_h1 for each nonzero (mu, nu) block. - """ normal_dir = face_idx % 3 surf_dirs = [d for d in range(3) if d != normal_dir] - R_n = self._surface_R_n[face_idx] # shape (nq1, nq2, 3, 3) + mu, nu = surf_dirs[0], surf_dirs[1] - for mu, nu, n_comp in [(0, 1, 2), (0, 2, 1), (1, 2, 0)]: - if mu == normal_dir or nu == normal_dir: - continue + fem_space_mu = self._tensor_fem_spaces[mu] + fem_space_nu = self._tensor_fem_spaces[nu] + + starts_mu = [int(s) for s in fem_space_mu.coeff_space.starts] + ends_mu = [int(e) for e in fem_space_mu.coeff_space.ends] + pads_mu = fem_space_mu.coeff_space.pads - fem_space_mu = self._tensor_fem_spaces[mu] - fem_space_nu = self._tensor_fem_spaces[nu] - - starts_mu = [int(s) for s in fem_space_mu.coeff_space.starts] - ends_mu = [int(e) for e in fem_space_mu.coeff_space.ends] - pads_mu = fem_space_mu.coeff_space.pads - - starts_nu = [int(s) for s in fem_space_nu.coeff_space.starts] - ends_nu = [int(e) for e in fem_space_nu.coeff_space.ends] - pads_nu = fem_space_nu.coeff_space.pads - - boundary_index_mu = starts_mu[normal_dir] if face_idx < 3 else ends_mu[normal_dir] - boundary_index_nu = starts_nu[normal_dir] if face_idx < 3 else ends_nu[normal_dir] - - mat_fun_mu_nu = R_n[..., mu, nu] # shape (nq1, nq2) - mat_fun_nu_mu = R_n[..., nu, mu] # = -mat_fun_mu_nu - - # assemble (mu, nu) block: row=mu, col=nu - self._assembly_kernel( - *self._surface_spans[face_idx][mu], - *fem_space_mu.degree, - *fem_space_nu.degree, - *starts_mu, - *pads_mu, - *self._surface_wts[face_idx][mu], - *self._surface_bases[face_idx][mu], - *self._surface_bases[face_idx][nu], - boundary_index_mu, - normal_dir, - mat_fun_mu_nu, - mat.blocks[mu][nu]._data, - ) - - # assemble (nu, mu) block: row=nu, col=mu - self._assembly_kernel( - *self._surface_spans[face_idx][nu], - *fem_space_nu.degree, - *fem_space_mu.degree, - *starts_nu, - *pads_nu, - *self._surface_wts[face_idx][nu], - *self._surface_bases[face_idx][nu], - *self._surface_bases[face_idx][mu], - boundary_index_nu, - normal_dir, - mat_fun_nu_mu, - mat.blocks[nu][mu]._data, - ) + starts_nu = [int(s) for s in fem_space_nu.coeff_space.starts] + ends_nu = [int(e) for e in fem_space_nu.coeff_space.ends] + pads_nu = fem_space_nu.coeff_space.pads + + boundary_index_mu = starts_mu[normal_dir] if face_idx < 3 else ends_mu[normal_dir] + boundary_index_nu = starts_nu[normal_dir] if face_idx < 3 else ends_nu[normal_dir] + + mat_fun_mu_nu = self._surface_R_n[face_idx][mu][..., mu, nu] + mat_fun_nu_mu = self._surface_R_n[face_idx][nu][..., nu, mu] + + self._assembly_kernel( + *self._surface_spans[face_idx][mu], + *fem_space_mu.degree, + *fem_space_nu.degree, + *starts_mu, + *pads_mu, + *self._surface_wts[face_idx][mu], + *self._surface_bases[face_idx][mu], + *self._surface_bases[face_idx][nu], + boundary_index_mu, + normal_dir, + mat_fun_mu_nu, + mat.blocks[mu][nu]._data, + ) + + self._assembly_kernel( + *self._surface_spans[face_idx][nu], + *fem_space_nu.degree, + *fem_space_mu.degree, + *starts_nu, + *pads_nu, + *self._surface_wts[face_idx][nu], + *self._surface_bases[face_idx][nu], + *self._surface_bases[face_idx][mu], + boundary_index_nu, + normal_dir, + mat_fun_nu_mu, + mat.blocks[nu][mu]._data, + ) def assemble(self, clear: bool = True): """Assembles the H(curl) boundary mass matrix.""" diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 39af836cc..b41d82d19 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -180,39 +180,101 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): assert xp.abs(numerical - exact) < 1e-2 +@pytest.mark.parametrize("num_elements", [[10, 10, 10]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) +@pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize("active_faces, u_idx, v_idx, exact", [ + ([True, False, False, False, False, False], 1, 2, 1.0), + ([False, True, False, False, False, False], 2, 0, 1.0), + ([False, False, True, False, False, False], 0, 1, 1.0), + ([False, False, False, True, False, False], 1, 2, -1.0), + ([False, False, False, False, True, False], 2, 0, -1.0), + ([False, False, False, False, False, True], 0, 1, -1.0), +]) +def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): + comm = MPI.COMM_WORLD + + grid = TensorProductGrid(num_elements=num_elements) + derham_opts = DerhamOptions(degree=degree, bcs=bcs) + derham = Derham(grid, derham_opts, comm=comm) + + domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + mass_ops = WeightedMassOperators(derham, domain) + + u_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == u_idx else xp.zeros_like(e1), + ] + + v_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == v_idx else xp.zeros_like(e1), + ] + + P = L2Projector("Hcurl", mass_ops) + u_h = P(u_funs) + v_h = P(v_funs) + + bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) + numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) + + print(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + + assert xp.abs(numerical - exact) < 1e-1 + + @pytest.mark.parametrize("num_elements", [[8, 9, 10]]) @pytest.mark.parametrize("degree", [[1, 2, 3]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -@pytest.mark.parametrize("active_faces, exact", [ - ([True, False, False, False, False, False], 1.0), - ([False, False, False, True, False, False], -1.0), +@pytest.mark.parametrize("active_faces, u_idx, v_idx, exact", [ + ([True, False, False, False, False, False], 1, 2, 12.0), + ([False, True, False, False, False, False], 2, 0, 6.0), + ([False, False, True, False, False, False], 0, 1, 8.0), + ([False, False, False, True, False, False], 1, 2, -12.0), + ([False, False, False, False, True, False], 2, 0, -6.0), + ([False, False, False, False, False, True ], 0, 1, -8.0), ]) -def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, exact): +def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): + """ + Tests the H(curl) boundary mass operator on a non-unit cuboid [-1,1] x [-1,3] x [0,3] + with constant unit vector fields u = e_{u_idx} and v = e_{v_idx}. + """ comm = MPI.COMM_WORLD grid = TensorProductGrid(num_elements=num_elements) derham_opts = DerhamOptions(degree=degree, bcs=bcs) derham = Derham(grid, derham_opts, comm=comm) - domain = domains.Cuboid(l1=0.0, r1=1.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + domain = domains.Cuboid(l1=-1.0, r1=1.0, l2=-1.0, r2=3.0, l3=0.0, r3=3.0) mass_ops = WeightedMassOperators(derham, domain) - u0 = lambda e1, e2, e3: xp.zeros_like(e1) - u1 = lambda e1, e2, e3: xp.ones_like(e1) - u2 = lambda e1, e2, e3: xp.zeros_like(e1) + u_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == u_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == u_idx else xp.zeros_like(e1), + ] - v0 = lambda e1, e2, e3: xp.zeros_like(e1) - v1 = lambda e1, e2, e3: xp.zeros_like(e1) - v2 = lambda e1, e2, e3: xp.ones_like(e1) + v_funs = [ + lambda e1, e2, e3: xp.ones_like(e1) if 0 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 1 == v_idx else xp.zeros_like(e1), + lambda e1, e2, e3: xp.ones_like(e1) if 2 == v_idx else xp.zeros_like(e1), + ] P = L2Projector("Hcurl", mass_ops) - u_h = P([u0, u1, u2]) - v_h = P([v0, v1, v2]) + u_h = P(u_funs) + v_h = P(v_funs) + + print(f"u_h P1 coeffs: {u_h.toarray()[:5]}") + u_h_l2 = P(u_funs) + print(f"u_h L2 coeffs: {u_h_l2.toarray()[:5]}") + print(f"ratio: {u_h.toarray()[u_h.toarray() != 0][:5] / u_h_l2.toarray()[u_h_l2.toarray() != 0][:5]}") bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) - logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + print(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") assert xp.abs(numerical - exact) < 1e-1 @@ -250,5 +312,15 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, e [2, 2, 2], (("free", "free"), ("free", "free"), ("free", "free")), [True, False, False, False, False, False], + 1, 2, 1.0 + ) + + test_boundary_mass_hcurl_cuboid_nontrivial( + [8, 9, 10], + [1, 2, 3], + (("free", "free"), ("free", "free"), ("free", "free")), + [True, False, False, False, False, False], + 1, 2, + 12.0, ) \ No newline at end of file From bfe8c6c915027c4f8af663f1f8e745a43717ac3b Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 4 Aug 2026 10:43:27 +0200 Subject: [PATCH 13/17] add docstrings to mass_kernels; update feectools commit --- feectools | 2 +- src/struphy/feec/mass_kernels.py | 401 +++++++++++++++++++++++++++++- src/struphy/feec/psydac_derham.py | 16 +- 3 files changed, 403 insertions(+), 16 deletions(-) diff --git a/feectools b/feectools index 8c88dec79..9abb31730 160000 --- a/feectools +++ b/feectools @@ -1 +1 @@ -Subproject commit 8c88dec79510b315024d4b7e0ccc28e76ad8c9e7 +Subproject commit 9abb31730b281d482bcfd8061bf4475c6cfb4ade diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 7d54d5574..e3c867f0f 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -21,13 +21,39 @@ def kernel_1d_mat( data: "float[:,:]", ): """ - Performs the integration of Lambda_i * mat_fun(eta1) * Lambda_j for the basis functions (i, j) available on the calling process. + Performs the integration of Lambda_(i1) * mat_fun(eta1) * Lambda_(j1) for the basis functions (i1, j1) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the codomain basis functions. + pj1 : int + Degree of the domain basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in data. + w1 : "float[:,:]" + Quadrature weights. The indexing is [global element, quadrature point]. + bi1 : "float[:,:,:,:]" + Values of codomain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + bj1 : "float[:,:,:,:]" + Values of domain basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:]" + Function under the integral evaluated at quadrature points (flattened). + data : "float[:,:]" + _data array of StencilMatrix to store the results. """ + # number of elements ne1 = spans1.size + # number of quadrature points in each element nq1 = shape(w1)[1] for iel1 in range(ne1): @@ -58,9 +84,29 @@ def kernel_1d_vec( data: "float[:]", ): """ - Performs the integration of Lambda_i * mat_fun(eta1) for the basis functions (i) available on the calling process. + Performs the integration of Lambda_(i1) * mat_fun(eta1) for the basis functions (i1) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in data. + w1 : "float[:,:]" + Quadrature weights. The indexing is [global element, quadrature point]. + bi1 : "float[:,:,:,:]" + Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:]" + Function under the integral evaluated at quadrature points (flattened). + data : "float[:]" + _data array of StencilVector to store the results. """ ne1 = spans1.size @@ -93,9 +139,28 @@ def kernel_1d_eval( values: "float[:]", ): """ - Evaluates sum_i [ coeffs_i * Lambda_i(quad_eta1) ] for all quadrature points on the calling process. + Evaluates sum_i1 [ coeffs_i1 * Lambda_i1(quad_eta1) ] for all quadrature points on the calling process. The results are written into values. + + Parameters + ---------- + spans1 : array[int] + Array of span indices; the span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells). + pi1 : int + Degree of the basis functions. + starts1 : int + Starting index on the current rank. + pads1 : int + Padding (=spline degree) for ghost regions in coeffs_data. + bi1 : "float[:,:,:,:]" + Values of basis functions. The indexing is [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:]" + Output array (flattened over elements and quadrature points) holding the evaluated function values; + it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. """ values[:] = 0.0 @@ -140,9 +205,35 @@ def kernel_2d_mat( data: "float[:,:,:,:]", ): """ - Performs the integration of Lambda_ij * mat_fun(eta1, eta2) * Lambda_lm for the basis functions (ij, lm) available on the calling process. + Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) * Lambda_(j1, j2) for the basis functions (i1, i2, j1, j2) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the codomain basis functions in direction 1 and 2. + pj1, pj2 : int + Degree of the domain basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in direction 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2 : "float[:,:,:,:]" + Values of domain basis functions in direction 1 and 2, same indexing convention as bi1, bi2. + mat_fun : "float[:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. + data : "float[:,:,:,:]" + _data array of StencilMatrix to store the results. """ ne1 = spans1.size @@ -195,9 +286,31 @@ def kernel_2d_vec( data: "float[:,:]", ): """ - Performs the integration of Lambda_ij * mat_fun(eta1, eta2) for the basis functions (ij) available on the calling process. + Performs the integration of Lambda_(i1, i2) * mat_fun(eta1, eta2) for the basis functions (i1, i2) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in direction 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in direction 1 and 2. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quadrature point in direction 1, flattened quadrature point in direction 2]. + data : "float[:,:]" + _data array of StencilVector to store the results. """ ne1 = spans1.size @@ -244,9 +357,29 @@ def kernel_2d_eval( values: "float[:,:]", ): """ - Evaluates sum_ij [ coeffs_ij * Lambda_ij(quad_eta1, quad_eta2) ] for all quadrature points on the calling process. + Evaluates sum_(i1, i2) [ coeffs_{i1,i2} * Lambda_{i1, i2}(quad_eta1, quad_eta2) ] for all quadrature points on the calling process. The results are written into values. + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in direction 1 and 2; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2 : int + Degree of the basis functions in direction 1 and 2. + starts1, starts2 : int + Starting index on the current rank, in direction 1 and 2. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in coeffs_data, in direction 1 and 2. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in direction 1 and 2. The indexing is + [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:,:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:,:]" + Output array (flattened over elements and quadrature points in each direction) holding the evaluated + function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. """ values[:, :] = 0.0 @@ -310,9 +443,35 @@ def kernel_3d_mat( data: "float[:,:,:,:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) * Lambda_lmn for the basis functions (ijk, lmn) available on the calling process. + Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) * Lambda_(j1,j2,j3) for the basis functions (i1,i2,i3, j1,j2,j3) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the codomain basis functions in direction 1, 2 and 3. + pj1, pj2, pj3 : int + Degree of the domain basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2, bj3 : "float[:,:,:,:]" + Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results. """ ne1 = spans1.size @@ -421,9 +580,31 @@ def kernel_3d_vec( data: "float[:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) for the basis functions (ijk) available on the calling process. + Performs the integration of Lambda_(i1,i2,i3) * mat_fun(eta1, eta2, eta3) for the basis functions (i1,i2,i3) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in data, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:]" + _data array of StencilVector to store the results. """ ne1 = spans1.size @@ -489,9 +670,29 @@ def kernel_3d_eval( values: "float[:,:,:]", ): """ - Evaluates sum_ijk [ coeffs_ijk * Lambda_ijk(quad_eta1, quad_eta2, quad_eta3) ] for all quadrature points on the calling process. + Evaluates sum_(i1,i2,i3) [ coeffs_{i1,i2,i3} * Lambda_{i1,i2,i3}(quad_eta1, quad_eta2, quad_eta3) ] for all quadrature points on the calling process. The results are written into values. + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions in coeffs_data, in direction 1, 2 and 3. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + coeffs_data : "float[:,:,:]" + _data array of StencilVector holding the spline coefficients of the function to be evaluated. + values : "float[:,:,:]" + Output array (flattened over elements and quadrature points in each direction) holding the evaluated + function values; it is set to zero at the start of the kernel, i.e. it is overwritten, not added to. """ values[:, :, :] = 0.0 @@ -570,10 +771,45 @@ def kernel_3d_matrixfree( data_in: "float[:,:,:]", ): """ - Performs the integration of Lambda_ijk * mat_fun(eta1, eta2, eta3) * f(eta1, eta2, eta3) for the basis functions (ijk) available on the calling process, + Performs the integration of Lambda_(i1, i2, i3) * mat_fun(eta1, eta2, eta3) * f(eta1, eta2, eta3) for the basis functions (i1, i2, i3) available on the calling process, where f is the spline function represented by the coefficients in data_in. - The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + The results are written into data_out (attention: data_out is NOT set to zero first, but the results are added to data_out). + This computes the action of the mass matrix on a vector without ever assembling the matrix itself. + + Parameters + ---------- + spansi1, spansi2, spansi3 : array[int] + Arrays of span indices in direction 1, 2 and 3 for the codomain ("i") basis functions; the span is the + index of the last non-vanishing spline on each grid element (cell). + spansj1, spansj2, spansj3 : array[int] + Arrays of span indices in direction 1, 2 and 3 for the domain ("j") basis functions. + pi1, pi2, pi3 : int + Degree of the codomain basis functions in direction 1, 2 and 3. + pj1, pj2, pj3 : int + Degree of the domain basis functions in direction 1, 2 and 3. + startsi1, startsi2, startsi3 : int + Starting index on the current rank for the codomain basis functions, in direction 1, 2 and 3. + startsj1, startsj2, startsj3 : int + Starting index on the current rank for the domain basis functions, in direction 1, 2 and 3. + padsi1, padsi2, padsi3 : int + Padding (=spline degree) for ghost regions in data_out, in direction 1, 2 and 3. + padsj1, padsj2, padsj3 : int + Padding (=spline degree) for ghost regions in data_in, in direction 1, 2 and 3. + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of codomain basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2, bj3 : "float[:,:,:,:]" + Values of domain basis functions in direction 1, 2 and 3, same indexing convention as bi1, bi2, bi3. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data_out : "float[:,:,:]" + _data array of StencilVector to store the results of the matrix-vector product. + data_in : "float[:,:,:]" + _data array of StencilVector holding the spline coefficients of the input function f. """ ne1 = spansi1.size @@ -691,6 +927,30 @@ def kernel_3d_diag( Computes the diagonal of a mass matrix, assuming that the domain and the codomain are the same. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2, spans3 : array[int] + Arrays of span indices in direction 1, 2 and 3; the span is the index of the last non-vanishing spline + on each grid element (cell). The length of each array is the number of elements (cells) in that direction. + pi1, pi2, pi3 : int + Degree of the basis functions in direction 1, 2 and 3. + starts1, starts2, starts3 : int + Starting index on the current rank, in direction 1, 2 and 3. + pads1, pads2, pads3 : int + Padding (=spline degree) for ghost regions, in direction 1, 2 and 3 (unused for data, which is a + StencilDiagonalMatrix and therefore has no padding, but kept for a uniform kernel signature). + w1, w2, w3 : "float[:,:]" + Quadrature weights in direction 1, 2 and 3. The indexing is [global element, quadrature point]. + bi1, bi2, bi3 : "float[:,:,:,:]" + Values of basis functions in direction 1, 2 and 3. The indexing is + [global element, local basis function, derivative, quadrature point]. + mat_fun : "float[:,:,:]" + Function under the integral evaluated at quadrature points (flattened in each direction). + The indexing is [flattened quad. point dir. 1, flattened quad. point dir. 2, flattened quad. point dir. 3]. + data : "float[:,:,:]" + _data array of StencilDiagonalMatrix to store the results. Periodic wrap-around (index -= nb) is applied + when a local index runs beyond the array bounds, since there are no ghost regions on this matrix type. """ ne1 = spans1.size @@ -788,6 +1048,45 @@ def surface_kernel_3d_vec( mat_fun: "float[:,:]", data: "float[:,:,:]", ): + """ + Performs the integration of Lambda_0ij * mat_fun(eta1, eta2) over the boundary surface at the fixed + (normal-direction) global index boundary_index, for the basis functions (ij) available on the calling + process in the two surface (tangential) directions. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in the two surface (tangential) directions; the span is the index of the last + non-vanishing spline on each grid element (cell) in that direction. + pi0 : int + Degree of the basis function in the normal direction (kept for a uniform kernel signature; not used + directly since the normal index is fixed to boundary_index). + pi1, pi2 : int + Degree of the basis functions in the two surface directions. + starts0 : int + Starting index on the current rank in the normal direction. + starts1, starts2 : int + Starting index on the current rank in the two surface directions. + pads0 : int + Padding (=spline degree) for ghost regions in data, in the normal direction. + pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, in the two surface directions. + w1, w2 : "float[:,:]" + Quadrature weights in the two surface directions. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of basis functions in the two surface directions. The indexing is + [global element, local basis function, derivative, quadrature point]. + boundary_index : int + Global index in the normal direction at which the boundary surface is located. + mat_fun : "float[:,:]" + Function under the integral evaluated at surface quadrature points (flattened in each surface direction). + data : "float[:,:,:]" + _data array of StencilVector to store the results; only the slice at the fixed normal index + (pads0 + i_local0) is written. + """ + ne1 = spans1.size ne2 = spans2.size @@ -849,6 +1148,48 @@ def surface_kernel_3d_mat_h1( mat_fun: "float[:,:]", data: "float[:,:,:,:,:,:]", ): + """ + Computes a boundary (surface) H1 mass matrix: the integration of Lambda_i * mat_fun(eta_s1, eta_s2) * Lambda_j + over the boundary surface at the fixed global index boundary_index in the normal_dir direction, for the + codomain ("i") and domain ("j") basis functions available on the calling process in the two tangential + directions orthogonal to normal_dir. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in the two tangential grid directions used for the surface quadrature (as + determined by normal_dir); the span is the index of the last non-vanishing spline on each grid element. + pi0, pi1, pi2 : int + Degree of the codomain basis functions along Cartesian axes 0, 1 and 2. + qi0, qi1, qi2 : int + Degree of the domain basis functions along Cartesian axes 0, 1 and 2. + starts0, starts1, starts2 : int + Starting index on the current rank along Cartesian axes 0, 1 and 2. + pads0, pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, along Cartesian axes 0, 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in the two tangential directions. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of codomain basis functions in the two tangential directions. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2 : "float[:,:,:,:]" + Values of domain basis functions in the two tangential directions, same indexing convention as bi1, bi2. + boundary_index : int + Global index along normal_dir at which the boundary surface is located. + normal_dir : int + Cartesian direction (0, 1 or 2) normal to the surface; the remaining two directions are the tangential + directions used for the surface integration. + mat_fun : "float[:,:]" + Function under the integral evaluated at surface quadrature points (flattened in each tangential direction). + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results. Rows are fixed at the boundary index in normal_dir; + columns are offset by (j_global - i_global) in the tangential directions and fixed to zero offset in + normal_dir, following the standard banded StencilMatrix storage convention. The axis order of the + stored indices depends on normal_dir (0, 1 or 2). + """ + ne1 = spans1.size ne2 = spans2.size @@ -968,4 +1309,42 @@ def surface_kernel_3d_mat_hcurl( n_cross_weight: "float[:,:]", data: "float[:,:,:,:,:,:]", ): + """ + Computes a boundary (surface) H(curl) mass matrix: the integration of (n x Lambda_i) * mat_fun(eta_s1, eta_s2) + * (n x Lambda_j), weighted by n_cross_weight, over the boundary surface at the fixed global index + boundary_index, for the codomain ("i") and domain ("j") tangential-trace basis functions available on the + calling process. Mirrors surface_kernel_3d_mat_h1, but for H(curl) spaces where only the tangential + components of the basis functions couple through the cross product with the surface normal n. + + The results are written into data (attention: data is NOT set to zero first, but the results are added to data). + + Parameters + ---------- + spans1, spans2 : array[int] + Arrays of span indices in the two tangential grid directions used for the surface quadrature; the span + is the index of the last non-vanishing spline on each grid element. + pi0, pi1, pi2 : int + Degree of the codomain basis functions along Cartesian axes 0, 1 and 2. + qi0, qi1, qi2 : int + Degree of the domain basis functions along Cartesian axes 0, 1 and 2. + starts0, starts1, starts2 : int + Starting index on the current rank along Cartesian axes 0, 1 and 2. + pads0, pads1, pads2 : int + Padding (=spline degree) for ghost regions in data, along Cartesian axes 0, 1 and 2. + w1, w2 : "float[:,:]" + Quadrature weights in the two tangential directions. The indexing is [global element, quadrature point]. + bi1, bi2 : "float[:,:,:,:]" + Values of codomain basis functions in the two tangential directions. The indexing is + [global element, local basis function, derivative, quadrature point]. + bj1, bj2 : "float[:,:,:,:]" + Values of domain basis functions in the two tangential directions, same indexing convention as bi1, bi2. + boundary_index : int + Global index in the normal direction at which the boundary surface is located. + n_cross_weight : "float[:,:]" + Weight from the cross product with the surface normal n (and mat_fun) evaluated at surface quadrature + points, indexed like mat_fun in surface_kernel_3d_mat_h1. + data : "float[:,:,:,:,:,:]" + _data array of StencilMatrix to store the results, following the same storage convention as in + surface_kernel_3d_mat_h1. + """ pass \ No newline at end of file diff --git a/src/struphy/feec/psydac_derham.py b/src/struphy/feec/psydac_derham.py index 82bc49377..631827e41 100644 --- a/src/struphy/feec/psydac_derham.py +++ b/src/struphy/feec/psydac_derham.py @@ -507,22 +507,30 @@ def proj_loc_grid_wts(self) -> tuple[tuple[xp.ndarray]] | None: @property def quad_grid_pts(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid points in each direction for each component of the vector space.""" + """Tuple of quadrature grid points in each direction for each component of the vector space. + The indexing is [global element, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_pts @property def quad_grid_wts(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid weights in each direction for each component of the vector space.""" + """Tuple of quadrature grid weights in each direction for each component of the vector space. + The indexing is [global element, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_wts @property def quad_grid_spans(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid basis function spans in each direction for each component of the vector space.""" + """Tuple of quadrature grid basis function spans in each direction for each component of the vector space. + The span is the index of the last non-vanishing spline on each grid element + (cell). The length of the returned array is the number of elements (cells).""" return self._quad_grid_spans @property def quad_grid_bases(self) -> tuple[tuple[xp.ndarray]]: - """Tuple of quadrature grid basis function values in each direction for each component of the vector space.""" + """Tuple of quadrature grid basis function values in each direction for each component of the vector space. + Indexing is [global element, basis function, derivative, local quadrature point]. + The length of the first dimension is the number of elements (cells).""" return self._quad_grid_bases From ea4ef1dd9394397a4cfe56acb08ff0f5b7af2fd3 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 4 Aug 2026 11:44:17 +0200 Subject: [PATCH 14/17] make column index more clear, consistent with other kernels --- src/struphy/feec/mass_kernels.py | 40 ++++++++++++++------------------ 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index e3c867f0f..1a5864fb8 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -1128,9 +1128,9 @@ def surface_kernel_3d_mat_h1( pi0: int, pi1: int, pi2: int, - qi0: int, - qi1: int, - qi2: int, + pj0: int, + pj1: int, + pj2: int, starts0: int, starts1: int, starts2: int, @@ -1162,13 +1162,13 @@ def surface_kernel_3d_mat_h1( Arrays of span indices in the two tangential grid directions used for the surface quadrature (as determined by normal_dir); the span is the index of the last non-vanishing spline on each grid element. pi0, pi1, pi2 : int - Degree of the codomain basis functions along Cartesian axes 0, 1 and 2. - qi0, qi1, qi2 : int - Degree of the domain basis functions along Cartesian axes 0, 1 and 2. + Degree of the codomain basis functions along logical axes 0, 1 and 2. + pj0, pj1, pj2 : int + Degree of the domain basis functions along logical axes 0, 1 and 2. starts0, starts1, starts2 : int - Starting index on the current rank along Cartesian axes 0, 1 and 2. + Starting index on the current rank along logical axes 0, 1 and 2. pads0, pads1, pads2 : int - Padding (=spline degree) for ghost regions in data, along Cartesian axes 0, 1 and 2. + Padding (=spline degree) for ghost regions in data, along logical axes 0, 1 and 2. w1, w2 : "float[:,:]" Quadrature weights in the two tangential directions. The indexing is [global element, quadrature point]. bi1, bi2 : "float[:,:,:,:]" @@ -1179,15 +1179,12 @@ def surface_kernel_3d_mat_h1( boundary_index : int Global index along normal_dir at which the boundary surface is located. normal_dir : int - Cartesian direction (0, 1 or 2) normal to the surface; the remaining two directions are the tangential + Logical direction (0, 1 or 2) normal to the surface; the remaining two directions are the tangential directions used for the surface integration. mat_fun : "float[:,:]" Function under the integral evaluated at surface quadrature points (flattened in each tangential direction). data : "float[:,:,:,:,:,:]" - _data array of StencilMatrix to store the results. Rows are fixed at the boundary index in normal_dir; - columns are offset by (j_global - i_global) in the tangential directions and fixed to zero offset in - normal_dir, following the standard banded StencilMatrix storage convention. The axis order of the - stored indices depends on normal_dir (0, 1 or 2). + _data array of StencilMatrix to store the results. """ ne1 = spans1.size @@ -1199,15 +1196,15 @@ def surface_kernel_3d_mat_h1( starts = [starts0, starts1, starts2] pads = [pads0, pads1, pads2] pi = [pi0, pi1, pi2] - qi = [qi0, qi1, qi2] + pj = [pj0, pj1, pj2] surf_dirs = [d for d in range(3) if d != normal_dir] pi_s1 = pi[surf_dirs[0]] pi_s2 = pi[surf_dirs[1]] - qi_s1 = qi[surf_dirs[0]] - qi_s2 = qi[surf_dirs[1]] + pj_s1 = pj[surf_dirs[0]] + pj_s2 = pj[surf_dirs[1]] starts_n = starts[normal_dir] starts_s1 = starts[surf_dirs[0]] @@ -1229,13 +1226,10 @@ def surface_kernel_3d_mat_h1( i_local1 = i_global1 - starts_s1 i_local2 = i_global2 - starts_s2 - for jl1 in range(qi_s1 + 1): - for jl2 in range(qi_s2 + 1): - j_global1 = spans1[iel1] - qi_s1 + jl1 - j_global2 = spans2[iel2] - qi_s2 + jl2 - - j_local1 = j_global1 - i_global1 - j_local2 = j_global2 - i_global2 + for jl1 in range(pj_s1 + 1): + for jl2 in range(pj_s2 + 1): + j_local1 = jl1 - il1 + j_local2 = jl2 - il2 value = 0.0 From 1b37d0c1faadb5634f6cf722d7de1137114f8688 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 4 Aug 2026 14:32:24 +0200 Subject: [PATCH 15/17] make boundary integrals work in parallel (MPI) --- src/struphy/feec/boundary_mass.py | 58 ++++++------- src/struphy/feec/mass_kernels.py | 22 ++--- src/struphy/feec/psydac_derham.py | 4 +- .../feec/tests/test_boundary_integrals.py | 82 ++++++++++++++----- 4 files changed, 95 insertions(+), 71 deletions(-) diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 5cf2cf652..0fbea0b5f 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -3,7 +3,7 @@ import cunumpy as xp from feectools.api.settings import PSYDAC_BACKEND_GPYCCEL -from feectools.linalg.block import BlockVector +from feectools.linalg.block import BlockLinearOperator, BlockVector from feectools.linalg.stencil import StencilMatrix, StencilVector from struphy.feec import mass_kernels @@ -13,7 +13,7 @@ from struphy.geometry.base import Domain from struphy.utils.pyccel import Pyccelkernel -from feectools.linalg.block import BlockLinearOperator, BlockVector +logger = logging.getLogger("struphy") class BoundaryIntegralOperators: @@ -34,7 +34,7 @@ def __init__( mass_ops: WeightedMassOperators, active_faces: list[bool] | None = None, ): - + self._mass_ops = mass_ops self._derham = mass_ops.derham self._domain = mass_ops.domain @@ -114,11 +114,7 @@ class BoundaryMassOperatorH1(LinOpWithTransp): Mass operators object, contains geometry and derham. """ - def __init__( - self, - mass_ops: WeightedMassOperators, - active_faces: list[bool] - ): + def __init__(self, mass_ops: WeightedMassOperators, active_faces: list[bool]): self._mass_ops = mass_ops self._derham = mass_ops.derham self._domain_obj = mass_ops.domain @@ -133,6 +129,7 @@ def __init__( self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + self._nbasis = self._derham.spline_attributes[self._space_key].nbasis # boundary and extraction operators self._V_extraction_op = self._derham.extraction_ops[self._space_key] @@ -253,26 +250,27 @@ def _assemble_face( ends = [int(end) for end in fem_space.coeff_space.ends] pads = fem_space.coeff_space.pads - boundary_index = starts[normal_dir] if face_idx < 3 else ends[normal_dir] - - fem_space = self._tensor_fem_spaces[0] - starts = [int(start) for start in fem_space.coeff_space.starts] - pads = fem_space.coeff_space.pads - - self._assembly_kernel( - *self._surface_spans[face_idx], - *fem_space.degree, - *fem_space.degree, - *starts, - *pads, - *self._surface_wts[face_idx], - *self._surface_bases[face_idx], - *self._surface_bases[face_idx], - boundary_index, - normal_dir, - self._surface_geom_weights[face_idx], - mat._data, - ) + boundary_index = 0 if face_idx < 3 else self._nbasis[0][normal_dir] - 1 + + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index=}, {starts=}, {ends=}, {pads=}") + + # only assemble if current rank is a true boundary (not an interior partition boundary) + if starts[normal_dir] == boundary_index or ends[normal_dir] == boundary_index: + logger.debug("Assembling face", face_idx) + self._assembly_kernel( + *self._surface_spans[face_idx], + *fem_space.degree, + *fem_space.degree, + *starts, + *pads, + *self._surface_wts[face_idx], + *self._surface_bases[face_idx], + *self._surface_bases[face_idx], + boundary_index, + normal_dir, + self._surface_geom_weights[face_idx], + mat._data, + ) def assemble( self, @@ -354,11 +352,9 @@ def __call__( self.assemble(clear=clear) return self - def toarray(self): return self._M0.toarray() - def tosparse(self): return self._M0.tosparse() @@ -617,4 +613,4 @@ def toarray(self): return self._M0.toarray() def tosparse(self): - return self._M0.tosparse() \ No newline at end of file + return self._M0.tosparse() diff --git a/src/struphy/feec/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 1a5864fb8..8a7244918 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -24,7 +24,7 @@ def kernel_1d_mat( Performs the integration of Lambda_(i1) * mat_fun(eta1) * Lambda_(j1) for the basis functions (i1, j1) available on the calling process. The results are written into data (attention: data is NOT set to zero first, but the results are added to data). - + Parameters ---------- spans1 : array[int] @@ -1109,15 +1109,9 @@ def surface_kernel_3d_vec( for q1 in range(nq1): for q2 in range(nq2): - wvol = ( - w1[iel1, q1] - * w2[iel2, q2] - * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - ) + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - value += ( - wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] - ) + value += wvol * bi1[iel1, il1, 0, q1] * bi2[iel2, il2, 0, q2] data[pads0 + i_local0, pads1 + i_local1, pads2 + i_local2] += value @@ -1184,7 +1178,7 @@ def surface_kernel_3d_mat_h1( mat_fun : "float[:,:]" Function under the integral evaluated at surface quadrature points (flattened in each tangential direction). data : "float[:,:,:,:,:,:]" - _data array of StencilMatrix to store the results. + _data array of StencilMatrix to store the results. """ ne1 = spans1.size @@ -1235,11 +1229,7 @@ def surface_kernel_3d_mat_h1( for q1 in range(nq1): for q2 in range(nq2): - wvol = ( - w1[iel1, q1] - * w2[iel2, q2] - * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] - ) + wvol = w1[iel1, q1] * w2[iel2, q2] * mat_fun[iel1 * nq1 + q1, iel2 * nq2 + q2] value += ( wvol @@ -1341,4 +1331,4 @@ def surface_kernel_3d_mat_hcurl( _data array of StencilMatrix to store the results, following the same storage convention as in surface_kernel_3d_mat_h1. """ - pass \ No newline at end of file + pass diff --git a/src/struphy/feec/psydac_derham.py b/src/struphy/feec/psydac_derham.py index 631827e41..e6ee50e60 100644 --- a/src/struphy/feec/psydac_derham.py +++ b/src/struphy/feec/psydac_derham.py @@ -515,7 +515,7 @@ def quad_grid_pts(self) -> tuple[tuple[xp.ndarray]]: @property def quad_grid_wts(self) -> tuple[tuple[xp.ndarray]]: """Tuple of quadrature grid weights in each direction for each component of the vector space. - The indexing is [global element, local quadrature point]. + The indexing is [global element, local quadrature point]. The length of the first dimension is the number of elements (cells).""" return self._quad_grid_wts @@ -529,7 +529,7 @@ def quad_grid_spans(self) -> tuple[tuple[xp.ndarray]]: @property def quad_grid_bases(self) -> tuple[tuple[xp.ndarray]]: """Tuple of quadrature grid basis function values in each direction for each component of the vector space. - Indexing is [global element, basis function, derivative, local quadrature point]. + Indexing is [global element, basis function, derivative, local quadrature point]. The length of the first dimension is the number of elements (cells).""" return self._quad_grid_bases diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index b41d82d19..2239bd140 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -15,14 +15,22 @@ logger = logging.getLogger("struphy") -@pytest.mark.parametrize("num_elements", [[8, 9, 10]],) +@pytest.mark.parametrize( + "num_elements", + [[8, 9, 10]], +) @pytest.mark.parametrize("degree", [[1, 2, 3]]) -@pytest.mark.parametrize("bcs", [(None, None, None), - (("free", "free"), None, None), - (None, ("free", "free"), None), - (None, None, ("free", "free")), - (("free", "free"), ("free", "free"), None), - (("free", "free"), ("free", "free"), ("free", "free")),]) +@pytest.mark.parametrize( + "bcs", + [ + (None, None, None), + (("free", "free"), None, None), + (None, ("free", "free"), None), + (None, None, ("free", "free")), + (("free", "free"), ("free", "free"), None), + (("free", "free"), ("free", "free"), ("free", "free")), + ], +) def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = 1 on the unit cube. @@ -55,8 +63,15 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) - numerical = xp.sum(v.toarray()) + numerical = xp.sum(coeffs) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") @@ -65,9 +80,14 @@ def test_boundary_mass_unit_cube_constant(num_elements, degree, bcs): @pytest.mark.parametrize("num_elements", [[8, 9, 10]]) @pytest.mark.parametrize("degree", [[1, 2, 3]]) -@pytest.mark.parametrize("bcs", [(("dirichlet", "free"), ("free", "free"), ("free", "free")), - (("free", "dirichlet"), ("free", "free"), ("free", "free")), - (("dirichlet", "dirichlet"), ("free", "free"), ("free", "free"))]) +@pytest.mark.parametrize( + "bcs", + [ + (("dirichlet", "free"), ("free", "free"), ("free", "free")), + (("free", "dirichlet"), ("free", "free"), ("free", "free")), + (("dirichlet", "dirichlet"), ("free", "free"), ("free", "free")), + ], +) def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): """ Tests the boundary mass operator for alpha = eta1 + eta2 + eta3 on the unit cube. @@ -97,8 +117,15 @@ def test_boundary_mass_unit_cube_nonconstant(num_elements, degree, bcs): bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() - numerical = xp.sum(v.toarray()) + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) + + numerical = xp.sum(coeffs) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") @@ -130,8 +157,15 @@ def test_boundary_mass_cuboid_nontrivial(num_elements, degree, bcs): bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) - numerical = xp.sum(v.toarray()) + numerical = xp.sum(coeffs) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") @@ -146,6 +180,7 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): Tests the boundary mass operator for alpha = exp(eta3) on a HollowCylinder. """ import math + comm = MPI.COMM_WORLD grid = TensorProductGrid(num_elements=num_elements) @@ -161,19 +196,22 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): alpha = lambda e1, e2, e3: xp.exp(e3) e = math.e - exact = xp.pi * ( - 2 * a1 * Lz * (e - 1) - + 2 * a2 * Lz * (e - 1) - + (a2**2 - a1**2) * (1 + e) - ) + exact = xp.pi * (2 * a1 * Lz * (e - 1) + 2 * a2 * Lz * (e - 1) + (a2**2 - a1**2) * (1 + e)) P = L2Projector("H1", mass_ops) alpha_h = P(alpha) bnd_ops = BoundaryIntegralOperators(mass_ops) v = bnd_ops.S0.dot(alpha_h) + arr = v.toarray() + + if comm is None: + coeffs = arr + else: + coeffs = xp.zeros_like(arr) + comm.Allreduce(arr, coeffs, op=MPI.SUM) - numerical = xp.sum(v.toarray()) + numerical = xp.sum(coeffs) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") @@ -279,9 +317,9 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active assert xp.abs(numerical - exact) < 1e-1 - if __name__ == "__main__": from struphy import set_logging_level + set_logging_level(logging.INFO) test_boundary_mass_unit_cube_constant( @@ -295,7 +333,7 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active [1, 2, 3], (("dirichlet", "free"), ("free", "free"), ("free", "free")), ) - + test_boundary_mass_cuboid_nontrivial( [8, 9, 10], [1, 2, 3], @@ -323,4 +361,4 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active [True, False, False, False, False, False], 1, 2, 12.0, - ) \ No newline at end of file + ) From f0f3aff2f66fc098e5d6d050b32b1d091a503024 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 4 Aug 2026 12:50:35 +0000 Subject: [PATCH 16/17] Add pullback to stretched cuboid H(curl) test. --- .../feec/tests/test_boundary_integrals.py | 47 +++++++++---------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 2239bd140..23dc26e20 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -258,13 +258,13 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) - print(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") assert xp.abs(numerical - exact) < 1e-1 -@pytest.mark.parametrize("num_elements", [[8, 9, 10]]) -@pytest.mark.parametrize("degree", [[1, 2, 3]]) +@pytest.mark.parametrize("num_elements", [[10, 10, 10]]) +@pytest.mark.parametrize("degree", [[2, 2, 2]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) @pytest.mark.parametrize("active_faces, u_idx, v_idx, exact", [ ([True, False, False, False, False, False], 1, 2, 12.0), @@ -288,33 +288,30 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active domain = domains.Cuboid(l1=-1.0, r1=1.0, l2=-1.0, r2=3.0, l3=0.0, r3=3.0) mass_ops = WeightedMassOperators(derham, domain) - u_funs = [ - lambda e1, e2, e3: xp.ones_like(e1) if 0 == u_idx else xp.zeros_like(e1), - lambda e1, e2, e3: xp.ones_like(e1) if 1 == u_idx else xp.zeros_like(e1), - lambda e1, e2, e3: xp.ones_like(e1) if 2 == u_idx else xp.zeros_like(e1), - ] - - v_funs = [ - lambda e1, e2, e3: xp.ones_like(e1) if 0 == v_idx else xp.zeros_like(e1), - lambda e1, e2, e3: xp.ones_like(e1) if 1 == v_idx else xp.zeros_like(e1), - lambda e1, e2, e3: xp.ones_like(e1) if 2 == v_idx else xp.zeros_like(e1), - ] + def make_pulled(domain, idx): + phys_funs = [ + lambda x, y, z: xp.ones_like(x) if 0 == idx else xp.zeros_like(x), + lambda x, y, z: xp.ones_like(x) if 1 == idx else xp.zeros_like(x), + lambda x, y, z: xp.ones_like(x) if 2 == idx else xp.zeros_like(x), + ] + def pulled(*etas): + return domain.pull(phys_funs, *etas, kind="1") + return [ + lambda *etas, p=pulled: p(*etas)[0], + lambda *etas, p=pulled: p(*etas)[1], + lambda *etas, p=pulled: p(*etas)[2], + ] P = L2Projector("Hcurl", mass_ops) - u_h = P(u_funs) - v_h = P(v_funs) - - print(f"u_h P1 coeffs: {u_h.toarray()[:5]}") - u_h_l2 = P(u_funs) - print(f"u_h L2 coeffs: {u_h_l2.toarray()[:5]}") - print(f"ratio: {u_h.toarray()[u_h.toarray() != 0][:5] / u_h_l2.toarray()[u_h_l2.toarray() != 0][:5]}") + u_h = P(make_pulled(domain, u_idx)) + v_h = P(make_pulled(domain, v_idx)) bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) - print(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") + logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") - assert xp.abs(numerical - exact) < 1e-1 + assert xp.abs(numerical - exact) < 1 if __name__ == "__main__": @@ -346,7 +343,7 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active ) test_boundary_mass_hcurl_per_face( - [16, 16, 16], + [10, 10, 10], [2, 2, 2], (("free", "free"), ("free", "free"), ("free", "free")), [True, False, False, False, False, False], @@ -355,7 +352,7 @@ def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active ) test_boundary_mass_hcurl_cuboid_nontrivial( - [8, 9, 10], + [10, 10, 10], [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), [True, False, False, False, False, False], From d340bd13d2993212a8e1cf4c8260b729bc4e33bd Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 4 Aug 2026 15:29:34 +0200 Subject: [PATCH 17/17] fix BoundaryMassOperatorHCurl to run in parallel (MPI); update feectools commt --- feectools | 2 +- src/struphy/feec/boundary_mass.py | 80 +++++++++++-------- .../feec/tests/test_boundary_integrals.py | 50 +++++++----- 3 files changed, 77 insertions(+), 55 deletions(-) diff --git a/feectools b/feectools index f091889ea..9abb31730 160000 --- a/feectools +++ b/feectools @@ -1 +1 @@ -Subproject commit f091889ea41bee9b7cc97dbec7b9f95ce1f4384c +Subproject commit 9abb31730b281d482bcfd8061bf4475c6cfb4ade diff --git a/src/struphy/feec/boundary_mass.py b/src/struphy/feec/boundary_mass.py index 0fbea0b5f..77064786b 100644 --- a/src/struphy/feec/boundary_mass.py +++ b/src/struphy/feec/boundary_mass.py @@ -398,6 +398,7 @@ def __init__( self._wts_l = self._derham.spline_attributes[self._space_key].quad_grid_wts self._bases_l = self._derham.spline_attributes[self._space_key].quad_grid_bases self._tensor_fem_spaces = self._derham.spline_attributes[self._space_key].tensor_spaces + self._nbasis = self._derham.spline_attributes[self._space_key].nbasis self._V_extraction_op = self._derham.extraction_ops[self._space_key] self._W_extraction_op = self._derham.extraction_ops[self._space_key] @@ -471,11 +472,11 @@ def __init__( # constant skew-symmetric cross-product matrix R_n such that R_n v = n_hat x v R_n_const = xp.zeros((3, 3)) R_n_const[0, 1] = -n_hat[2] - R_n_const[0, 2] = n_hat[1] - R_n_const[1, 0] = n_hat[2] + R_n_const[0, 2] = n_hat[1] + R_n_const[1, 0] = n_hat[2] R_n_const[1, 2] = -n_hat[0] R_n_const[2, 0] = -n_hat[1] - R_n_const[2, 1] = n_hat[0] + R_n_const[2, 1] = n_hat[0] # store R_n per component mu on its own quadrature grid shape surface_R_n_per_mu = [None, None, None] @@ -516,7 +517,11 @@ def codomain(self): def dtype(self): return self._dtype - def _assemble_face(self, face_idx: int, mat): + def _assemble_face( + self, + face_idx: int, + mat: BlockLinearOperator, + ): normal_dir = face_idx % 3 surf_dirs = [d for d in range(3) if d != normal_dir] @@ -533,41 +538,48 @@ def _assemble_face(self, face_idx: int, mat): ends_nu = [int(e) for e in fem_space_nu.coeff_space.ends] pads_nu = fem_space_nu.coeff_space.pads - boundary_index_mu = starts_mu[normal_dir] if face_idx < 3 else ends_mu[normal_dir] - boundary_index_nu = starts_nu[normal_dir] if face_idx < 3 else ends_nu[normal_dir] + boundary_index_mu = 0 if face_idx < 3 else self._nbasis[mu][normal_dir] - 1 + boundary_index_nu = 0 if face_idx < 3 else self._nbasis[nu][normal_dir] - 1 + + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index_mu=}, {starts_mu=}, {ends_mu=}, {pads_mu=}") + logger.debug(f"{normal_dir=}, {face_idx=} {boundary_index_nu=}, {starts_nu=}, {ends_nu=}, {pads_nu=}") mat_fun_mu_nu = self._surface_R_n[face_idx][mu][..., mu, nu] mat_fun_nu_mu = self._surface_R_n[face_idx][nu][..., nu, mu] - self._assembly_kernel( - *self._surface_spans[face_idx][mu], - *fem_space_mu.degree, - *fem_space_nu.degree, - *starts_mu, - *pads_mu, - *self._surface_wts[face_idx][mu], - *self._surface_bases[face_idx][mu], - *self._surface_bases[face_idx][nu], - boundary_index_mu, - normal_dir, - mat_fun_mu_nu, - mat.blocks[mu][nu]._data, - ) + if starts_mu[normal_dir] == boundary_index_mu or ends_mu[normal_dir] == boundary_index_mu: + logger.debug(f"Assembling face {face_idx} for block ({mu},{nu})") + self._assembly_kernel( + *self._surface_spans[face_idx][mu], + *fem_space_mu.degree, + *fem_space_nu.degree, + *starts_mu, + *pads_mu, + *self._surface_wts[face_idx][mu], + *self._surface_bases[face_idx][mu], + *self._surface_bases[face_idx][nu], + boundary_index_mu, + normal_dir, + mat_fun_mu_nu, + mat.blocks[mu][nu]._data, + ) - self._assembly_kernel( - *self._surface_spans[face_idx][nu], - *fem_space_nu.degree, - *fem_space_mu.degree, - *starts_nu, - *pads_nu, - *self._surface_wts[face_idx][nu], - *self._surface_bases[face_idx][nu], - *self._surface_bases[face_idx][mu], - boundary_index_nu, - normal_dir, - mat_fun_nu_mu, - mat.blocks[nu][mu]._data, - ) + if starts_nu[normal_dir] == boundary_index_nu or ends_nu[normal_dir] == boundary_index_nu: + logger.debug(f"Assembling face {face_idx} for block ({nu},{mu})") + self._assembly_kernel( + *self._surface_spans[face_idx][nu], + *fem_space_nu.degree, + *fem_space_mu.degree, + *starts_nu, + *pads_nu, + *self._surface_wts[face_idx][nu], + *self._surface_bases[face_idx][nu], + *self._surface_bases[face_idx][mu], + boundary_index_nu, + normal_dir, + mat_fun_nu_mu, + mat.blocks[nu][mu]._data, + ) def assemble(self, clear: bool = True): """Assembles the H(curl) boundary mass matrix.""" diff --git a/src/struphy/feec/tests/test_boundary_integrals.py b/src/struphy/feec/tests/test_boundary_integrals.py index 23dc26e20..e79294a18 100644 --- a/src/struphy/feec/tests/test_boundary_integrals.py +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -221,14 +221,17 @@ def test_boundary_mass_hollow_cylinder_nonconstant(num_elements, degree, bcs): @pytest.mark.parametrize("num_elements", [[10, 10, 10]]) @pytest.mark.parametrize("degree", [[2, 2, 2]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -@pytest.mark.parametrize("active_faces, u_idx, v_idx, exact", [ - ([True, False, False, False, False, False], 1, 2, 1.0), - ([False, True, False, False, False, False], 2, 0, 1.0), - ([False, False, True, False, False, False], 0, 1, 1.0), - ([False, False, False, True, False, False], 1, 2, -1.0), - ([False, False, False, False, True, False], 2, 0, -1.0), - ([False, False, False, False, False, True], 0, 1, -1.0), -]) +@pytest.mark.parametrize( + "active_faces, u_idx, v_idx, exact", + [ + ([True, False, False, False, False, False], 1, 2, 1.0), + ([False, True, False, False, False, False], 2, 0, 1.0), + ([False, False, True, False, False, False], 0, 1, 1.0), + ([False, False, False, True, False, False], 1, 2, -1.0), + ([False, False, False, False, True, False], 2, 0, -1.0), + ([False, False, False, False, False, True], 0, 1, -1.0), + ], +) def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): comm = MPI.COMM_WORLD @@ -256,7 +259,7 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u v_h = P(v_funs) bnd_ops = BoundaryIntegralOperators(mass_ops, active_faces=active_faces) - numerical = xp.dot(v_h.toarray(), bnd_ops.S1.dot(u_h).toarray()) + numerical = bnd_ops.S1.dot_inner(u_h, v_h) logger.info(f"numerical = {numerical}, exact = {exact}, error = {xp.abs(numerical - exact)}") @@ -266,14 +269,17 @@ def test_boundary_mass_hcurl_per_face(num_elements, degree, bcs, active_faces, u @pytest.mark.parametrize("num_elements", [[10, 10, 10]]) @pytest.mark.parametrize("degree", [[2, 2, 2]]) @pytest.mark.parametrize("bcs", [(("free", "free"), ("free", "free"), ("free", "free"))]) -@pytest.mark.parametrize("active_faces, u_idx, v_idx, exact", [ - ([True, False, False, False, False, False], 1, 2, 12.0), - ([False, True, False, False, False, False], 2, 0, 6.0), - ([False, False, True, False, False, False], 0, 1, 8.0), - ([False, False, False, True, False, False], 1, 2, -12.0), - ([False, False, False, False, True, False], 2, 0, -6.0), - ([False, False, False, False, False, True ], 0, 1, -8.0), -]) +@pytest.mark.parametrize( + "active_faces, u_idx, v_idx, exact", + [ + ([True, False, False, False, False, False], 1, 2, 12.0), + ([False, True, False, False, False, False], 2, 0, 6.0), + ([False, False, True, False, False, False], 0, 1, 8.0), + ([False, False, False, True, False, False], 1, 2, -12.0), + ([False, False, False, False, True, False], 2, 0, -6.0), + ([False, False, False, False, False, True], 0, 1, -8.0), + ], +) def test_boundary_mass_hcurl_cuboid_nontrivial(num_elements, degree, bcs, active_faces, u_idx, v_idx, exact): """ Tests the H(curl) boundary mass operator on a non-unit cuboid [-1,1] x [-1,3] x [0,3] @@ -294,8 +300,10 @@ def make_pulled(domain, idx): lambda x, y, z: xp.ones_like(x) if 1 == idx else xp.zeros_like(x), lambda x, y, z: xp.ones_like(x) if 2 == idx else xp.zeros_like(x), ] + def pulled(*etas): return domain.pull(phys_funs, *etas, kind="1") + return [ lambda *etas, p=pulled: p(*etas)[0], lambda *etas, p=pulled: p(*etas)[1], @@ -347,8 +355,9 @@ def pulled(*etas): [2, 2, 2], (("free", "free"), ("free", "free"), ("free", "free")), [True, False, False, False, False, False], - 1, 2, - 1.0 + 1, + 2, + 1.0, ) test_boundary_mass_hcurl_cuboid_nontrivial( @@ -356,6 +365,7 @@ def pulled(*etas): [1, 2, 3], (("free", "free"), ("free", "free"), ("free", "free")), [True, False, False, False, False, False], - 1, 2, + 1, + 2, 12.0, )