Merge pull request #56 from pyccel/devel-periodic

Remove obsolete code generation for periodic case.

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():