diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index b94013cac..260bb0381 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -46,7 +46,7 @@ def update_plot(ax, t, sol_ex, sol_num): #============================================================================== def run_maxwell_1d(*, L, eps, ncells, degree, periodic, Cp, nsteps, tend, splitting_order, plot_interval, diagnostics_interval, - bc_mode, tol, verbose): + bc_mode, tol, verbose, mult=1): import numpy as np import matplotlib.pyplot as plt @@ -132,7 +132,7 @@ class CollelaMapping1D(Mapping): # Discrete physical domain and discrete DeRham sequence domain_h = discretize(domain, ncells=[ncells], periodic=[periodic], comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=[degree]) + derham_h = discretize(derham, domain_h, degree=[degree], multiplicity=[mult]) # Discrete bilinear forms a0_h = discretize(a0, domain_h, (derham_h.V0, derham_h.V0), backend=PSYDAC_BACKEND_PYTHON) @@ -496,6 +496,33 @@ def test_maxwell_1d_periodic(): assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL + +def test_maxwell_1d_periodic_mult(): + + namespace = run_maxwell_1d( + L = 1.0, + eps = 0.5, + ncells = 30, + degree = 3, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + bc_mode = None, + verbose = False, + mult = 2 + ) + + TOL = 1e-6 + 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 def test_maxwell_1d_dirichlet_strong(): diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 1da94f125..5c27b678e 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -185,7 +185,7 @@ def update_plot(fig, t, x, y, field_h, field_ex): def run_maxwell_2d_TE(*, use_spline_mapping, eps, ncells, degree, periodic, Cp, nsteps, tend, - splitting_order, plot_interval, diagnostics_interval, tol, verbose): + splitting_order, plot_interval, diagnostics_interval, tol, verbose, mult=1): import os @@ -290,7 +290,7 @@ class CollelaMapping2D(Mapping): #-------------------------------------------------------------------------- if use_spline_mapping: domain_h = discretize(domain, filename=filename, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h) + derham_h = discretize(derham, domain_h, multiplicity = [mult,mult]) periodic_list = mapping.get_callable_mapping().space.periodic degree_list = mapping.get_callable_mapping().space.degree @@ -311,7 +311,7 @@ class CollelaMapping2D(Mapping): else: # Discrete physical domain and discrete DeRham sequence domain_h = discretize(domain, ncells=[ncells, ncells], periodic=[periodic, periodic], comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=[degree, degree]) + derham_h = discretize(derham, domain_h, degree=[degree, degree], multiplicity = [mult,mult]) # Discrete bilinear forms a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) @@ -710,6 +710,91 @@ def test_maxwell_2d_periodic(): assert abs(namespace['error_Ey'] - ref['error_Ey']) / ref['error_Ey'] <= TOL assert abs(namespace['error_Bz'] - ref['error_Bz']) / ref['error_Bz'] <= TOL +def test_maxwell_2d_multiplicity(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 10, + degree = 5, + periodic = False, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult = 2 + ) + + TOL = 1e-5 + ref = dict(error_l2_Ex = 4.350041934920621e-04, + error_l2_Ey = 4.350041934920621e-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 + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL + +def test_maxwell_2d_periodic_multiplicity(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 30, + degree = 3, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult =2 + ) + + TOL = 1e-6 + 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 + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL + + +def test_maxwell_2d_periodic_multiplicity_equal_deg(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 10, + degree = 2, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult =2 + ) + + TOL = 1e-6 + 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 + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL + def test_maxwell_2d_dirichlet(): diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index af65969e4..7de5f0e74 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -86,7 +86,7 @@ def evaluation_all_times(fields, x, y, z): return ak_value #================================================================================== -def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter): +def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=1): #------------------------------------------------------------------------------ # Symbolic objects: SymPDE @@ -113,7 +113,7 @@ def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, pe #------------------------------------------------------------------------------ domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=degree) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity = [mult,mult,mult]) a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) a2_h = discretize(a2, domain_h, (derham_h.V2, derham_h.V2), backend=PSYDAC_BACKEND_GPYCCEL) @@ -181,7 +181,7 @@ def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, pe return error #================================================================================== -def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter): +def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=1): #------------------------------------------------------------------------------ # Symbolic objects: SymPDE @@ -208,7 +208,7 @@ def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, #------------------------------------------------------------------------------ domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=degree) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity = [mult,mult,mult]) a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) a2_h = discretize(a2, domain_h, (derham_h.V2, derham_h.V2), backend=PSYDAC_BACKEND_GPYCCEL) @@ -356,6 +356,46 @@ class CollelaMapping3D(Mapping): error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter) assert abs(error - 0.24586986658559362) < 1e-9 + +#------------------------------------------------------------------------------ +def test_maxwell_3d_2_mult(): + class CollelaMapping3D(Mapping): + + _expressions = {'x': 'k1*(x1 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'y': 'k2*(x2 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'z': 'k3*x3'} + + _ldim = 3 + _pdim = 3 + + M = CollelaMapping3D('M', k1=1, k2=1, k3=1, eps=0.1) + logical_domain = Cube('C', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)) + + # exact solution + e_ex_0 = lambda t, x, y, z: 0 + e_ex_1 = lambda t, x, y, z: -np.cos(2*np.pi*t-2*np.pi*z) + e_ex_2 = lambda t, x, y, z: 0 + + e_ex = (e_ex_0, e_ex_1, e_ex_2) + + b_ex_0 = lambda t, x, y, z : np.cos(2*np.pi*t-2*np.pi*z) + b_ex_1 = lambda t, x, y, z : 0 + b_ex_2 = lambda t, x, y, z : 0 + + b_ex = (b_ex_0, b_ex_1, b_ex_2) + + #space parameters + ncells = [7, 7, 7] + degree = [2, 2, 2] + periodic = [True, True, True] + + #time parameters + dt = 0.5*1/max(ncells) + niter = 2 + 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.24749763720543216) < 1e-9 #============================================================================== # CLEAN UP SYMPY NAMESPACE @@ -368,4 +408,8 @@ def teardown_module(): def teardown_function(): from sympy.core import cache cache.clear_cache() + +if __name__ == '__main__' : + test_maxwell_3d_2_mult() + diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 4bbb9a217..3567b1f1a 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -301,7 +301,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): +def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, multiplicity = 1): """Computes the collocation matrix If called with normalization='M', this uses M-splines instead of B-splines. @@ -326,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 ------- @@ -342,19 +346,20 @@ 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: assert out.shape == ((xgrid.shape[0], nb)) and out.dtype == np.dtype('float') bool_normalization = normalization == "M" + multiplicity = int(multiplicity) - collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, out) + collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, out, multiplicity=multiplicity) return out #============================================================================== -def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_boundary=True, out=None): +def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multiplicity=1, check_boundary=True, out=None): """Computes the histopolation matrix. If called with normalization='M', this uses M-splines instead of B-splines. @@ -375,6 +380,10 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_bo 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``. @@ -418,23 +427,23 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_bo knots = np.ascontiguousarray(knots, dtype=float) xgrid = np.ascontiguousarray(xgrid, dtype=float) - elevated_knots = elevate_knots(knots, degree, periodic) + elevated_knots = elevate_knots(knots, degree, periodic, multiplicity=multiplicity) normalization = normalization == "M" 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') - - histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, check_boundary, elevated_knots, out) + multiplicity = int(multiplicity) + histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, check_boundary, elevated_knots, out, multiplicity = multiplicity) return out #============================================================================== @@ -472,7 +481,7 @@ def breakpoints(knots, degree, tol=1e-15, out=None): return out[:i_final] #============================================================================== -def greville(knots, degree, periodic, out=None): +def greville(knots, degree, periodic, out=None, multiplicity=1): """ Compute coordinates of all Greville points. @@ -490,6 +499,10 @@ def greville(knots, degree, periodic, 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 ------- @@ -499,9 +512,10 @@ def greville(knots, degree, periodic, out=None): """ 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) - greville_p(knots, degree, periodic, out) + multiplicity = int(multiplicity) + greville_p(knots, degree, periodic, out, multiplicity) return out #=============================================================================== @@ -560,9 +574,12 @@ def elements_spans(knots, degree, out=None): def make_knots(breaks, degree, periodic, multiplicity=1, out=None): """ Create spline knots from breakpoints, with appropriate boundary conditions. - Let p be spline degree. If domain is periodic, knot sequence is extended - by periodicity so that first p basis functions are identical to last p. - Otherwise, knot sequence is clamped (i.e. endpoints are repeated p times). + + If domain is periodic, knot sequence is extended by periodicity to have a + total of (n_cells-1)*mult+2p+2 knots (all break points are repeated mult + time and we add p+1-mult knots by periodicity at each side). + + Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). Parameters ---------- diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index dce0d059b..a2d9baf51 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -34,8 +34,9 @@ def find_span_p(knots: 'float[:]', degree: int, x: float): .. [1] L. Piegl and W. Tiller. The NURBS Book, 2nd ed., Springer-Verlag Berlin Heidelberg GmbH, 1997. """ - # Knot index at left/right boundary + # last knot on the left boundary low = degree + # first knot on the right boundary high = len(knots)-1-degree # Check if point is exactly on left/right boundary, or outside domain @@ -383,7 +384,7 @@ 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[:,:]'): + 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`. @@ -410,11 +411,16 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali 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. + """ # 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) @@ -460,7 +466,7 @@ 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[:,:]'): + 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. @@ -488,6 +494,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 ----- @@ -497,7 +507,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) @@ -545,7 +555,8 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma False, False, xgrid_new[:actual_len], - colloc) + colloc, + multiplicity = multiplicity) m = colloc.shape[0] - 1 n = colloc.shape[1] - 1 @@ -687,7 +698,7 @@ def breakpoints_p(knots: 'float[:]', degree: int, out: 'float[:]', tol: float = # ============================================================================= -def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]'): +def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]', multiplicity: int=1): """ Compute coordinates of all Greville points. @@ -705,15 +716,19 @@ def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]'): 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 == 0: + if p == multiplicity-1: for i in range(n): - out[i] = sum(T[i:i + 2]) / 2 + out[i] = sum(T[i:i + p + 2]) / (p + 2) else: for i in range(1, 1+n): out[i - 1] = sum(T[i:i + p]) / p @@ -786,9 +801,12 @@ def elements_spans_p(knots: 'float[:]', degree: int, out: 'int[:]'): def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:]', multiplicity: int = 1): """ Create spline knots from breakpoints, with appropriate boundary conditions. - Let p be spline degree. If domain is periodic, knot sequence is extended - by periodicity so that first p basis functions are identical to last p. - Otherwise, knot sequence is clamped (i.e. endpoints are repeated p times). + + If domain is periodic, knot sequence is extended by periodicity to have a + total of (n_cells-1)*mult+2p+2 knots (all break points are repeated mult + time and we add p+1-mult knots by periodicity at each side). + + Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). Parameters ---------- @@ -811,20 +829,23 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] It should be of the appropriate shape and dtype. """ ncells = len(breaks) - 1 - for i in range(1, ncells): - out[degree + 1 + (i - 1) * multiplicity:degree + 1 + i * multiplicity] = breaks[i] - - out[degree] = breaks[0] - out[len(out) - degree - 1] = breaks[-1] + + for i in range(0, ncells+1): + out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] + + len_out = len(out) if periodic: period = breaks[-1]-breaks[0] - out[:degree] = breaks[ncells - degree:ncells] - period - out[len(out) - degree:] = breaks[1:degree + 1] + period + out[: degree + 1 - multiplicity] = out[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period + out[len_out - (degree + 1 - multiplicity) :] = out[degree + 1:2*(degree + 1)- multiplicity] + period + + + # else: - out[0:degree + 1] = breaks[0] - out[len(out) - degree - 1:] = breaks[-1] + out[0:degree + 1 - multiplicity] = breaks[0] + out[len_out - degree - 1 + multiplicity:] = breaks[-1] # ============================================================================= @@ -866,8 +887,8 @@ def elevate_knots_p(knots: 'float[:]', degree: int, periodic: bool, out: 'float[ if periodic: T, p = knots, degree period = T[len(knots) -1 - p] - T[p] - left = T[len(knots) -2 - 2 * p] - period - right = T[2 * p + 1] + period + left = T[len(knots) -2 - 2 * p + multiplicity-1] - period + right = T[2 * p + 2 - multiplicity] + period out[0] = left out[-1] = right diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index f1bc9129b..5d9cdaad2 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -375,21 +375,24 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): if periodic: assert len(breaks) > degree - - p = degree - T = np.zeros( multiplicity*len(breaks[1:-1])+2+2*p ) - - T[p+1:-p-1] = np.repeat(breaks[1:-1], multiplicity) - T[p] = breaks[ 0] - T[-p-1] = breaks[-1] - + + 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[0:p] = [xi-period for xi in breaks[-p-1:-1 ]] - T[-p:] = [xi+period for xi in breaks[ 1: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[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 @@ -398,10 +401,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:]] @@ -491,6 +494,7 @@ 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())) + def test_find_span(knots, degree, x): expected = find_span_true(knots, degree, x) out = find_span(knots, degree, x) diff --git a/psydac/feec/derivatives.py b/psydac/feec/derivatives.py index 28fec0c11..94dd2643d 100644 --- a/psydac/feec/derivatives.py +++ b/psydac/feec/derivatives.py @@ -101,8 +101,8 @@ def __init__(self, V, W, diffdir, *, negative=False, transposed=False): self._codomain = W # the local area in the codomain without padding - self._idslice = tuple([slice(pad, e-s+1+pad) for pad, s, e - in zip(self._codomain.pads, self._codomain.starts, self._codomain.ends)]) + self._idslice = tuple([slice(pad*m, e-s+1+pad*m) for pad, s, e, m + in zip(self._codomain.pads, self._codomain.starts, self._codomain.ends, self._codomain.shifts)]) # prepare the slices (they are of the right size then, we checked this already) # identity slice @@ -112,12 +112,13 @@ def __init__(self, V, W, diffdir, *, negative=False, transposed=False): diff_pad = self._codomain.pads[self._diffdir] diff_s = self._codomain.starts[self._diffdir] diff_e = self._codomain.ends[self._diffdir] + diff_m = self._codomain.shifts[self._diffdir] # the diffslice depends on the transposition if self._transposed: - diff_partslice = slice(diff_pad-1, diff_e-diff_s+1+diff_pad-1) + diff_partslice = slice(diff_m*diff_pad-1, diff_e-diff_s+1+diff_pad*diff_m-1) else: - diff_partslice = slice(diff_pad+1, diff_e-diff_s+1+diff_pad+1) + diff_partslice = slice(diff_m*diff_pad+1, diff_e-diff_s+1+diff_m*diff_pad+1) diffslice = tuple([diff_partslice if i==self._diffdir else idslice[i] for i in range(self._domain.ndim)]) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 9106a9cbb..2584c391b 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -170,7 +170,7 @@ def __init__(self, space, nquads = None): solverblocks += [KroneckerLinearSolver(tensorspaces[i].vector_space, tensorspaces[i].vector_space, solvercells)] - dataslice = tuple(slice(p, -p) for p in tensorspaces[i].vector_space.pads) + dataslice = tuple(slice(p*m, -p*m) for p, m in zip(tensorspaces[i].vector_space.pads,tensorspaces[i].vector_space.shifts)) dofs[i] = rhsblocks[i]._data[dataslice] # finish arguments and create a lambda diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 44f7b91cc..10fdf699a 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -12,6 +12,8 @@ from psydac.feec.derivatives import ScalarCurl_2D, VectorCurl_2D, Curl_3D from psydac.feec.derivatives import Divergence_2D, Divergence_3D +from psydac.feec.global_projectors import Projector_H1 + from psydac.ddm.cart import DomainDecomposition from mpi4py import MPI @@ -348,13 +350,14 @@ def test_directional_derivative_operator_3d_par(domain, ncells, degree, periodic @pytest.mark.parametrize('degree', [2, 3, 4, 5]) @pytest.mark.parametrize('periodic', [True, False]) @pytest.mark.parametrize('seed', [1,3]) +@pytest.mark.parametrize('multiplicity', [1,2]) -def test_Derivative_1D(domain, ncells, degree, periodic, seed): +def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests np.random.seed(seed) breaks = np.linspace(*domain, num=ncells+1) - knots = make_knots(breaks, degree, periodic) + knots = make_knots(breaks, degree, periodic, multiplicity=multiplicity) # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') @@ -365,12 +368,13 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed): # L2 space (1-forms) V1 = V0.reduce_degree(axes=[0], basis='M') + # Create random field in V0 + u0 = FemField(V0) + # Linear operator: 1D derivative grad = Derivative_1D(V0, V1) # Create random field in V0 - u0 = FemField(V0) - s, = V0.vector_space.starts e, = V0.vector_space.ends @@ -395,8 +399,9 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) +@pytest.mark.parametrize('multiplicity', [(1, 1), (1, 2), (2, 2)]) -def test_Gradient_2D(domain, ncells, degree, periodic, seed): +def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests np.random.seed(seed) @@ -404,8 +409,8 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed): breaks = [np.linspace(*lims, num=n+1) for lims, n in zip(domain, ncells)] # H1 space (0-forms) - Nx, Ny = [SplineSpace(degree=d, grid=g, periodic=p, basis='B') \ - for d, g, p in zip(degree, breaks, periodic)] + Nx, Ny = [SplineSpace(degree=d, grid=g, periodic=p, basis='B', multiplicity=m) \ + for d, g, p, m in zip(degree, breaks, periodic, multiplicity)] domain_decomposition = DomainDecomposition(ncells, periodic) V0 = TensorFemSpace(domain_decomposition, Nx, Ny) @@ -456,8 +461,9 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed): (False, True, False), (False, False, True)]) @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): +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 @@ -467,9 +473,11 @@ def test_Gradient_3D(domain, ncells, degree, periodic, seed): # Compute breakpoints along each direction breaks = [np.linspace(*lims, num=n+1) for lims, n in zip(domain, ncells)] + #change multiplicity if higher than degree to avoid problems (case p1e-15)[0] + indices = np.where(np.diff(knots[degree:len(knots)-degree])>1e-15)[0] if len(indices)>0: multiplicity = np.diff(indices).max(initial=1) @@ -107,7 +107,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 @@ -131,10 +131,9 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._nbasis = nbasis self._breaks = grid self._ncells = len(grid) - 1 - self._greville = greville(knots, degree, periodic) - self._ext_greville = greville(elevate_knots(knots, degree, periodic), degree+1, periodic) + 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) @@ -178,7 +177,8 @@ def init_interpolation( self, dtype=float ): degree = self.degree, periodic = self.periodic, normalization = self.basis, - xgrid = self.greville + xgrid = self.greville, + multiplicity = self.multiplicity ) if self.periodic: @@ -212,8 +212,10 @@ def init_histopolation( self, dtype=float): degree = self.degree, periodic = self.periodic, normalization = self.basis, - xgrid = self.ext_greville + xgrid = self.ext_greville, + multiplicity = self._multiplicity ) + self.hmat= imat if self.periodic: # Convert to CSC format and compute sparse LU decomposition