add new tests and last example

psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py

@@ -23,7 +23,7 @@
 from psydac.feec.multipatch.operators                   import HodgeOperator
 from psydac.feec.multipatch.plotting_utilities          import plot_field
 from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain
-from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution
+from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution_OBSOLETE
 from psydac.feec.multipatch.utilities                   import time_count
 from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
 
@@ -165,7 +165,7 @@ def lift_u_bc(u_bc):
     # (not all the returned functions are useful here)
     N_diag = 200
     method = 'conga'
-    f_scal, f_vect, u_bc, ph_ref, uh_ref, p_ex, u_ex, phi, grad_phi = get_source_and_solution(
+    f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE(
         source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,
         refsol_params=[N_diag, method, source_proj],
     )

psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py

@@ -24,7 +24,7 @@
 from psydac.feec.multipatch.operators                   import HodgeOperator
 from psydac.feec.multipatch.plotting_utilities          import plot_field
 from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain
-from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution
+from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution_OBSOLETE
 from psydac.feec.multipatch.utilities                   import time_count
 from psydac.linalg.utilities                            import array_to_psydac
 from psydac.fem.basic                                   import FemField
@@ -227,9 +227,8 @@ def lift_u_bc(u_bc):
     print('getting the source and ref solution...')
     N_diag = 200
     method = 'conga'
-    f_scal, f_vect, u_bc, ph_ref, uh_ref, p_ex, u_ex, phi, grad_phi = get_source_and_solution(
+    f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE(
         source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,
-        refsol_params=[N_diag, method, source_proj],
     )
 
     # compute approximate source f_h

psydac/feec/multipatch/examples_nc/h1_source_pbms_nc.py

@@ -0,0 +1,278 @@
+# coding: utf-8
+
+from mpi4py import MPI
+
+import os
+import numpy as np
+from collections import OrderedDict
+
+from sympy import lambdify
+from scipy.sparse.linalg import spsolve
+
+from sympde.expr.expr import LinearForm
+from sympde.expr.expr import integral, Norm
+from sympde.topology  import Derham
+from sympde.topology import element_of
+
+
+from psydac.api.settings        import PSYDAC_BACKENDS
+from psydac.feec.multipatch.api import discretize
+from psydac.feec.pull_push      import pull_2d_h1
+
+from psydac.feec.multipatch.fem_linear_operators        import IdLinearOperator
+from psydac.feec.multipatch.operators                   import HodgeOperator
+from psydac.feec.multipatch.plotting_utilities          import plot_field
+from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain
+from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution_OBSOLETE
+from psydac.feec.multipatch.utilities                   import time_count
+from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
+from psydac.api.postprocessing import OutputManager, PostProcessManager
+
+from psydac.linalg.utilities import array_to_psydac
+from psydac.fem.basic        import FemField
+
+def solve_h1_source_pbm_nc(
+        nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='P_L2', source_type='manu_poisson',
+        eta=-10., mu=1., gamma_h=10.,
+        plot_source=False, plot_dir=None, hide_plots=True
+):
+    """
+    solver for the problem: find u in H^1, such that
+
+      A u = f             on \Omega
+        u = u_bc          on \partial \Omega
+
+    where the operator
+
+      A u := eta * u  -  mu * div grad u
+
+    is discretized as  Ah: V0h -> V0h  in a broken-FEEC approach involving a discrete sequence on a 2D multipatch domain \Omega,
+
+      V0h  --grad->  V1h  -—curl-> V2h
+
+    Examples:
+
+      - Helmholtz equation with
+          eta = -omega**2
+          mu  = 1
+
+      - Poisson equation with Laplace operator L = A,
+          eta = 0
+          mu  = 1
+
+    :param nc: nb of cells per dimension, in each patch
+    :param deg: coordinate degree in each patch
+    :param gamma_h: jump penalization parameter
+    :param source_proj: approximation operator for the source, possible values are 'P_geom' or 'P_L2'
+    :param source_type: must be implemented in get_source_and_solution()
+    """
+
+    ncells = nc
+    degree = [deg,deg]
+
+    # if backend_language is None:
+    #     if domain_name in ['pretzel', 'pretzel_f'] and nc > 8:
+    #         backend_language='numba'
+    #     else:
+    #         backend_language='python'
+    # print('[note: using '+backend_language+ ' backends in discretize functions]')
+
+    print('---------------------------------------------------------------------------------------------------------')
+    print('Starting solve_h1_source_pbm function with: ')
+    print(' ncells = {}'.format(ncells))
+    print(' degree = {}'.format(degree))
+    print(' domain_name = {}'.format(domain_name))
+    print(' source_proj = {}'.format(source_proj))
+    print(' backend_language = {}'.format(backend_language))
+    print('---------------------------------------------------------------------------------------------------------')
+
+    print('building the multipatch domain...')
+    domain = build_multipatch_domain(domain_name=domain_name)
+    mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior])
+    mappings_list = list(mappings.values())
+    ncells_h = {patch.name: [ncells[i], ncells[i]] for (i,patch) in enumerate(domain.interior)}
+    domain_h = discretize(domain, ncells=ncells_h)
+
+    print('building the symbolic and discrete deRham sequences...')
+    derham  = Derham(domain, ["H1", "Hcurl", "L2"])
+    derham_h = discretize(derham, domain_h, degree=degree)
+
+    # multi-patch (broken) spaces
+    V0h = derham_h.V0
+    V1h = derham_h.V1
+    V2h = derham_h.V2
+    print('dim(V0h) = {}'.format(V0h.nbasis))
+    print('dim(V1h) = {}'.format(V1h.nbasis))
+    print('dim(V2h) = {}'.format(V2h.nbasis))
+
+    print('broken differential operators...')
+    # broken (patch-wise) differential operators
+    bD0, bD1 = derham_h.broken_derivatives_as_operators
+    bD0_m = bD0.to_sparse_matrix()
+    # bD1_m = bD1.to_sparse_matrix()
+
+    print('building the discrete operators:')
+    print('commuting projection operators...')
+    nquads = [4*(d + 1) for d in degree]
+    P0, P1, P2 = derham_h.projectors(nquads=nquads)
+
+    I0 = IdLinearOperator(V0h)
+    I0_m = I0.to_sparse_matrix()
+
+    print('Hodge operators...')
+    # multi-patch (broken) linear operators / matrices
+    H0 = HodgeOperator(V0h, domain_h, backend_language=backend_language)
+    H1 = HodgeOperator(V1h, domain_h, backend_language=backend_language)
+
+    dH0_m = H0.get_dual_Hodge_sparse_matrix()  # = mass matrix of V0
+    H0_m  = H0.to_sparse_matrix()              # = inverse mass matrix of V0
+    dH1_m = H1.get_dual_Hodge_sparse_matrix()  # = mass matrix of V1
+    # H1_m  = H1.to_sparse_matrix()              # = inverse mass matrix of V1
+
+    print('conforming projection operators...')
+    # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence)
+    cP0_m = construct_V0_conforming_projection(V0h, domain_h, hom_bc=True)
+    # cP1_m = construct_V1_conforming_projection(V1h, domain_h, hom_bc=True)
+
+    if not os.path.exists(plot_dir):
+        os.makedirs(plot_dir)
+
+    def lift_u_bc(u_bc):
+        if u_bc is not None:
+            print('lifting the boundary condition in V0h...  [warning: Not Tested Yet!]')
+            # note: for simplicity we apply the full P1 on u_bc, but we only need to set the boundary dofs
+            u_bc = lambdify(domain.coordinates, u_bc)
+            u_bc_log = [pull_2d_h1(u_bc, m.get_callable_mapping()) for m in mappings_list]
+            # it's a bit weird to apply P1 on the list of (pulled back) logical fields -- why not just apply it on u_bc ?
+            uh_bc = P0(u_bc_log)
+            ubc_c = uh_bc.coeffs.toarray()
+            # removing internal dofs (otherwise ubc_c may already be a very good approximation of uh_c ...)
+            ubc_c = ubc_c - cP0_m.dot(ubc_c)
+        else:
+            ubc_c = None
+        return ubc_c
+
+    # Conga (projection-based) stiffness matrices:
+    # div grad:
+    pre_DG_m = - bD0_m.transpose() @ dH1_m @ bD0_m
+
+    # jump penalization:
+    jump_penal_m = I0_m - cP0_m
+    JP0_m = jump_penal_m.transpose() * dH0_m * jump_penal_m
+
+    pre_A_m = cP0_m.transpose() @ ( eta * dH0_m - mu * pre_DG_m )  # useful for the boundary condition (if present)
+    A_m = pre_A_m @ cP0_m + gamma_h * JP0_m
+
+    print('getting the source and ref solution...')
+    # (not all the returned functions are useful here)
+    N_diag = 200
+    method = 'conga'
+    f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE(
+        source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,
+        refsol_params=[N_diag, method, source_proj],
+    )
+
+    # compute approximate source f_h
+    b_c = f_c = None
+    if source_proj == 'P_geom':
+        print('projecting the source with commuting projection P0...')
+        f = lambdify(domain.coordinates, f_scal)
+        f_log = [pull_2d_h1(f, m.get_callable_mapping()) for m in mappings_list]
+        f_h = P0(f_log)
+        f_c = f_h.coeffs.toarray()
+        b_c = dH0_m.dot(f_c)
+
+    elif source_proj == 'P_L2':
+        print('projecting the source with L2 projection...')
+        v  = element_of(V0h.symbolic_space, name='v')
+        expr = f_scal * v
+        l = LinearForm(v, integral(domain, expr))
+        lh = discretize(l, domain_h, V0h)
+        b  = lh.assemble()
+        b_c = b.toarray()
+        if plot_source:
+            f_c = H0_m.dot(b_c)
+    else:
+        raise ValueError(source_proj)
+
+    if plot_source:
+        plot_field(numpy_coeffs=f_c, Vh=V0h, space_kind='h1', domain=domain, title='f_h with P = '+source_proj, filename=plot_dir+'/fh_'+source_proj+'.png', hide_plot=hide_plots)
+
+    ubc_c = lift_u_bc(u_bc)
+
+    if ubc_c is not None:
+        # modified source for the homogeneous pbm
+        print('modifying the source with lifted bc solution...')
+        b_c = b_c - pre_A_m.dot(ubc_c)
+
+    # direct solve with scipy spsolve
+    print('solving source problem with scipy.spsolve...')
+    uh_c = spsolve(A_m, b_c)
+
+    # project the homogeneous solution on the conforming problem space
+    print('projecting the homogeneous solution on the conforming problem space...')
+    uh_c = cP0_m.dot(uh_c)
+
+    if ubc_c is not None:
+        # adding the lifted boundary condition
+        print('adding the lifted boundary condition...')
+        uh_c += ubc_c
+
+    print('getting and plotting the FEM solution from numpy coefs array...')
+    title = r'solution $\phi_h$ (amplitude)'
+    params_str = 'eta={}_mu={}_gamma_h={}'.format(eta, mu, gamma_h)
+    plot_field(numpy_coeffs=uh_c, Vh=V0h, space_kind='h1', domain=domain, title=title, filename=plot_dir+params_str+'_phi_h.png', hide_plot=hide_plots)
+
+
+    if u_ex:
+        u         = element_of(V0h.symbolic_space, name='u')
+        l2norm    = Norm(u - u_ex, domain, kind='l2')
+        l2norm_h  = discretize(l2norm, domain_h, V0h)
+        uh_c      = array_to_psydac(uh_c, V0h.vector_space)
+        l2_error  = l2norm_h.assemble(u=FemField(V0h, coeffs=uh_c))
+        return l2_error
+
+if __name__ == '__main__':
+
+    t_stamp_full = time_count()
+
+    quick_run = True
+    # quick_run = False
+
+    omega = np.sqrt(170) # source
+    roundoff = 1e4
+    eta = int(-omega**2 * roundoff)/roundoff
+    # print(eta)
+    # source_type = 'elliptic_J'
+    source_type = 'manu_poisson'
+
+    # if quick_run:
+    #     domain_name = 'curved_L_shape'
+    #     nc = 4
+    #     deg = 2
+    # else:
+    #     nc = 8
+    #     deg = 4
+
+    domain_name = 'pretzel_f'
+    # domain_name = 'curved_L_shape'
+    nc = np.array([8, 8, 16, 16, 8, 4, 4, 4, 4, 4, 2, 2, 4, 16, 16, 8, 2, 2, 2])
+    deg = 2
+
+    # nc = 2
+    # deg = 2
+
+    run_dir = '{}_{}_nc={}_deg={}/'.format(domain_name, source_type, nc, deg)
+    solve_h1_source_pbm_nc(
+        nc=nc, deg=deg,
+        eta=eta,
+        mu=1, #1,
+        domain_name=domain_name,
+        source_type=source_type,
+        backend_language='pyccel-gcc',
+        plot_source=True,
+        plot_dir='./plots/h1_tests_source_february/'+run_dir,
+        hide_plots=True,
+    )
+
+    time_count(t_stamp_full, msg='full program')

psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py

@@ -30,7 +30,7 @@
 from psydac.feec.multipatch.non_matching_multipatch_domain_utilities import create_square_domain
 from psydac.api.postprocessing import OutputManager, PostProcessManager
 
-def hcurl_solve_eigen_pbm_multipatch_dg(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0,
+def hcurl_solve_eigen_pbm_dg(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0,
                           generalized_pbm=False, sigma=None, ref_sigmas=[], nb_eigs_solve=8, nb_eigs_plot=5, skip_eigs_threshold=1e-7,
                           plot_dir=None, hide_plots=True, m_load_dir="",):
 

psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py

@@ -33,9 +33,9 @@
 from psydac.api.postprocessing import OutputManager, PostProcessManager
 
 
-def hcurl_solve_eigen_pbm_multipatch_nc(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0,
+def hcurl_solve_eigen_pbm_nc(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0,
                           generalized_pbm=False, sigma=None, ref_sigmas=[], nb_eigs_solve=8, nb_eigs_plot=5, skip_eigs_threshold=1e-7,
-                          plot_dir=None, hide_plots=True, m_load_dir="",):
+                          plot_dir=None, hide_plots=True, m_load_dir=None,):
 
     diags = {}
 

psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py

@@ -2,8 +2,8 @@
 import numpy as np
 
 
-from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_multipatch_nc
-from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_multipatch_dg
+from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_nc
+from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_dg
 
 from psydac.feec.multipatch.utilities                   import time_count, get_run_dir, get_plot_dir, get_mat_dir, get_sol_dir, diag_fn
 from psydac.feec.multipatch.utils_conga_2d              import write_diags_to_file
@@ -235,7 +235,7 @@
 #   - we look for nb_eigs_solve eigenvalues close to sigma (skip zero eigenvalues if skip_zero_eigs==True)
 #   - we plot nb_eigs_plot eigenvectors
 if method == 'feec':
-    diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_nc(
+    diags, eigenvalues = hcurl_solve_eigen_pbm_nc(
         ncells=ncells, degree=degree,
         gamma_h=gamma_h,
         generalized_pbm=generalized_pbm,
@@ -253,7 +253,7 @@
         m_load_dir=m_load_dir,
     )
 elif method == 'dg': 
-    diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_dg(
+    diags, eigenvalues = hcurl_solve_eigen_pbm_dg(
     ncells=ncells, degree=degree,
     gamma_h=gamma_h,
     generalized_pbm=generalized_pbm,

psydac/feec/multipatch/examples_nc/hcurl_source_pbms_nc.py

@@ -29,6 +29,7 @@
 from psydac.feec.multipatch.utilities                   import time_count #, export_sol, import_sol
 from psydac.linalg.utilities                            import array_to_psydac
 from psydac.fem.basic                                   import FemField
+from psydac.feec.multipatch.examples.ppc_test_cases     import get_source_and_solution_OBSOLETE
 
 from psydac.feec.multipatch.non_matching_operators import construct_V0_conforming_projection, construct_V1_conforming_projection
 from psydac.api.postprocessing import OutputManager, PostProcessManager
@@ -39,8 +40,8 @@ def solve_hcurl_source_pbm_nc(
         project_sol=False,
         plot_source=False, plot_dir=None, hide_plots=True, skip_plot_titles=False,
         cb_min_sol=None, cb_max_sol=None,
-        m_load_dir="", sol_filename="", sol_ref_filename="",
-        ref_nc=None, ref_deg=None,
+        m_load_dir=None, sol_filename="", sol_ref_filename="",
+        ref_nc=None, ref_deg=None, test=False
 ):
     """
     solver for the problem: find u in H(curl), such that
@@ -246,9 +247,14 @@ def lift_u_bc(u_bc):
     t_stamp = time_count(t_stamp)
     print()
     print(' -- getting source --')
-    f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl(
+    if source_type == 'manu_maxwell':
+        f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE(
         source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,
-    )
+        )
+    else:
+        f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl(
+        source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,
+        )
     # compute approximate source f_h
     t_stamp = time_count(t_stamp)
     tilde_f_c = f_c = None
@@ -352,5 +358,14 @@ def lift_u_bc(u_bc):
 
     time_count(t_stamp)
 
+    if test:
+        u         = element_of(V1h.symbolic_space, name='u')
+        l2norm    = Norm(Matrix([u[0] - u_ex[0],u[1] - u_ex[1]]), domain, kind='l2')
+        l2norm_h  = discretize(l2norm, domain_h, V1h)
+        uh_c      = array_to_psydac(uh_c, V1h.vector_space)
+        l2_error  = l2norm_h.assemble(u=FemField(V1h, coeffs=uh_c))
+        print(l2_error)
+        return l2_error
+
 
     return diags

psydac/feec/multipatch/examples_nc/hcurl_source_testcase.py

@@ -16,7 +16,6 @@
 test_case = 'maxwell_hom_eta=170'   # used in paper
 # test_case = 'maxwell_inhom'   # used in paper
 
-compute_ref_sol = False  # (not needed for inhomogeneous test-case, as exact solution is known)
 
 #
 # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py

@@ -3,6 +3,11 @@
 import numpy as np
 
 from psydac.feec.multipatch.examples.hcurl_source_pbms_conga_2d import solve_hcurl_source_pbm
+from psydac.feec.multipatch.examples_nc.hcurl_source_pbms_nc import solve_hcurl_source_pbm_nc
+
+from psydac.feec.multipatch.examples.hcurl_eigen_pbms_conga_2d import hcurl_solve_eigen_pbm
+from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_nc
+from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_dg
 
 def test_time_harmonic_maxwell_pretzel_f():
     nc,deg      = 10,2
@@ -24,6 +29,161 @@ def test_time_harmonic_maxwell_pretzel_f():
 
     assert abs(l2_error - 0.06246693595198972)<1e-10
 
+def test_time_harmonic_maxwell_pretzel_f_nc():
+    deg      = 2
+    nc = np.array([20, 20, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 10, 10])
+
+    source_type = 'manu_maxwell'
+    domain_name = 'pretzel_f'
+    source_proj = 'tilde_Pi'
+
+    omega = np.sqrt(170) # source
+    roundoff = 1e4
+    eta = int(-omega**2 * roundoff)/roundoff
+
+    l2_error = solve_hcurl_source_pbm_nc(
+        nc=nc, deg=deg,
+        eta=eta,
+        nu=0,
+        mu=1, 
+        domain_name=domain_name,
+        source_type=source_type,
+        source_proj=source_proj,
+        plot_dir='./plots/th_maxell_nc',
+        backend_language='pyccel-gcc', 
+        test=True)
+
+    assert abs(l2_error - 0.04753982587323614)<1e-10
+
+def test_maxwell_eigen_curved_L_shape():
+    domain_name = 'curved_L_shape' 
+
+    nc = 10
+    deg = 2
+
+    ref_sigmas = [
+        0.181857115231E+01,
+        0.349057623279E+01,
+        0.100656015004E+02,
+        0.101118862307E+02,
+        0.124355372484E+02,
+        ]
+    sigma = 7
+    nb_eigs_solve = 7
+    nb_eigs_plot = 7
+    skip_eigs_threshold = 1e-7
+
+    eigenvalues, eigenvectors = hcurl_solve_eigen_pbm(
+        nc=nc, deg=deg,
+        gamma_h=0,
+        nu=0,
+        mu=1,
+        sigma=sigma,
+        skip_eigs_threshold=skip_eigs_threshold,
+        nb_eigs=nb_eigs_solve,
+        nb_eigs_plot=nb_eigs_plot,
+        domain_name=domain_name,
+        backend_language='pyccel-gcc',
+        plot_dir='./plots/eigen_maxell',
+    )
+
+    error = 0
+    n_errs = min(len(ref_sigmas), len(eigenvalues))
+    for k in range(n_errs):
+        error += (eigenvalues[k]-ref_sigmas[k])**2
+    error = np.sqrt(error)
+
+    assert abs(error - 0.023413963252245817)<1e-10
+
+def test_maxwell_eigen_curved_L_shape_nc():
+    domain_name = 'curved_L_shape' 
+    domain=[[1, 3],[0, np.pi/4]] 
+
+    ncells = np.array([[None, 10],
+                   [10,   20]])
+
+    degree = [2, 2]
+
+    ref_sigmas = [
+        0.181857115231E+01,
+        0.349057623279E+01,
+        0.100656015004E+02,
+        0.101118862307E+02,
+        0.124355372484E+02,
+        ]
+    sigma = 7
+    nb_eigs_solve = 7
+    nb_eigs_plot = 7
+    skip_eigs_threshold = 1e-7
+
+    diags, eigenvalues = hcurl_solve_eigen_pbm_nc(
+        ncells=ncells, degree=degree,
+        gamma_h=0,
+        generalized_pbm=True,
+        nu=0,
+        mu=1,
+        sigma=sigma,
+        ref_sigmas=ref_sigmas,
+        skip_eigs_threshold=skip_eigs_threshold,
+        nb_eigs_solve=nb_eigs_solve,
+        nb_eigs_plot=nb_eigs_plot,
+        domain_name=domain_name, domain=domain,
+        backend_language='pyccel-gcc',
+        plot_dir='./plots/eigen_maxell_nc',
+    )
+
+    error = 0
+    n_errs = min(len(ref_sigmas), len(eigenvalues))
+    for k in range(n_errs):
+        error += (eigenvalues[k]-ref_sigmas[k])**2
+    error = np.sqrt(error)
+
+    assert abs(error - 0.004289103786542442)<1e-10
+
+def test_maxwell_eigen_curved_L_shape_dg():
+    domain_name = 'curved_L_shape' 
+    domain=[[1, 3],[0, np.pi/4]] 
+
+    ncells = np.array([[None, 10],
+                   [10,   20]])
+
+    degree = [2, 2]
+
+    ref_sigmas = [
+        0.181857115231E+01,
+        0.349057623279E+01,
+        0.100656015004E+02,
+        0.101118862307E+02,
+        0.124355372484E+02,
+        ]
+    sigma = 7
+    nb_eigs_solve = 7
+    nb_eigs_plot = 7
+    skip_eigs_threshold = 1e-7
+
+    diags, eigenvalues = hcurl_solve_eigen_pbm_dg(
+        ncells=ncells, degree=degree,
+        gamma_h=0,
+        generalized_pbm=True,
+        nu=0,
+        mu=1,
+        sigma=sigma,
+        ref_sigmas=ref_sigmas,
+        skip_eigs_threshold=skip_eigs_threshold,
+        nb_eigs_solve=nb_eigs_solve,
+        nb_eigs_plot=nb_eigs_plot,
+        domain_name=domain_name, domain=domain,
+        backend_language='pyccel-gcc',
+        plot_dir='./plots/eigen_maxell_dg',
+    )
+
+    error = 0
+    n_errs = min(len(ref_sigmas), len(eigenvalues))
+    for k in range(n_errs):
+        error += (eigenvalues[k]-ref_sigmas[k])**2
+    error = np.sqrt(error)
+
+    assert abs(error - 0.004208158031148591)<1e-10
 
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
@@ -35,4 +195,4 @@ def teardown_module():
 
 def teardown_function():
     from sympy.core import cache
-    cache.clear_cache()
+    cache.clear_cache()

psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py

@@ -1,4 +1,7 @@
+import numpy as np
+
 from psydac.feec.multipatch.examples.h1_source_pbms_conga_2d import solve_h1_source_pbm
+from psydac.feec.multipatch.examples_nc.h1_source_pbms_nc import solve_h1_source_pbm_nc
 
 def test_poisson_pretzel_f():
 
@@ -16,9 +19,27 @@ def test_poisson_pretzel_f():
                 backend_language='pyccel-gcc',
                 plot_source=False,
                 plot_dir='./plots/h1_tests_source_february/'+run_dir)
-    print(l2_error)
-    assert abs(l2_error-8.054935880021907e-05)<1e-10
 
+    assert abs(l2_error-0.11860734907095004)<1e-10
+
+def test_poisson_pretzel_f_nc():
+
+    source_type = 'manu_poisson_2'
+    domain_name = 'pretzel_f'
+    nc  = np.array([20, 20, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 10, 10])
+    deg = 2
+    run_dir = '{}_{}_nc={}_deg={}/'.format(domain_name, source_type, nc, deg)
+    l2_error = solve_h1_source_pbm_nc(
+                nc=nc, deg=deg,
+                eta=0,
+                mu=1,
+                domain_name=domain_name,
+                source_type=source_type,
+                backend_language='pyccel-gcc',
+                plot_source=False,
+                plot_dir='./plots/h1_tests_source_february/'+run_dir)
+
+    assert abs(l2_error-0.04324704991715671)<1e-10
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
 #==============================================================================
@@ -29,7 +50,4 @@ def teardown_module():
 
 def teardown_function():
     from sympy.core import cache
-    cache.clear_cache()
-
-if __name__ == '__main__':
-    test_poisson_pretzel_f()
+    cache.clear_cache()