diff --git a/src/struphy/feec/boundary_integrals.py b/src/struphy/feec/boundary_integrals.py new file mode 100644 index 000000000..64c82a0c4 --- /dev/null +++ b/src/struphy/feec/boundary_integrals.py @@ -0,0 +1,539 @@ +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 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/mass_kernels.py b/src/struphy/feec/mass_kernels.py index 7b4f09720..796f6096d 100644 --- a/src/struphy/feec/mass_kernels.py +++ b/src/struphy/feec/mass_kernels.py @@ -766,3 +766,193 @@ 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 + + +def surface_kernel_3d_mat_h1( + 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 + + +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 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..568d71901 --- /dev/null +++ b/src/struphy/feec/tests/test_boundary_integrals.py @@ -0,0 +1,170 @@ +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_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 + +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_mass_unit_cube_constant(num_elements, degree, bcs): + """ + Tests the boundary mass operator for alpha = 1 on the unit cube. + """ + 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_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.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]]) +@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. + """ + 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_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.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_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.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_ops = BoundaryOperators(mass_ops) + v = bnd_ops.S0.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) + + 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