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 diff --git a/feectools/core/bsplines.py b/feectools/core/bsplines.py index 900178871..93f9da008 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. @@ -435,17 +439,17 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli 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, 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..9c55f115f 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) @@ -559,6 +566,12 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali 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 @@ -569,8 +582,23 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali actual_j = (spans[i] - degree + j) % nb out[i, actual_j] = basis[i, j] else: + k = 0 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"{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] - 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) @@ -582,9 +610,25 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali 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] - 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"{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): @@ -594,8 +638,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. @@ -677,13 +729,15 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma # 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, 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 574e4cbb1..76b06e21e 100644 --- a/feectools/fem/splines.py +++ b/feectools/fem/splines.py @@ -111,6 +111,19 @@ 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] + + 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': scaling_array = 1 / basis_integrals(knots, degree) @@ -128,8 +141,8 @@ 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._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) + self._greville = grev_tmp + 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) @@ -178,6 +191,7 @@ def init_interpolation( self, dtype=float ): periodic = self.periodic, normalization = self.basis, xgrid = self.greville, + dirichlet = self.dirichlet, multiplicity = self.multiplicity ) @@ -213,6 +227,7 @@ def init_histopolation( self, dtype=float): periodic = self.periodic, normalization = self.basis, xgrid = self.ext_greville, + dirichlet = self.dirichlet, multiplicity = self._multiplicity )