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 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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