Remove obsolete code generation for periodic case.

This is not needed anymore, because the code that is generated for the
non-periodic case also works for a periodic domain. Moreover, the old
special code was inefficient, not well tested, and it was not working
if the spline degree was not the same along all directions.

To verify that the periodic case is treated properly, a unit test is added:
2D Poisson's equation on a unit circle, with periodic BCs along the theta
direction, homogeneous Dirichlet BC at r=1, and different spline degrees
along r and theta.

psydac/api/ast/fem.py

@@ -327,10 +327,6 @@ def fields(self):
     def fields_coeffs(self):
         return self._fields_coeffs
 
-    @property
-    def fields_tmp_coeffs(self):
-        return self._fields_tmp_coeffs
-
     @property
     def vector_fields(self):
         return self._vector_fields
@@ -679,8 +675,6 @@ def _initialize(self):
                                           dtype='real', rank=dim, cls=IndexedVariable)
         fields_val    = variables(['{}_values'.format(f) for f in fields_logical_str],
                                           dtype='real', rank=dim, cls=IndexedVariable)
-        fields_tmp_coeffs = variables(['tmp_coeff_{}'.format(f) for f in field_atoms],
-                                              dtype='real', rank=dim, cls=IndexedVariable)
 
         vector_fields        = symbols(vector_fields_str)
         vector_fields_logical = symbols(vector_fields_logical_str)
@@ -747,7 +741,6 @@ def _initialize(self):
         self._fields = fields
         self._fields_logical = fields_logical
         self._fields_coeffs = fields_coeffs
-        self._fields_tmp_coeffs = fields_tmp_coeffs
         self._vector_fields = vector_fields
         self._vector_fields_logical = vector_fields_logical
         self._vector_fields_coeffs = vector_fields_coeffs
@@ -1065,43 +1058,28 @@ def _initialize(self):
 #==============================================================================
 class Assembly(SplBasic):
 
-    def __new__(cls, kernel, name=None, discrete_space=None, periodic=None,
-                comm=None, mapping=None, is_rational_mapping=None, backend=None):
+    def __new__(cls, kernel, name=None, discrete_space=None, comm=None,
+                mapping=None, is_rational_mapping=None, backend=None):
 
+        # ... Check arguments
         if not isinstance(kernel, Kernel):
             raise TypeError('> Expecting a kernel')
 
+        if isinstance(discrete_space, (tuple, list)):
+            space = discrete_space[0]
+        else:
+            space = discrete_space
+
+        if not isinstance(space, (SplineSpace, TensorFemSpace, ProductFemSpace)):
+            raise NotImplementedError('Only Spline, Tensor and Product spaces are available')
+        # ...
+
         obj = SplBasic.__new__(cls, kernel.tag, name=name,
                                prefix='assembly', mapping=mapping,
                                is_rational_mapping=is_rational_mapping)
 
-
-        # ... get periodicity of the space
-        dim    = kernel.weak_form.ldim
-        periodic = [False for i in range(0, dim)]
-        if not( discrete_space is None ):
-            if isinstance(discrete_space, (tuple, list)):
-                space = discrete_space[0]
-
-            else:
-                space = discrete_space
-
-            if isinstance(space, SplineSpace):
-                periodic = [space.periodic]
-
-            elif isinstance(space, TensorFemSpace):
-                periodic = space.periodic
-
-            elif isinstance(space, ProductFemSpace):
-                periodic = space.spaces[0].periodic
-
-            else:
-                raise NotImplementedError('Only Spline, Tensor and Product spaces are available')
-        # ...
-
         obj._kernel = kernel
         obj._discrete_space = discrete_space
-        obj._periodic = periodic
         obj._comm = comm
         obj._discrete_boundary = kernel.discrete_boundary
         obj._backend = backend
@@ -1124,10 +1102,6 @@ def kernel(self):
     def discrete_space(self):
         return self._discrete_space
 
-    @property
-    def periodic(self):
-        return self._periodic
-
     @property
     def comm(self):
         return self._comm
@@ -1157,7 +1131,6 @@ def _initialize(self):
         form   = self.weak_form
         fields = kernel.fields
         fields_coeffs = kernel.fields_coeffs
-        fields_tmp_coeffs = kernel.fields_tmp_coeffs
         vector_fields = kernel.vector_fields
         vector_fields_coeffs = kernel.vector_fields_coeffs
         zero_terms = kernel.zero_terms
@@ -1326,56 +1299,11 @@ def _initialize(self):
         # ... kernel call
         mats = tuple(element_matrices.values())
 
-        if not( self.comm is None ) and any(self.periodic) :
-            # ...
-            stmts = []
-            for i,il,p,span,n,per in zip( indices_i,
-                                          indices_il,
-                                          test_degrees,
-                                          indices_span,
-                                          npts,
-                                          self.periodic ):
-
-                if not per:
-                    stmts += [Assign(i, span-p+il)]
-
-                else:
-                    stmts += [Assign(i, Mod(span-p+il, n))]
-
-#            _indices_i = [i for i,s,p in zip(indices_i, starts, test_degrees)]
-            _indices_i = [i-s+p for i,s,p in zip(indices_i, starts, test_degrees)]
-            for x,tmp_x in zip(fields_coeffs, fields_tmp_coeffs):
-#                stmts += [Print([_indices_i, '    ', indices_i, starts])]
-                stmts += [Assign(tmp_x[indices_il], x[_indices_i])]
-
-            ranges = [Range(0, test_degrees[i]+1) for i in range(dim)]
-            for x,rx in list(zip(indices_il, ranges))[::-1]:
-                stmts = [For(x, rx, stmts)]
-
-            body += stmts
-            # ...
-
-            # ...
-            f_coeffs = tuple(fields_tmp_coeffs)
-            # ...
-
-            # ... TODO add tmp for vector fields and mapping
-            gslices = [Slice(sp-s-d+p,sp+p+1-s) for sp,d,p,s in zip(indices_span[::ln],
-                                                           test_degrees[::ln],
-                                                           test_pads[::ln],
-                                                           starts)]
-            vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs])
-            m_coeffs = tuple([f[gslices] for f in kernel.mapping_coeffs])
-            # ...
-
-        else:
-            gslices = [Slice(sp-s-d+p,sp+p+1-s) for sp,d,p,s in zip(indices_span[::ln],
-                                                           test_degrees[::ln],
-                                                           test_pads[::ln],
-                                                           starts)]
-            f_coeffs = tuple([f[gslices] for f in fields_coeffs])
-            vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs])
-            m_coeffs = tuple([f[gslices] for f in kernel.mapping_coeffs])
+        gslices = [Slice(sp-s-d+p, sp+p+1-s) for sp,d,p,s in
+                   zip(indices_span[::ln], test_degrees[::ln], test_pads[::ln], starts)]
+        f_coeffs  = tuple([f[gslices] for f in fields_coeffs])
+        vf_coeffs = tuple([f[gslices] for f in vector_fields_coeffs])
+        m_coeffs  = tuple([f[gslices] for f in kernel.mapping_coeffs])
 
         if is_bilinear:
             if not unique_scalar_space:
@@ -1543,12 +1471,6 @@ def _initialize(self):
             fields_shape = tuple(FunctionCall('len',[p[0,Slice(None,None)]]) for p in points)
             for F_value in self.kernel.vector_fields_val:
                 prelude += [Assign(F_value, Zeros(fields_shape))]
-        # ...
-        if not( self.comm is None ) and any(self.periodic) :
-            for v in self.kernel.fields_tmp_coeffs:
-                stmt = Assign(v, Zeros(orders))
-                prelude += [stmt]
-        # ...
 
         # ...
         if self.debug:

psydac/api/tests/test_api_2d_scalar_mapping.py

@@ -30,13 +30,20 @@
 # ...
 
 #==============================================================================
-def run_poisson_2d_dir(filename, solution, f, comm=None):
+def run_poisson_2d_dir(filename, solution, f, dir_boundary=None, comm=None):
 
     # ... abstract model
     domain = Domain.from_file(filename)
 
     V = ScalarFunctionSpace('V', domain)
 
+    if dir_boundary is None:
+        B_dirichlet = domain.boundary
+    elif len(dir_boundary) == 1:
+        B_dirichlet = domain.get_boundary(**dir_boundary[0])
+    else:
+        B_dirichlet = Union(*[domain.get_boundary(**kw) for kw in dir_boundary])
+
     x,y = domain.coordinates
 
     F = element_of(V, name='F')
@@ -56,7 +63,7 @@ def run_poisson_2d_dir(filename, solution, f, comm=None):
     l2norm = Norm(error, domain, kind='l2')
     h1norm = Norm(error, domain, kind='h1')
 
-    bc = EssentialBC(u, 0, domain.boundary)
+    bc = EssentialBC(u, 0, B_dirichlet)
     equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc)
     # ...
 
@@ -353,7 +360,6 @@ def test_api_poisson_2d_dir_collela():
     assert( abs(h1_error - expected_h1_error) < 1.e-7)
 
 #==============================================================================
-# TODO h1 norm not working => regression due to the last changes in sympde
 def test_api_poisson_2d_dir_quart_circle():
     filename = os.path.join(mesh_dir, 'quart_circle.h5')
 
@@ -370,7 +376,25 @@ def test_api_poisson_2d_dir_quart_circle():
     expected_h1_error =  0.009473407914765117
 
     assert( abs(l2_error - expected_l2_error) < 1.e-7)
-#    assert( abs(h1_error - expected_h1_error) < 1.e-7)
+    assert( abs(h1_error - expected_h1_error) < 1.e-7)
+
+#==============================================================================
+def test_api_poisson_2d_dir_circle():
+    filename = os.path.join(mesh_dir, 'circle.h5')
+
+    from sympy.abc import x,y
+
+    solution = (1 - (x**2 + y**2)) * cos(2*pi*x) * cos(2*pi*y)
+    f        = -laplace(solution)
+
+    dir_boundary = [{'axis': 0, 'ext': 1}]
+    l2_error, h1_error = run_poisson_2d_dir(filename, solution, f, dir_boundary)
+
+    expected_l2_error = 0.0015245737751297718
+    expected_h1_error = 0.06653900724243668
+
+    assert( abs(l2_error - expected_l2_error) < 1.e-7)
+    assert( abs(h1_error - expected_h1_error) < 1.e-7)
 
 #==============================================================================
 def test_api_poisson_2d_dirneu_identity_1():