From cd576f7daa30bae68e17cf8542ff21135d73f9fa Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 22 Jan 2026 15:06:47 +0100 Subject: [PATCH 1/5] remove --force-reinstall --no-cache-dir from psydac install --- .github/workflows/test-struphy.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test-struphy.yml b/.github/workflows/test-struphy.yml index 10608aa09..dc729f183 100644 --- a/.github/workflows/test-struphy.yml +++ b/.github/workflows/test-struphy.yml @@ -78,7 +78,7 @@ jobs: psydac-accelerate --cleanup --yes pip show feectools pip uninstall feectools -y - python -m pip install ${{ env.PSYDAC_DIR }} --force-reinstall --no-cache-dir + python -m pip install ${{ env.PSYDAC_DIR }} echo "feectools location after installing feectools from the current branch (should be the same as in the previous step)" pip show feectools From 40dd9df6ddbf204ba026e7f62979bf85cc12da5c Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 22 Jan 2026 15:12:41 +0100 Subject: [PATCH 2/5] adapt greville points for dirichlet bcs --- feectools/fem/splines.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/feectools/fem/splines.py b/feectools/fem/splines.py index 574e4cbb1..75046587b 100644 --- a/feectools/fem/splines.py +++ b/feectools/fem/splines.py @@ -111,6 +111,13 @@ def __init__(self, degree, knots=None, grid=None, multiplicity=None, parent_mult if dirichlet[1]: defect += 1 nbasis = len(knots) - degree - 1 - defect + # greville points + grev_tmp = greville(knots, degree, periodic, multiplicity = multiplicity) + if dirichlet[0]: + grev_tmp = grev_tmp[1:] + if dirichlet[1]: + grev_tmp = grev_tmp[:-1] + # Coefficients to convert B-splines to M-splines (if needed) if basis == 'M': scaling_array = 1 / basis_integrals(knots, degree) @@ -128,7 +135,7 @@ def __init__(self, degree, knots=None, grid=None, multiplicity=None, parent_mult self._nbasis = nbasis self._breaks = grid self._ncells = len(grid) - 1 - self._greville = greville(knots, degree, periodic, multiplicity = multiplicity) + self._greville = grev_tmp self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) self._scaling_array = scaling_array self._parent_multiplicity = parent_multiplicity @@ -178,6 +185,7 @@ def init_interpolation( self, dtype=float ): periodic = self.periodic, normalization = self.basis, xgrid = self.greville, + dirichlet = self.dirichlet, multiplicity = self.multiplicity ) From df13013c0be1d9a0705805a878756d83a3e90604 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Jan 2026 07:08:51 +0100 Subject: [PATCH 3/5] pass dirichlet to collocation matrix --- feectools/core/bsplines.py | 18 +++--- feectools/core/bsplines_kernels.py | 59 ++++++++++++++++--- feectools/feec/global_geometric_projectors.py | 6 ++ feectools/fem/splines.py | 1 + 4 files changed, 68 insertions(+), 16 deletions(-) diff --git a/feectools/core/bsplines.py b/feectools/core/bsplines.py index 900178871..00111a872 100644 --- a/feectools/core/bsplines.py +++ b/feectools/core/bsplines.py @@ -300,7 +300,7 @@ def basis_funs_all_ders(knots, degree, x, span, n, normalization='B', out=None): return out #============================================================================== -def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, multiplicity = 1): +def collocation_matrix(knots, degree, periodic, normalization, xgrid, dirichlet, out=None, multiplicity = 1): """Computes the collocation matrix If called with normalization='M', this uses M-splines instead of B-splines. @@ -321,6 +321,9 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, xgrid : array_like Evaluation points. + + dirichlet: tuple[bool] + Whether to have hom. Dirichlet conditions left and/or right. out : array, optional If provided, the result will be inserted into this array. @@ -345,10 +348,11 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, knots = np.ascontiguousarray(knots, dtype=float) xgrid = np.ascontiguousarray(xgrid, dtype=float) + nb = len(xgrid) if out is None: - nb = len(knots) - degree - 1 - if periodic: - nb -= degree + 1 - multiplicity + # nb = len(knots) - degree - 1 + # if periodic: + # nb -= degree + 1 - multiplicity out = np.zeros((xgrid.shape[0], nb), dtype=float) else: @@ -357,11 +361,11 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, bool_normalization = normalization == "M" multiplicity = int(multiplicity) - collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, out, multiplicity=multiplicity) + collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, dirichlet, out, multiplicity=multiplicity) return out #============================================================================== -def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multiplicity=1, check_boundary=True, out=None): +def histopolation_matrix(knots, degree, periodic, normalization, xgrid, dirichlet, multiplicity=1, check_boundary=True, out=None): """Computes the histopolation matrix. If called with normalization='M', this uses M-splines instead of B-splines. @@ -445,7 +449,7 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli assert out.shape == (len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1) assert out.dtype == np.dtype('float') multiplicity = int(multiplicity) - histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, check_boundary, elevated_knots, out, multiplicity = multiplicity) + histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, dirichlet, check_boundary, elevated_knots, out, multiplicity = multiplicity) return out #============================================================================== diff --git a/feectools/core/bsplines_kernels.py b/feectools/core/bsplines_kernels.py index c8ba08e92..81b38abbb 100644 --- a/feectools/core/bsplines_kernels.py +++ b/feectools/core/bsplines_kernels.py @@ -512,8 +512,14 @@ def basis_integrals_p(knots: 'float[:]', degree: int, out: 'float[:]'): # ============================================================================= -def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - out: 'float[:,:]', multiplicity : int = 1): +def collocation_matrix_p(knots: 'float[:]', + degree: int, + periodic: bool, + normalization: bool, + xgrid: 'float[:]', + dirichlet: tuple[bool], + out: 'float[:,:]', + multiplicity : int = 1): """ Compute the collocation matrix :math:`C_ij = B_j(x_i)`, which contains the values of each B-spline basis function :math:`B_j` at all locations :math:`x_i`. @@ -547,9 +553,10 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali """ # Number of basis functions (in periodic case remove degree repeated elements) - nb = len(knots)-degree-1 - if periodic: - nb -= degree + 1 - multiplicity + nb = len(xgrid) + # nb = len(knots)-degree-1 + # if periodic: + # nb -= degree + 1 - multiplicity # Number of evaluation points nx = len(xgrid) @@ -570,7 +577,19 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali out[i, actual_j] = basis[i, j] else: for i in range(nx): - out[i, spans[i] - degree:spans[i] + 1] = basis[i, :] + print(f"\n{i = }") + print("not normalization") + print(f"{dirichlet = }") + print(f"{spans = }") + print(f"{degree = }") + print(f"{spans[i] - dirichlet[0] - degree = }") + print(f"{spans[i] - dirichlet[0] + 1 = }") + print(f"{basis[i, :] = }") + if dirichlet[1] and i==nx - 1: + out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0]] = basis[i, :-1] + else: + out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0] + 1] = basis[i, :] + print(f"{out = }") else: integrals = np.zeros(knots.shape[0] - degree - 1) basis_integrals_p(knots, degree, integrals) @@ -584,7 +603,20 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali else: for i in range(nx): local_scaling = scaling[spans[i] - degree:spans[i] + 1] - out[i, spans[i] - degree:spans[i] + 1] = basis[i, :] * local_scaling[:] + print(f"\n{i = }") + print("else") + print(f"{dirichlet = }") + print(f"{spans = }") + print(f"{degree = }") + print(f"{spans[i] - dirichlet[0] - degree = }") + print(f"{spans[i] - dirichlet[0] + 1 = }") + print(f"{basis[i, :] = }") + print(f"{local_scaling[:] = }") + # if dirichlet[1] and i==nx - 1: + # out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0]] = basis[i, :-1] * local_scaling[:] + # else: + out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0] + 1] = basis[i, :] * local_scaling[:] + print(f"{out = }") # Mitigate round-off errors for x in range(nx): @@ -594,8 +626,16 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali # ============================================================================= -def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - check_boundary: bool, elevated_knots: 'float[:]', out: 'float[:,:]', multiplicity : int = 1): +def histopolation_matrix_p(knots: 'float[:]', + degree: int, + periodic: bool, + normalization: bool, + xgrid: 'float[:]', + dirichlet: tuple[bool], + check_boundary: bool, + elevated_knots: 'float[:]', + out: 'float[:,:]', + multiplicity : int = 1): """Computes the histopolation matrix. If called with normalization=True, this uses M-splines instead of B-splines. @@ -684,6 +724,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma False, False, xgrid_new[:actual_len], + dirichlet, colloc, multiplicity = multiplicity) diff --git a/feectools/feec/global_geometric_projectors.py b/feectools/feec/global_geometric_projectors.py index dac6a4312..7449a3268 100644 --- a/feectools/feec/global_geometric_projectors.py +++ b/feectools/feec/global_geometric_projectors.py @@ -175,6 +175,12 @@ def __init__(self, space, nquads = None): M._data[row_i_loc + m*p, (col_i + p - row_i)%V.imat.shape[1]] = V.imat[row_i, col_i] # check if stencil matrix was built correctly + print(f"\n{M.toarray()[s:e + 1] = }") + print(f"\n{V.imat[s:e + 1] = }") + print(f"{V = }") + print(f"\n{M.toarray().shape = }") + print(f"\n{V.imat.shape = }") + assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1]) # TODO Fix toarray() for multiplicity m > 1 matrixcells += [M.copy()] diff --git a/feectools/fem/splines.py b/feectools/fem/splines.py index 75046587b..ff1015680 100644 --- a/feectools/fem/splines.py +++ b/feectools/fem/splines.py @@ -221,6 +221,7 @@ def init_histopolation( self, dtype=float): periodic = self.periodic, normalization = self.basis, xgrid = self.ext_greville, + dirichlet = self.dirichlet, multiplicity = self._multiplicity ) From 5da77fcf7089ba92648657eb840fa9dee474ee4f Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Jan 2026 13:54:08 +0100 Subject: [PATCH 4/5] adapt bsplines to smaller spaces for dirichlet bcs --- feectools/core/bsplines.py | 8 ++-- feectools/core/bsplines_kernels.py | 69 ++++++++++++++++++------------ 2 files changed, 45 insertions(+), 32 deletions(-) diff --git a/feectools/core/bsplines.py b/feectools/core/bsplines.py index 00111a872..93f9da008 100644 --- a/feectools/core/bsplines.py +++ b/feectools/core/bsplines.py @@ -439,14 +439,14 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, dirichle if out is None: if periodic: - out = np.zeros((len(xgrid), len(knots) - 2 * degree - 2 + multiplicity), dtype=float) + out = np.zeros((len(xgrid), len(xgrid)), dtype=float) else: - out = np.zeros((len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1), dtype=float) + out = np.zeros((len(xgrid) - 1, len(xgrid) - 1), dtype=float) else: if periodic: - assert out.shape == (len(xgrid), len(knots) - 2 * degree - 2 + multiplicity) + assert out.shape == (len(xgrid), len(xgrid)) else: - assert out.shape == (len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1) + assert out.shape == (len(xgrid) - 1, len(xgrid) - 1) assert out.dtype == np.dtype('float') multiplicity = int(multiplicity) histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, dirichlet, check_boundary, elevated_knots, out, multiplicity = multiplicity) diff --git a/feectools/core/bsplines_kernels.py b/feectools/core/bsplines_kernels.py index 81b38abbb..9c55f115f 100644 --- a/feectools/core/bsplines_kernels.py +++ b/feectools/core/bsplines_kernels.py @@ -566,6 +566,12 @@ def collocation_matrix_p(knots: 'float[:]', find_spans_p(knots, degree, xgrid, spans) basis_funs_array_p(knots, degree, xgrid, spans, basis) + # adapt indexing for hom. dirichlet + if spans[0] == degree: + adapt_col_index = False + else: + adapt_col_index = dirichlet[0] + # Fill in non-zero matrix values # Rescaling of B-splines, to get M-splines if needed @@ -576,20 +582,23 @@ def collocation_matrix_p(knots: 'float[:]', actual_j = (spans[i] - degree + j) % nb out[i, actual_j] = basis[i, j] else: + k = 0 for i in range(nx): - print(f"\n{i = }") - print("not normalization") - print(f"{dirichlet = }") - print(f"{spans = }") - print(f"{degree = }") - print(f"{spans[i] - dirichlet[0] - degree = }") - print(f"{spans[i] - dirichlet[0] + 1 = }") - print(f"{basis[i, :] = }") - if dirichlet[1] and i==nx - 1: - out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0]] = basis[i, :-1] + # print(f"\n{i = }") + # print("not normalization") + # print(f"{dirichlet = }") + # print(f"{spans = }") + # print(f"{degree = }") + # print(f"{adapt_col_index = }") + # print(f"{spans[i] - adapt_col_index - degree = }") + # print(f"{spans[i] - adapt_col_index + 1 = }") + # print(f"{basis[i, :] = }") + if dirichlet[1] and i>=nx - degree: + k += 1 + out[i, spans[i] - adapt_col_index - degree:spans[i] - adapt_col_index + 1 - k] = basis[i, :-k] else: - out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0] + 1] = basis[i, :] - print(f"{out = }") + out[i, spans[i] - adapt_col_index - degree:spans[i] - adapt_col_index + 1] = basis[i, :] + # print(f"{out = }") else: integrals = np.zeros(knots.shape[0] - degree - 1) basis_integrals_p(knots, degree, integrals) @@ -601,22 +610,25 @@ def collocation_matrix_p(knots: 'float[:]', out[i, actual_j] = basis[i, j] * scaling[spans[i] - degree + j] else: + k = 0 for i in range(nx): local_scaling = scaling[spans[i] - degree:spans[i] + 1] - print(f"\n{i = }") - print("else") - print(f"{dirichlet = }") - print(f"{spans = }") - print(f"{degree = }") - print(f"{spans[i] - dirichlet[0] - degree = }") - print(f"{spans[i] - dirichlet[0] + 1 = }") - print(f"{basis[i, :] = }") - print(f"{local_scaling[:] = }") - # if dirichlet[1] and i==nx - 1: - # out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0]] = basis[i, :-1] * local_scaling[:] - # else: - out[i, spans[i] - dirichlet[0] - degree:spans[i] - dirichlet[0] + 1] = basis[i, :] * local_scaling[:] - print(f"{out = }") + # print(f"\n{i = }") + # print("else") + # print(f"{dirichlet = }") + # print(f"{spans = }") + # print(f"{degree = }") + # print(f"{adapt_col_index = }") + # print(f"{spans[i] - adapt_col_index - degree = }") + # print(f"{spans[i] - adapt_col_index + 1 = }") + # print(f"{basis[i, :] = }") + # print(f"{local_scaling[:] = }") + if dirichlet[1] and i>=nx - degree: + k += 1 + out[i, spans[i] - adapt_col_index - degree:spans[i] - adapt_col_index + 1 - k] = basis[i, :-k] * local_scaling[:-k] + else: + out[i, spans[i] - adapt_col_index - degree:spans[i] - adapt_col_index + 1] = basis[i, :] * local_scaling[:] + # print(f"{out = }") # Mitigate round-off errors for x in range(nx): @@ -717,8 +729,9 @@ def histopolation_matrix_p(knots: 'float[:]', # NOTES: # . cannot use M-splines in analytical formula for histopolation matrix # . always use non-periodic splines to avoid circulant matrix structure - nb_elevated = len(elevated_knots) - (degree + 1) - 1 - colloc = np.zeros((actual_len, nb_elevated)) + # nb_elevated = len(elevated_knots) - (degree + 1) - 1 + # colloc = np.zeros((actual_len, nb_elevated)) + colloc = np.zeros((actual_len, actual_len)) collocation_matrix_p(elevated_knots, degree + 1, False, From 2060b077976a567819b27c1c1aa9d37f1f78e431 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 23 Jan 2026 13:56:56 +0100 Subject: [PATCH 5/5] adapt external greville points --- feectools/fem/splines.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/feectools/fem/splines.py b/feectools/fem/splines.py index ff1015680..76b06e21e 100644 --- a/feectools/fem/splines.py +++ b/feectools/fem/splines.py @@ -117,6 +117,12 @@ def __init__(self, degree, knots=None, grid=None, multiplicity=None, parent_mult grev_tmp = grev_tmp[1:] if dirichlet[1]: grev_tmp = grev_tmp[:-1] + + ext_grev_tmp = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) + if dirichlet[0]: + ext_grev_tmp = ext_grev_tmp[1:] + if dirichlet[1]: + ext_grev_tmp = ext_grev_tmp[:-1] # Coefficients to convert B-splines to M-splines (if needed) if basis == 'M': @@ -136,7 +142,7 @@ def __init__(self, degree, knots=None, grid=None, multiplicity=None, parent_mult self._breaks = grid self._ncells = len(grid) - 1 self._greville = grev_tmp - self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) + self._ext_greville = ext_grev_tmp self._scaling_array = scaling_array self._parent_multiplicity = parent_multiplicity self._histopolation_grid = unroll_edges(self.domain, self.ext_greville)