Refactor derham in struphy

.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
 

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
 
 #==============================================================================

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)
 

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()]

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
         )