Non Matching Multipatch

psydac/api/fem.py

@@ -727,7 +727,9 @@ def allocate_matrices(self, backend=None):
                         elif use_prolongation:
                             Ps  = [knot_insertion_projection_operator(trs, trs.get_refined_space(ncells)) for trs in trial_fem_space.spaces]
                             P   = BlockLinearOperator(trial_fem_space.vector_space, trial_fem_space.get_refined_space(ncells).vector_space)
-                            for ni,Pi in enumerate(Ps):P[ni,ni] = Pi
+                            for ni,Pi in enumerate(Ps):
+                                P[ni,ni] = Pi
+
                             mat = ComposedLinearOperator(trial_space, test_space, mat, P)
 
                     self._matrix[i,j] = mat

psydac/api/postprocessing.py

@@ -953,12 +953,9 @@ def _reconstruct_spaces(self):
         for subdomain_names, space_dict in subdomains_to_spaces.items():
             if space_dict == {}:
                 continue
-            ncells_dict = {interior_name: interior_names_to_ncells[interior_name] for interior_name in subdomain_names}
-            # No need for a a dict until PR about non-conforming meshes is merged
-            # Check for conformity
-            ncells =  list(ncells_dict.values())[0]
-            assert all(ncells_patch == ncells for ncells_patch in ncells_dict.values())
 
+            ncells = {interior_name: interior_names_to_ncells[interior_name] for interior_name in subdomain_names}
+
             subdomain    = domain.get_subdomain(subdomain_names)
             space_name_0 = list(space_dict.keys())[0]
             periodic     = space_dict[space_name_0][2].get('periodic', None)

psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py

@@ -23,8 +23,9 @@
 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_scalar_conforming_projection, construct_vector_conforming_projection
 
 from psydac.linalg.utilities import array_to_psydac
 from psydac.fem.basic        import FemField
@@ -118,17 +119,15 @@ def solve_h1_source_pbm(
     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
+    H0_m  = H0.to_sparse_matrix()              # = mass matrix of V0
+    dH0_m = H0.get_dual_Hodge_sparse_matrix()  # = inverse mass matrix of V0
+    H1_m  = H1.to_sparse_matrix()              # = mass matrix of V1
+    dH1_m = H1.get_dual_Hodge_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 = 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_scalar_conforming_projection(V0h, hom_bc=[True,True])
+    # cP1_m = construct_vector_conforming_projection(V1h, domain_h, hom_bc=True)
 
     if not os.path.exists(plot_dir):
         os.makedirs(plot_dir)
@@ -138,7 +137,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()
@@ -150,20 +149,20 @@ def lift_u_bc(u_bc):
 
     # Conga (projection-based) stiffness matrices:
     # div grad:
-    pre_DG_m = - bD0_m.transpose() @ dH1_m @ bD0_m
+    pre_DG_m = - bD0_m.transpose() @ H1_m @ bD0_m
 
     # jump penalization:
     jump_penal_m = I0_m - cP0_m
-    JP0_m = jump_penal_m.transpose() * dH0_m * jump_penal_m
+    JP0_m = jump_penal_m.transpose() * H0_m * jump_penal_m
 
-    pre_A_m = cP0_m.transpose() @ ( eta * dH0_m - mu * pre_DG_m )  # useful for the boundary condition (if present)
+    pre_A_m = cP0_m.transpose() @ ( eta * H0_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, 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],
     )
@@ -173,26 +172,26 @@ 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)
+        b_c = H0_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, backend=PSYDAC_BACKENDS[backend_language])
+        lh = discretize(l, domain_h, V0h)
         b  = lh.assemble()
         b_c = b.toarray()
         if plot_source:
-            f_c = H0_m.dot(b_c)
+            f_c = dH0_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)
+        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)
 
@@ -265,6 +264,7 @@ def lift_u_bc(u_bc):
         mu=1, #1,
         domain_name=domain_name,
         source_type=source_type,
+        source_proj = 'P_geom',
         backend_language='pyccel-gcc',
         plot_source=True,
         plot_dir='./plots/h1_tests_source_february/'+run_dir,

psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py

@@ -18,10 +18,11 @@
 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_scalar_conforming_projection, construct_vector_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,
+def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language='python', mu=1, nu=0, gamma_h=10,
                           sigma=None, nb_eigs=4, nb_eigs_plot=4,
-                          plot_dir=None, hide_plots=True, m_load_dir="",):
+                          plot_dir=None, hide_plots=True, m_load_dir="",skip_eigs_threshold = 1e-7,):
     """
     solver for the eigenvalue problem: find lambda in R and u in H0(curl), such that
 
@@ -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
@@ -94,18 +95,17 @@ def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language
     H1 = HodgeOperator(V1h, domain_h, backend_language=backend_language, load_dir=m_load_dir, load_space_index=1)
     H2 = HodgeOperator(V2h, domain_h, backend_language=backend_language, load_dir=m_load_dir, load_space_index=2)
 
-    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
-    dH2_m = H2.get_dual_Hodge_sparse_matrix()  # = mass matrix of V2
+    H0_m  = H0.to_sparse_matrix()                # = mass matrix of V0
+    dH0_m = H0.get_dual_Hodge_sparse_matrix()    # = inverse mass matrix of V0
+    H1_m  = H1.to_sparse_matrix()                # = mass matrix of V1
+    dH1_m = H1.get_dual_Hodge_sparse_matrix()    # = inverse mass matrix of V1
+    H2_m = H2.to_sparse_matrix()                 # = mass matrix of V2
+    # dH2_m = H2.get_dual_Hodge_sparse_matrix()  # = inverse mass matrix of V2
 
     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_scalar_conforming_projection(V0h, hom_bc=[True, True])
+    cP1_m = construct_vector_conforming_projection(V1h, hom_bc=[True, True])
 
     print('broken differential operators...')
     bD0, bD1 = derham_h.broken_derivatives_as_operators
@@ -118,33 +118,59 @@ def hcurl_solve_eigen_pbm(nc=4, deg=4, domain_name='pretzel_f', backend_language
     # Conga (projection-based) stiffness matrices
     # curl curl:
     print('curl-curl stiffness matrix...')
-    pre_CC_m = bD1_m.transpose() @ dH2_m @ bD1_m
+    pre_CC_m = bD1_m.transpose() @ H2_m @ bD1_m
     CC_m = cP1_m.transpose() @ pre_CC_m @ cP1_m  # Conga stiffness matrix
 
     # grad div:
     print('grad-div stiffness matrix...')
-    pre_GD_m = - dH1_m @ bD0_m @ cP0_m @ H0_m @ cP0_m.transpose() @ bD0_m.transpose() @ dH1_m
+    pre_GD_m = - H1_m @ bD0_m @ cP0_m @ dH0_m @ cP0_m.transpose() @ bD0_m.transpose() @ H1_m
     GD_m = cP1_m.transpose() @ pre_GD_m @ cP1_m  # Conga stiffness matrix
 
     # jump penalization in V1h:
     jump_penal_m = I1_m - cP1_m
-    JP_m = jump_penal_m.transpose() * dH1_m * jump_penal_m
+    JP_m = jump_penal_m.transpose() * H1_m * jump_penal_m
 
     print('computing the full operator matrix...')
     print('mu = {}'.format(mu))
     print('nu = {}'.format(nu))
     A_m = mu * CC_m - nu * GD_m + gamma_h * JP_m
 
-    eigenvalues, eigenvectors = get_eigenvalues(nb_eigs, sigma, A_m, dH1_m)
+    if False: #gneralized problen
+        print('adding jump stabilization to RHS of generalized eigenproblem...')
+        B_m = cP1_m.transpose() @ H1_m @ cP1_m + JS_m
+    else:
+        B_m = H1_m
+
+    print('solving matrix eigenproblem...')
+    all_eigenvalues, all_eigenvectors_transp = get_eigenvalues(nb_eigs, sigma, A_m, B_m)
+    #Eigenvalue processing
+
+    zero_eigenvalues = []
+    if skip_eigs_threshold is not None:
+        eigenvalues = []
+        eigenvectors = []
+        for val, vect in zip(all_eigenvalues, all_eigenvectors_transp.T):
+            if abs(val) < skip_eigs_threshold: 
+                zero_eigenvalues.append(val)
+                # we skip the eigenvector
+            else:
+                eigenvalues.append(val)
+                eigenvectors.append(vect)
+    else:
+        eigenvalues = all_eigenvalues
+        eigenvectors = all_eigenvectors_transp.T
 
+
+
     # plot first eigenvalues
 
-    for i in range(min(nb_eigs_plot, nb_eigs)):
+    for i in range(min(nb_eigs_plot, len(eigenvalues))):
 
-        print('looking at emode i = {}... '.format(i))
         lambda_i  = eigenvalues[i]
-        emode_i = np.real(eigenvectors[:,i])
-        norm_emode_i = np.dot(emode_i,dH1_m.dot(emode_i))
+        print('looking at emode i = {}: {}... '.format(i, lambda_i))
+
+        emode_i = np.real(eigenvectors[i])
+        norm_emode_i = np.dot(emode_i,H1_m.dot(emode_i))
         print('norm of computed eigenmode: ', norm_emode_i)
         eh_c = emode_i/norm_emode_i  # numpy coeffs of the normalized eigenmode
         plot_field(numpy_coeffs=eh_c, Vh=V1h, space_kind='hcurl', domain=domain, title='mode e_{}, lambda_{}={}'.format(i,i,lambda_i),
@@ -203,8 +229,8 @@ def get_eigenvalues(nb_eigs, sigma, A_m, M_m):
 
     t_stamp_full = time_count()
 
-    quick_run = True
-    # quick_run = False
+    # quick_run = True
+    quick_run = False
 
     if quick_run:
         domain_name = 'curved_L_shape'
@@ -214,23 +240,32 @@ 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
+
+    sigma = 7
+    nb_eigs_solve = 7
+    nb_eigs_plot = 7
+    skip_eigs_threshold = 1e-7
 
     m_load_dir = 'matrices_{}_nc={}_deg={}/'.format(domain_name, nc, deg)
     run_dir = 'eigenpbm_{}_nc={}_deg={}/'.format(domain_name, nc, deg)
     hcurl_solve_eigen_pbm(
         nc=nc, deg=deg,
         nu=0,
         mu=1, #1,
-        sigma=1,
         domain_name=domain_name,
         backend_language='pyccel-gcc',
         plot_dir='./plots/tests_source_february/'+run_dir,
         hide_plots=True,
-        m_load_dir=m_load_dir,
+        m_load_dir=m_load_dir, 
+        gamma_h=0,
+        sigma=sigma, 
+        nb_eigs=nb_eigs_solve, 
+        nb_eigs_plot=nb_eigs_plot,
+        skip_eigs_threshold=skip_eigs_threshold,
     )
 
     time_count(t_stamp_full, msg='full program')