make multipatch examples run

psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py

@@ -25,6 +25,7 @@
 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.utilities                   import time_count
+from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
 
 from psydac.linalg.utilities import array_to_psydac
 from psydac.fem.basic        import FemField
@@ -92,7 +93,7 @@ def solve_h1_source_pbm(
 
     print('building the symbolic and discrete deRham sequences...')
     derham  = Derham(domain, ["H1", "Hcurl", "L2"])
-    derham_h = discretize(derham, domain_h, degree=degree, backend=PSYDAC_BACKENDS[backend_language])
+    derham_h = discretize(derham, domain_h, degree=degree)
 
     # multi-patch (broken) spaces
     V0h = derham_h.V0
@@ -128,10 +129,8 @@ def solve_h1_source_pbm(
 
     print('conforming projection operators...')
     # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence)
-    cP0 = derham_h.conforming_projection(space='V0', hom_bc=True, backend_language=backend_language)
-    cP0_m = cP0.to_sparse_matrix()
-    # cP1 = derham_h.conforming_projection(space='V1', hom_bc=True, backend_language=backend_language)
-    # cP1_m = cP1.to_sparse_matrix()
+    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)
@@ -141,7 +140,7 @@ def lift_u_bc(u_bc):
             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) for m in mappings_list]
+            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()
@@ -176,7 +175,7 @@ def lift_u_bc(u_bc):
     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) for m in mappings_list]
+        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)
@@ -186,7 +185,7 @@ def lift_u_bc(u_bc):
         v  = element_of(V0h.symbolic_space, name='v')
         expr = f_scal * v
         l = LinearForm(v, integral(domain, expr))
-        lh = discretize(l, domain_h, V0h, backend=PSYDAC_BACKENDS[backend_language])
+        lh = discretize(l, domain_h, V0h)
         b  = lh.assemble()
         b_c = b.toarray()
         if plot_source:

psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py

@@ -18,6 +18,7 @@
 from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain
 from psydac.feec.multipatch.plotting_utilities          import plot_field
 from psydac.feec.multipatch.utilities                   import time_count
+from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
 
 def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language='python', mu=1, nu=1, gamma_h=10,
                           sigma=None, nb_eigs=4, nb_eigs_plot=4,
@@ -71,7 +72,7 @@ def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language
 
     print('building symbolic and discrete derham sequences...')
     derham  = Derham(domain, ["H1", "Hcurl", "L2"])
-    derham_h = discretize(derham, domain_h, degree=degree, backend=PSYDAC_BACKENDS[backend_language])
+    derham_h = discretize(derham, domain_h, degree=degree)
 
     V0h = derham_h.V0
     V1h = derham_h.V1
@@ -102,10 +103,8 @@ def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language
 
     print('conforming projection operators...')
     # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence)
-    cP0 = derham_h.conforming_projection(space='V0', hom_bc=True, backend_language=backend_language, load_dir=m_load_dir)
-    cP1 = derham_h.conforming_projection(space='V1', hom_bc=True, backend_language=backend_language, load_dir=m_load_dir)
-    cP0_m = cP0.to_sparse_matrix()
-    cP1_m = cP1.to_sparse_matrix()
+    cP0_m = construct_V0_conforming_projection(V0h, domain_h, True)
+    cP1_m = construct_V1_conforming_projection(V1h, domain_h, True)
 
     print('broken differential operators...')
     bD0, bD1 = derham_h.broken_derivatives_as_operators
@@ -214,10 +213,10 @@ def get_eigenvalues(nb_eigs, sigma, A_m, M_m):
         nc = 8
         deg = 4
 
-    domain_name = 'pretzel_f'
-    # domain_name = 'curved_L_shape'
+    #domain_name = 'pretzel_f'
+    domain_name = 'curved_L_shape'
     nc = 10
-    deg = 2
+    deg = 3
 
     m_load_dir = 'matrices_{}_nc={}_deg={}/'.format(domain_name, nc, deg)
     run_dir = 'eigenpbm_{}_nc={}_deg={}/'.format(domain_name, nc, deg)

psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py

@@ -29,6 +29,8 @@
 from psydac.linalg.utilities                            import array_to_psydac
 from psydac.fem.basic                                   import FemField
 
+from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
+
 def solve_hcurl_source_pbm(
         nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='P_geom', source_type='manu_J',
         eta=-10., mu=1., nu=1., gamma_h=10.,
@@ -107,7 +109,7 @@ def solve_hcurl_source_pbm(
 
     t_stamp = time_count(t_stamp)
     print('building discrete derham sequence...')
-    derham_h = discretize(derham, domain_h, degree=degree, backend=PSYDAC_BACKENDS[backend_language])
+    derham_h = discretize(derham, domain_h, degree=degree)
 
     t_stamp = time_count(t_stamp)
     print('building commuting projection operators...')
@@ -162,10 +164,8 @@ def solve_hcurl_source_pbm(
     t_stamp = time_count(t_stamp)
     print('building the conforming Projection operators and matrices...')
     # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence)
-    cP0 = derham_h.conforming_projection(space='V0', hom_bc=True, backend_language=backend_language, load_dir=m_load_dir)
-    cP1 = derham_h.conforming_projection(space='V1', hom_bc=True, backend_language=backend_language, load_dir=m_load_dir)
-    cP0_m = cP0.to_sparse_matrix()
-    cP1_m = cP1.to_sparse_matrix()
+    cP0_m = construct_V0_conforming_projection(V0h, domain_h, hom_bc=True)
+    cP1_m = construct_V1_conforming_projection(V1h, domain_h, hom_bc=True)
 
     t_stamp = time_count(t_stamp)
     print('building the broken differential operators and matrices...')
@@ -183,7 +183,7 @@ def lift_u_bc(u_bc):
             # note: for simplicity we apply the full P1 on u_bc, but we only need to set the boundary dofs
             u_bc_x = lambdify(domain.coordinates, u_bc[0])
             u_bc_y = lambdify(domain.coordinates, u_bc[1])
-            u_bc_log = [pull_2d_hcurl([u_bc_x, u_bc_y], m) for m in mappings_list]
+            u_bc_log = [pull_2d_hcurl([u_bc_x, u_bc_y], 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 = P1(u_bc_log)
             ubc_c = uh_bc.coeffs.toarray()
@@ -240,7 +240,7 @@ def lift_u_bc(u_bc):
         print('projecting the source with commuting projection...')
         f_x = lambdify(domain.coordinates, f_vect[0])
         f_y = lambdify(domain.coordinates, f_vect[1])
-        f_log = [pull_2d_hcurl([f_x, f_y], m) for m in mappings_list]
+        f_log = [pull_2d_hcurl([f_x, f_y], m.get_callable_mapping()) for m in mappings_list]
         f_h = P1(f_log)
         f_c = f_h.coeffs.toarray()
         b_c = dH1_m.dot(f_c)
@@ -251,7 +251,7 @@ def lift_u_bc(u_bc):
         v  = element_of(V1h.symbolic_space, name='v')
         expr = dot(f_vect,v)
         l = LinearForm(v, integral(domain, expr))
-        lh = discretize(l, domain_h, V1h, backend=PSYDAC_BACKENDS[backend_language])
+        lh = discretize(l, domain_h, V1h)
         b  = lh.assemble()
         b_c = b.toarray()
         if plot_source:

psydac/feec/multipatch/examples/mixed_source_pbms_conga_2d.py

@@ -28,6 +28,8 @@
 from psydac.feec.multipatch.examples.hcurl_eigen_pbms_conga_2d import get_eigenvalues
 from psydac.feec.multipatch.utilities                          import time_count
 
+from psydac.feec.multipatch.non_matching_operators      import construct_V0_conforming_projection, construct_V1_conforming_projection
+
 def solve_magnetostatic_pbm(
         nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='P_L2_wcurl_J',
         source_type='dipole_J', bc_type='metallic',
@@ -120,7 +122,7 @@ def solve_magnetostatic_pbm(
 
     print('building symbolic and discrete derham sequences...')
     derham  = Derham(domain, ["H1", "Hcurl", "L2"])
-    derham_h = discretize(derham, domain_h, degree=degree, backend=PSYDAC_BACKENDS[backend_language])
+    derham_h = discretize(derham, domain_h, degree=degree)
 
     V0h = derham_h.V0
     V1h = derham_h.V1
@@ -137,18 +139,18 @@ def solve_magnetostatic_pbm(
     # these physical projection operators should probably be in the interface...
     def P0_phys(f_phys):
         f = lambdify(domain.coordinates, f_phys)
-        f_log = [pull_2d_h1(f, m) for m in mappings_list]
+        f_log = [pull_2d_h1(f, m.get_callable_mapping()) for m in mappings_list]
         return P0(f_log)
 
     def P1_phys(f_phys):
         f_x = lambdify(domain.coordinates, f_phys[0])
         f_y = lambdify(domain.coordinates, f_phys[1])
-        f_log = [pull_2d_hcurl([f_x, f_y], m) for m in mappings_list]
+        f_log = [pull_2d_hcurl([f_x, f_y], m.get_callable_mapping()) for m in mappings_list]
         return P1(f_log)
 
     def P2_phys(f_phys):
         f = lambdify(domain.coordinates, f_phys)
-        f_log = [pull_2d_l2(f, m) for m in mappings_list]
+        f_log = [pull_2d_l2(f, m.get_callable_mapping()) for m in mappings_list]
         return P2(f_log)
 
     I0_m = IdLinearOperator(V0h).to_sparse_matrix()
@@ -175,10 +177,9 @@ def P2_phys(f_phys):
 
     print('conforming projection operators...')
     # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence)
-    cP0 = derham_h.conforming_projection(space='V0', hom_bc=hom_bc, backend_language=backend_language, load_dir=m_load_dir)
-    cP1 = derham_h.conforming_projection(space='V1', hom_bc=hom_bc, backend_language=backend_language, load_dir=m_load_dir)
-    cP0_m = cP0.to_sparse_matrix()
-    cP1_m = cP1.to_sparse_matrix()
+    cP0_m = construct_V0_conforming_projection(V0h, domain_h, hom_bc=True)
+    cP1_m = construct_V1_conforming_projection(V1h, domain_h, hom_bc=True)
+
 
     print('broken differential operators...')
     bD0, bD1 = derham_h.broken_derivatives_as_operators

psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py

@@ -16,7 +16,7 @@ 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
 
 #==============================================================================
@@ -30,3 +30,6 @@ def teardown_module():
 def teardown_function():
     from sympy.core import cache
     cache.clear_cache()
+
+if __name__ == '__main__':
+    test_poisson_pretzel_f()