From 20cc5b190153ada2a3b8f2be3b071ba2942050dc Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 20 Oct 2023 16:44:11 +0200 Subject: [PATCH] solved the problem, update tests and add tests in feec api --- psydac/api/fem.py | 4 +- psydac/api/grid.py | 6 +- psydac/api/tests/test_api_feec_1d.py | 8 +- psydac/api/tests/test_api_feec_2d.py | 24 +++--- psydac/api/tests/test_api_feec_3d.py | 2 +- psydac/core/bsplines.py | 42 ++++----- psydac/core/bsplines_kernels.py | 37 ++++---- psydac/core/tests/test_bsplines.py | 14 +-- psydac/core/tests/test_bsplines_pyccel.py | 67 +++++++-------- .../tests/test_differentiation_matrices.py | 59 +++++++------ psydac/feec/tests/test_global_projectors.py | 85 +++++++++++-------- psydac/fem/splines.py | 19 +++-- psydac/fem/tensor.py | 8 +- 13 files changed, 193 insertions(+), 182 deletions(-) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 20de7baaf..cf593d140 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -610,7 +610,7 @@ def construct_arguments(self, with_openmp=False): space = spaces[i].spaces[axis] points_i = points[i][axis] periodic = space.periodic - local_span = find_span(space.knots, space.degree, points_i[0, 0], periodic) + local_span = find_span(space.knots, space.degree, points_i[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points_i[0, 0], local_span, nderiv, space.basis) map_basis[i][axis] = map_basis[i][axis].copy() @@ -1214,7 +1214,7 @@ def construct_arguments(self, with_openmp=False): nderiv = self.max_nderiv space = space.spaces[axis] points = points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) + local_span = find_span(space.knots, space.degree, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) map_basis[axis] = map_basis[axis].copy() diff --git a/psydac/api/grid.py b/psydac/api/grid.py index 1771bc131..f04a1e01f 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -206,7 +206,7 @@ def __init__( self, V, nderiv , trial=False, grid=None): for i,Vi in enumerate(V): space = Vi.spaces[axis] points = grid.points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) + local_span = find_span(space.knots, space.degree, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) @@ -230,13 +230,11 @@ def space(self): # TODO have a parallel version of this function, as done for fem def create_collocation_basis( glob_points, space, nderiv=1 ): - #V.C 12/10/23 : Is this function really used somewhere? What is it suppose to do? T = space.knots # knots sequence p = space.degree # spline degree n = space.nbasis # total number of control points grid = space.breaks # breakpoints nc = space.ncells # number of cells in domain (nc=len(grid)-1) - perio = space.periodic # periodicity of the space #------------------------------------------- # GLOBAL GRID @@ -248,7 +246,7 @@ def create_collocation_basis( glob_points, space, nderiv=1 ): # glob_basis = np.zeros( (p+1,nderiv+1,nq) ) # TODO use this for local basis fct glob_basis = np.zeros( (n+p,nderiv+1,nq) ) # n+p for ghosts for iq,xq in enumerate(glob_points): - span = find_span( T, p, xq, perio ) + span = find_span( T, p, xq) glob_spans[iq] = span ders = basis_funs_all_ders( T, p, xq, span, nderiv ) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index 93f69a45d..3949454db 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -415,8 +415,8 @@ def discrete_energies(self, e, b): errB = lambda x1: (push_1d_l2(B, x1, F) - B_ex(t, *F(x1)))**2 * np.sqrt(F.metric_det(x1)) error_l2_E = np.sqrt(derham_h.V1.integral(errE)) error_l2_B = np.sqrt(derham_h.V0.integral(errB)) - print('L2 norm of error on E(t,x) at final time: {:.2e}'.format(error_l2_E)) - print('L2 norm of error on B(t,x) at final time: {:.2e}'.format(error_l2_B)) + print('L2 norm of error on E(t,x) at final time: {:.12e}'.format(error_l2_E)) + print('L2 norm of error on B(t,x) at final time: {:.12e}'.format(error_l2_B)) if diagnostics_interval: @@ -518,8 +518,8 @@ def test_maxwell_1d_periodic_mult(): ) TOL = 1e-6 - ref = dict(error_E = 4.54259221e-04, - error_B = 3.71519383e-04) + ref = dict(error_E = 4.24689338e-04, + error_B = 4.03195792e-04) assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index c0c8f0c0d..2db063d02 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -730,9 +730,9 @@ def test_maxwell_2d_multiplicity(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 4.0298626142e-02, - error_l2_Ey = 4.0298626142e-02, - error_l2_Bz = 3.0892842064e-01) + ref = dict(error_l2_Ex = 4.35005708e-04, + error_l2_Ey = 4.35005708e-04, + error_l2_Bz = 3.76106860e-03) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -746,8 +746,7 @@ def test_maxwell_2d_periodic_multiplicity(): ncells = 30, degree = 3, periodic = True, - Cp = 0.5, #it's seems that with higher multiplicity we need - #to be more carefull with the CFL (with C=5. error is ) + Cp = 0.5, nsteps = 1, tend = None, splitting_order = 2, @@ -759,9 +758,9 @@ def test_maxwell_2d_periodic_multiplicity(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 8.2398553953e-02, - error_l2_Ey = 8.2398553953e-02, - error_l2_Bz = 2.6713526134e-01) + ref = dict(error_l2_Ex = 1.78557685e-04, + error_l2_Ey = 1.78557685e-04, + error_l2_Bz = 1.40582413e-04) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -788,9 +787,9 @@ def test_maxwell_2d_periodic_multiplicity_equal_deg(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 7.43229667e-02, - error_l2_Ey = 7.43229667e-02, - error_l2_Bz = 2.59863979e-01) + ref = dict(error_l2_Ex = 2.50585008e-02, + error_l2_Ey = 2.50585008e-02, + error_l2_Bz = 1.58438290e-02) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -913,9 +912,6 @@ def test_maxwell_2d_dirichlet_par(): # SCRIPT CAPABILITIES #============================================================================== if __name__ == '__main__': - - test_maxwell_2d_multiplicity() - exit() import argparse diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index 1bbf53318..ac7db93c0 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -395,7 +395,7 @@ class CollelaMapping3D(Mapping): T = dt*niter error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=2) - assert abs(error - 0.9020674104318802) < 1e-9 + assert abs(error - 0.24749763720543216) < 1e-9 #============================================================================== # CLEAN UP SYMPY NAMESPACE diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 2fa85aa05..13fc22026 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -57,7 +57,7 @@ #============================================================================== -def find_span(knots, degree, x, periodic, multiplicity = 1): +def find_span(knots, degree, x): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -74,15 +74,7 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): Polynomial degree of B-splines. x : float - Location of interest. - - periodic : bool - Indicating if the domain is periodic. - - multiplicity : int - Multplicity in the knot sequence. One assume that all the knots have - multiplicity multiplicity, and in the none periodic case that the - boundary knots have multiplicity p+1 (has created by make_knots) + Location of interest.. Returns ------- @@ -91,12 +83,10 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): """ x = float(x) knots = np.ascontiguousarray(knots, dtype=float) - multiplicity = int(multiplicity) - periodic = bool(periodic) - return find_span_p(knots, degree, x, periodic, multiplicity = multiplicity) + return find_span_p(knots, degree, x) #============================================================================== -def find_spans(knots, degree, x, periodic, out=None): +def find_spans(knots, degree, x, out=None): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -131,7 +121,7 @@ def find_spans(knots, degree, x, periodic, out=None): else: assert out.shape == x.shape and out.dtype == np.dtype('int') - find_spans_p(knots, degree, x, periodic, out) + find_spans_p(knots, degree, x, out) return out #============================================================================== @@ -336,6 +326,10 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, out : array, optional If provided, the result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Returns ------- @@ -352,7 +346,7 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, if out is None: nb = len(knots) - degree - 1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity out = np.zeros((xgrid.shape[0], nb), dtype=float) else: @@ -386,6 +380,10 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli xgrid : array_like Grid points. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. check_boundary : bool, default=True If true and ``periodic``, will check the boundaries of ``xgrid``. @@ -435,12 +433,12 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli if out is None: if periodic: - out = np.zeros((len(xgrid), len(knots) - 2 * degree - 1), dtype=float) + out = np.zeros((len(xgrid), len(knots) - 2 * degree - 2 + multiplicity), dtype=float) else: out = np.zeros((len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1), dtype=float) else: if periodic: - assert out.shape == (len(xgrid), len(knots) - 2 * degree - 1) + assert out.shape == (len(xgrid), len(knots) - 2 * degree - 2 + multiplicity) else: assert out.shape == (len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1) assert out.dtype == np.dtype('float') @@ -501,6 +499,10 @@ def greville(knots, degree, periodic, out=None, multiplicity=1): out : array, optional If provided, the result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Returns ------- @@ -510,7 +512,7 @@ def greville(knots, degree, periodic, out=None, multiplicity=1): """ knots = np.ascontiguousarray(knots, dtype=float) if out is None: - n = len(knots) - 2 * degree - 1 if periodic else len(knots) - degree - 1 + n = len(knots) - 2 * degree - 2 + multiplicity if periodic else len(knots) - degree - 1 out = np.zeros(n) multiplicity = int(multiplicity) greville_p(knots, degree, periodic, out, multiplicity) @@ -1011,7 +1013,7 @@ def alpha_function(i, k, t, n, p, knots): mat = np.zeros((n+1,n)) - left = find_span( knots, p, t, False ) + left = find_span( knots, p, t ) # ... j = 0 diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index d5c60be6a..b41d3305f 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -5,7 +5,7 @@ import numpy as np # ============================================================================= -def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multiplicity : int = 1): +def find_span_p(knots: 'float[:]', degree: int, x: float): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -24,14 +24,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi x : float Location of interest. - periodic : bool - Indicating if the domain is periodic. - - multiplicity : int - Multplicity in the knot sequence. One assume that all the knots have - multiplicity multiplicity, and in the none periodic case that the - boundary knots have multiplicity p+1 (has created by make_knots) - Returns ------- span : int @@ -43,9 +35,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi Springer-Verlag Berlin Heidelberg GmbH, 1997. """ # Knot index at left/right boundary - #In the periodic case the p+1 knots is the last knot equal to the left - #boundary, while the first knots being the right boundary is multiplicity+p - #before the last knot of the sequence (see make knots) low = degree high = len(knots)-1-degree @@ -66,7 +55,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi # ============================================================================= -def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', periodic : bool, multiplicity : int = 1): +def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]'): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -92,7 +81,7 @@ def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', p n = x.shape[0] for i in range(n): - out[i] = find_span_p(knots, degree, x[i], periodic, multiplicity=multiplicity) + out[i] = find_span_p(knots, degree, x[i]) # ============================================================================= @@ -422,11 +411,15 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali The result will be inserted into this array. It should be of the appropriate shape and dtype. + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. + """ # Number of basis functions (in periodic case remove degree repeated elements) nb = len(knots)-degree-1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity # Number of evaluation points nx = len(xgrid) @@ -434,7 +427,7 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans, periodic, multiplicity = multiplicity) + find_spans_p(knots, degree, xgrid, spans) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values @@ -500,6 +493,10 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma out : array The result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Notes ----- @@ -509,7 +506,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma """ nb = len(knots) - degree - 1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity # Number of evaluation points nx = len(xgrid) @@ -718,10 +715,14 @@ def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]', m out : array The result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. """ T = knots p = degree - n = len(T)-2*p-1 if periodic else len(T)-p-1 + n = len(T)-2*p-2 + multiplicity if periodic else len(T)-p-1 # Compute greville abscissas as average of p consecutive knot values if p == multiplicity-1: diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index c539bbd9a..4b3f0b987 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -32,9 +32,9 @@ def test_find_span( lims, nc, p, eps=1e-12 ): knots = np.r_[ [grid[0]]*p, grid, [grid[-1]]*p ] for i,xi in enumerate( grid ): - assert find_span( knots, p, x=xi-eps, periodic=False ) == p + max( 0, i-1 ) - assert find_span( knots, p, x=xi, periodic=False ) == p + min( i, nc-1 ) - assert find_span( knots, p, x=xi+eps, periodic=False ) == p + min( i, nc-1 ) + assert find_span( knots, p, x=xi-eps ) == p + max( 0, i-1 ) + assert find_span( knots, p, x=xi ) == p + min( i, nc-1 ) + assert find_span( knots, p, x=xi+eps ) == p + min( i, nc-1 ) #============================================================================== @pytest.mark.parametrize( 'lims', ([0,1], [-2,3]) ) @@ -48,7 +48,7 @@ def test_basis_funs( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) basis = basis_funs( knots, p, x, span ) assert len( basis ) == p+1 assert np.all( basis >= 0 ) @@ -66,7 +66,7 @@ def test_basis_funs_1st_der( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) ders = basis_funs_1st_der( knots, p, x, span ) assert len( ders ) == p+1 assert abs( sum( ders ) ) < tol @@ -86,7 +86,7 @@ def test_basis_funs_all_ders( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) ders = basis_funs_all_ders( knots, p, x, span, n ) # Test output array @@ -195,7 +195,7 @@ def test_cell_index(i_grid, expected): yy = np.zeros( (len(xx), nb) ) zz = np.zeros( (len(xx), nb) ) for i,x in enumerate( xx ): - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) yy[i,span-p:span+1] = basis_funs ( knots, p, x, span ) zz[i,span-p:span+1] = basis_funs_1st_der( knots, p, x, span ) diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index aeb9f68cd..f243f39ae 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -21,13 +21,11 @@ # "True" Functions ############################################################################### -def find_span_true( knots, degree, x, periodic, multiplicity=1 ): +def find_span_true( knots, degree, x): # Knot index at left/right boundary low = degree - if periodic : - high = len(knots)-multiplicity-degree - else : - high = len(knots)-1-degree + high = len(knots)-1-degree + # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low if x >= knots[high]: return high-1 @@ -211,7 +209,7 @@ def collocation_matrix_true(knots, degree, periodic, normalization, xgrid): # Fill in non-zero matrix values for i,x in enumerate( xgrid ): - span = find_span_true( knots, degree, x, periodic ) + span = find_span( knots, degree, x ) basis = basis_funs_true( knots, degree, x, span ) mat[i,js(span)] = normalize(basis, span) @@ -377,27 +375,29 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): if periodic: assert len(breaks) > degree - - p = degree - if periodic : - T = np.zeros(multiplicity * len(breaks[1:]) + 1 + 2 * degree) - else : - T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) - + + T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) + + ncells = len(breaks) - 1 + + for i in range(0, ncells+1): + T[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] + + len_out = len(T) + if periodic: period = breaks[-1]-breaks[0] - T[p+1:-p] = np.repeat(breaks[1:], multiplicity) - T[0:p+1] = [xi-period for xi in T[-2*p-1:-p ]] - T[-p:] = [xi+period for xi in T[ p+1:2*p+1]] + T[: degree + 1 - multiplicity] = T[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period + T[len_out - (degree + 1 - multiplicity) :] = T[degree + 1:2*(degree + 1)- multiplicity] + period + + + # else: - T[p+1:-p-1] = np.repeat(breaks[1:-1], multiplicity) - T[p] = breaks[ 0] - T[-p-1] = breaks[-1] - T[0:p+1] = breaks[ 0] - T[-p-1:] = breaks[-1] + T[0:degree + 1 - multiplicity] = breaks[0] + T[len_out - degree - 1 + multiplicity:] = breaks[-1] return T @@ -406,10 +406,10 @@ def elevate_knots_true(knots, degree, periodic, multiplicity=1, tol=1e-15): knots = np.array(knots) if periodic: - [T, p] = knots, degree - period = T[-1-p] - T[p] - left = [T[-1-p-(p+1)] - period] - right = [T[ p+(p+1)] + period] + T, p = knots, degree + period = T[len(knots) -1 - p] - T[p] + left = [T[len(knots) -2 - 2 * p + multiplicity-1] - period] + right = [T[2 * p + 2 - multiplicity] + period] else: left = [knots[0],*knots[:degree+1]] right = [knots[-1],*knots[-degree-1:]] @@ -461,7 +461,7 @@ def basis_ders_on_quad_grid_true(knots, degree, quad_grid, nders, normalization) for ie in range(ne): xx = quad_grid[ie, :] for iq, xq in enumerate(xx): - span = find_span_true(knots, degree, xq, False) + span = find_span_true(knots, degree, xq) ders = basis_funs_all_ders_true(knots, degree, xq, span, nders) if normalization == 'M': ders *= scaling[None, span-degree:span+1] @@ -499,11 +499,10 @@ def basis_integrals_true(knots, degree): (np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0]), 2), (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) -@pytest.mark.parametrize('periodic', [True, False]) -@pytest.mark.parametrize('multiplicity', [1, 2, 3]) -def test_find_span(knots, degree, x, periodic, multiplicity): - expected = find_span_true(knots, degree, x, periodic, multiplicity=multiplicity) - out = find_span(knots, degree, x, periodic, multiplicity=multiplicity) + +def test_find_span(knots, degree, x): + expected = find_span_true(knots, degree, x) + out = find_span(knots, degree, x) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) @@ -521,7 +520,7 @@ def test_find_span(knots, degree, x, periodic, multiplicity): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs(knots, degree, x): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_true(knots, degree, x, span) out = basis_funs(knots, degree, x, span) @@ -541,7 +540,7 @@ def test_basis_funs(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs_1st_der(knots, degree, x): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_1st_der_true(knots, degree, x, span) out = basis_funs_1st_der(knots, degree, x, span) @@ -563,7 +562,7 @@ def test_basis_funs_1st_der(knots, degree, x): @pytest.mark.parametrize('n', (2, 3, 4, 5)) @pytest.mark.parametrize('normalization', ('B', 'M')) def test_basis_funs_all_ders(knots, degree, x, n, normalization): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_all_ders_true(knots, degree, x, span, n, normalization) out = basis_funs_all_ders(knots, degree, x, span, n, normalization) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 3feccd55d..35e017421 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -461,21 +461,18 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): @pytest.mark.parametrize('seed', [1,3]) @pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) - def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): if any([ncells[d] <= degree[d] and periodic[d] for d in range(3)]): return - #change mulitplicity if higher than degree to avoid problems (case p1e-15)[0] - if len(indices)>0: - multiplicity = np.diff(indices).max(initial=1) - else: - multiplicity = max(1,len(knots[degree+1:-degree-1])) + if multiplicity is None: #V.C 20/10/23 : why computing the multiplicity again? + if len(indices)>0: + multiplicity = np.diff(indices).max(initial=1) + else: + multiplicity = max(1,len(knots[degree+1:-degree-1])) + if parent_multiplicity is None: parent_multiplicity = multiplicity @@ -107,7 +109,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul # Number of basis function in space (= cardinality) if periodic: - nbasis = len(knots) - 2*degree - 1 + nbasis = len(knots) - 2*degree - 2 + multiplicity else: defect = 0 if dirichlet[0]: defect += 1 @@ -134,7 +136,6 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul 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._scaling_array = scaling_array - self._parent_multiplicity = parent_multiplicity self._histopolation_grid = unroll_edges(self.domain, self.ext_greville) @@ -181,7 +182,7 @@ def init_interpolation( self, dtype=float ): xgrid = self.greville, multiplicity = self.multiplicity ) - + if self.periodic: # Convert to CSC format and compute sparse LU decomposition self._interpolator = SparseSolver( csc_matrix( imat ) ) @@ -290,7 +291,7 @@ def eval_field(self, field, *eta , weights=None): eta = eta[0] - span = find_span( self.knots, self.degree, eta, self.periodic) + span = find_span( self.knots, self.degree, eta) basis_array = basis_funs( self.knots, self.degree, eta, span) index = slice(span-self.degree, span + 1) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 6a9b0b86b..4ea6a9b74 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -193,9 +193,7 @@ def eval_field( self, field, *eta, weights=None): for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): knots = space.knots degree = space.degree - multiplicity = space.multiplicity - periodic = space.periodic - span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) + span = find_span( knots, degree, x ) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# @@ -622,9 +620,7 @@ def eval_field_gradient( self, field, *eta , weights=None): knots = space.knots degree = space.degree - multiplicity = space.multiplicity - periodic = space.periodic - span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) + span = find_span( knots, degree, x ) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------#