Non Matching Multipatch (#320)

## List of changes - Add non matching conforming projection operators w/ homogeneous BCs for V0 and V1 - Update multipatch examples to newer Psydac version - Add `__init__.py` to tests - Add non matching examples from previous branch - Add more tests - Change in naming convention: primal Hodge operator corresponds to mass matrix --------- Co-authored-by: Yaman Güçlü <[email protected]> Co-authored-by: Martin Campos Pinto <[email protected]> Co-authored-by: jowezarek <[email protected]> Co-authored-by: e-moral-sanchez <[email protected]> Co-authored-by: Elena Moral Sánchez <[email protected]> Co-authored-by: kvrigor <[email protected]>
pyccel · Aug 2, 2024 · d3651dd · d3651dd
1 parent 655cb8a
commit d3651dd
Show file tree

Hide file tree

Showing 28 changed files with 7,108 additions and 1,762 deletions.
diff --git a/.gitignore b/.gitignore
@@ -5,6 +5,7 @@
 *.lock
 __psydac__/
 __*pyccel__/
+docs/source/modules/STUBDIR/*
 
 build
 *build*

diff --git a/docs/source/modules/feec.multipatch.rst b/docs/source/modules/feec.multipatch.rst
@@ -9,5 +9,9 @@ feec.multipatch
 
     multipatch.api
     multipatch.fem_linear_operators
+    multipatch.multipatch_domain_utilities
+    multipatch.non_matching_operators
     multipatch.operators
+    multipatch.plotting_utilities
     multipatch.utilities
+    multipatch.utils_conga_2d
diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py
@@ -102,7 +102,7 @@ def get_max_degree_of_one_space(Vh):
     Vh : FemSpace
         The finite element space under investigation.
 
-    Results
+    Returns
     -------
     list[int]
         The maximum polynomial degre of Vh with respect to each coordinate.
@@ -133,7 +133,7 @@ def get_max_degree(*spaces):
     *spaces : tuple[FemSpace]
         The finite element spaces under investigation.
 
-    Results
+    Returns
     -------
     list[int]
         The maximum polynomial degree across all spaces, with respect to each

diff --git a/psydac/api/fem.py b/psydac/api/fem.py
@@ -767,7 +767,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

diff --git a/psydac/api/postprocessing.py b/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)

diff --git a/psydac/api/settings.py b/psydac/api/settings.py
@@ -67,4 +67,4 @@
     'pyccel-intel' : PSYDAC_BACKEND_IPYCCEL,
     'pyccel-pgi'   : PSYDAC_BACKEND_PGPYCCEL,
     'pyccel-nvidia': PSYDAC_BACKEND_NVPYCCEL,
-}
+}
diff --git a/psydac/api/tests/build_domain.py b/psydac/api/tests/build_domain.py
@@ -1,10 +1,16 @@
 # coding: utf-8
 
+# todo: this file has a lot of redundant code with psydac/feec/multipatch/multipatch_domain_utilities.py
+# it should probably be removed in the future
+
 import numpy as np
 
 from sympde.topology import Square, Domain
 from sympde.topology import IdentityMapping, PolarMapping, AffineMapping, Mapping
 
+# remove after sympde PR #155 is merged and call Domain.join instead
+from psydac.feec.multipatch.multipatch_domain_utilities import sympde_Domain_join
+
 #==============================================================================
 # small extension to SymPDE:
 class TransposedPolarMapping(Mapping):
@@ -20,6 +26,7 @@ class TransposedPolarMapping(Mapping):
     _ldim        = 2
     _pdim        = 2
 
+# todo: remove this
 def create_domain(patches, interfaces, name):
     connectivity = []
     patches_interiors = [D.interior for D in patches]
@@ -54,6 +61,8 @@ def flip_axis(name='no_name', c1=0., c2=0.):
     )
 
 #==============================================================================
+
+# todo: use build_multipatch_domain instead
 def build_pretzel(domain_name='pretzel', r_min=None, r_max=None):
     """
     design pretzel-like domain
@@ -189,23 +198,29 @@ def build_pretzel(domain_name='pretzel', r_min=None, r_max=None):
                     domain_14,
                     ])
 
-    interfaces = [
-        [domain_1.get_boundary(axis=1, ext=+1), domain_5.get_boundary(axis=1, ext=-1),  1],
-        [domain_5.get_boundary(axis=1, ext=+1), domain_6.get_boundary(axis=1, ext=1),   1],
-        [domain_6.get_boundary(axis=1, ext=-1), domain_2.get_boundary(axis=1, ext=-1),  1],
-        [domain_2.get_boundary(axis=1, ext=+1), domain_7.get_boundary(axis=1, ext=-1),  1],
-        [domain_7.get_boundary(axis=1, ext=+1), domain_3.get_boundary(axis=1, ext=-1),  1],
-        [domain_3.get_boundary(axis=1, ext=+1), domain_9.get_boundary(axis=1, ext=-1),  1],
-        [domain_9.get_boundary(axis=1, ext=+1), domain_4.get_boundary(axis=1, ext=-1),  1],
-        [domain_4.get_boundary(axis=1, ext=+1), domain_12.get_boundary(axis=1, ext=1),  1],
-        [domain_12.get_boundary(axis=1, ext=-1), domain_1.get_boundary(axis=1, ext=-1), 1],
-        [domain_6.get_boundary(axis=0, ext=-1), domain_13.get_boundary(axis=0, ext=1),  1],
-        [domain_7.get_boundary(axis=0, ext=-1), domain_13.get_boundary(axis=0, ext=-1), 1],
-        [domain_5.get_boundary(axis=0, ext=-1), domain_14.get_boundary(axis=0, ext=-1), 1],
-        [domain_12.get_boundary(axis=0, ext=-1), domain_14.get_boundary(axis=0, ext=+1),1],
+    axis_0 = 0
+    axis_1 = 1
+    ext_0 = -1
+    ext_1 = +1
+
+    connectivity = [
+        [(domain_1,  axis_1, ext_1), (domain_5,  axis_1, ext_0), 1],
+        [(domain_5,  axis_1, ext_1), (domain_6,  axis_1, ext_1), 1],
+        [(domain_6,  axis_1, ext_0), (domain_2,  axis_1, ext_0), 1],
+        [(domain_2,  axis_1, ext_1), (domain_7,  axis_1, ext_0), 1],
+        [(domain_7,  axis_1, ext_1), (domain_3,  axis_1, ext_0), 1],
+        [(domain_3,  axis_1, ext_1), (domain_9,  axis_1, ext_0), 1],
+        [(domain_9,  axis_1, ext_1), (domain_4,  axis_1, ext_0), 1],
+        [(domain_4,  axis_1, ext_1), (domain_12, axis_1, ext_1), 1],
+        [(domain_12, axis_1, ext_0), (domain_1,  axis_1, ext_0), 1],
+        [(domain_6,  axis_0, ext_0), (domain_13, axis_0, ext_1), 1],
+        [(domain_7,  axis_0, ext_0), (domain_13, axis_0, ext_0), 1],
+        [(domain_5,  axis_0, ext_0), (domain_14, axis_0, ext_0), 1],
+        [(domain_12, axis_0, ext_0), (domain_14, axis_0, ext_1), 1],
         ]
 
+    # domain = Domain.join(patches, connectivity, name=domain_name)
+    domain = sympde_Domain_join(patches, connectivity, name=domain_name)
 
-    domain = create_domain(patches, interfaces, domain_name)
     return domain