Merge branch 'devel' into non_matching_multipatch

psydac/api/ast/nodes.py

@@ -9,7 +9,7 @@
 from sympy import Mul,Integer
 from sympy.core.singleton     import Singleton
 from sympy.core.containers    import Tuple
-from sympy.core.compatibility import with_metaclass
+from sympde.old_sympy_utilities import with_metaclass
 
 from sympde.topology import element_of
 from sympde.topology import ScalarFunction, VectorFunction

psydac/api/ast/parser.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 
+from sympy import S
 from sympy import IndexedBase, Indexed
 from sympy import Mul, Matrix, Expr
 from sympy import Add, And, StrictLessThan, Eq
@@ -253,9 +254,9 @@ def _visit_VectorAssign(self, expr, **kwargs):
         lhs = self._visit(expr.lhs)
         rhs = self._visit(expr.rhs)
         if expr.op is None:
-            return [Assign(l,r) for l,r in zip(lhs, rhs) if l is not None and r is not None]
+            return [Assign(l,r) for l,r in zip(lhs, rhs) if l is not None and r is not None and l is not S.Zero]
         else:
-            return [AugAssign(l,expr.op, r) for l,r in zip(lhs, rhs) if l is not None and r is not None]
+            return [AugAssign(l,expr.op, r) for l,r in zip(lhs, rhs) if l is not None and r is not None and l is not S.Zero]
     # ....................................................
     def _visit_Assign(self, expr, **kwargs):
 
@@ -1383,7 +1384,7 @@ def _visit_ElementOf(self, expr, **kwargs):
             targets = self._visit_BlockStencilMatrixLocalBasis(target)
             for i in range(targets.shape[0]):
                 for j in range(targets.shape[1]):
-                    if targets[i,j] is None:
+                    if targets[i,j] is S.Zero:
                         continue
                     if trials[j] in pads.trials_multiplicity:
                         trials_m  = pads.trials_multiplicity[trials[j]]
@@ -1414,7 +1415,7 @@ def _visit_ElementOf(self, expr, **kwargs):
             indices = list(rows)
             for i in range(targets.shape[0]):
                 for j in range(targets.shape[1]):
-                    if targets[i,j] is None:
+                    if targets[i,j] is S.Zero:
                         continue
                     targets[i,j] = targets[i,j][indices]
             return targets
@@ -1445,7 +1446,6 @@ def _visit_BlockStencilMatrixLocalBasis(self, expr, **kwargs):
         for i, v in enumerate(tests):
             for j, u in enumerate(trials):
                 if expr.expr[i, j] == 0:
-                    targets[i, j] = None
                     continue
                 mat = StencilMatrixLocalBasis(u=u, v=v, pads=pads[i, j], tag=tag, dtype=dtype)
                 mat = self._visit_StencilMatrixLocalBasis(mat, **kwargs)
@@ -1464,7 +1464,6 @@ def _visit_BlockStencilMatrixGlobalBasis(self, expr, **kwargs):
         for i, v in enumerate(tests):
             for j, u in enumerate(trials):
                 if expr.expr[i, j] == 0:
-                    targets[i, j] = None
                     continue
                 mat = StencilMatrixGlobalBasis(u=u, v=v, pads=pads, tag=tag, dtype=dtype)
                 mat = self._visit_StencilMatrixGlobalBasis(mat, **kwargs)
@@ -1480,7 +1479,6 @@ def _visit_BlockStencilVectorLocalBasis(self, expr, **kwargs):
         targets = Matrix.zeros(len(tests), 1)
         for i, v in enumerate(tests):
             if expr.expr[i, 0] == 0:
-                targets[i, 0] = None
                 continue
             mat = StencilVectorLocalBasis(v, pads, tag, dtype)
             mat = self._visit_StencilVectorLocalBasis(mat, **kwargs)
@@ -1496,7 +1494,6 @@ def _visit_BlockStencilVectorGlobalBasis(self, expr, **kwargs):
         targets = Matrix.zeros(len(tests), 1)
         for i,v in enumerate(tests):
             if expr.expr[i, 0] == 0:
-                targets[i, 0] = None
                 continue
             mat = StencilVectorGlobalBasis(v, pads, tag, dtype)
             mat = self._visit_StencilVectorGlobalBasis(mat, **kwargs)
@@ -1513,7 +1510,6 @@ def _visit_BlockScalarLocalBasis(self, expr, **kwargs):
         for i,v in enumerate(tests):
             for j,u in enumerate(trials):
                 if expr.expr[i, j] == 0:
-                    targets[i, j] = None
                     continue
                 var = ScalarLocalBasis(u, v, tag, dtype)
                 var = self._visit_ScalarLocalBasis(var, **kwargs)

psydac/api/ast/tests/poisson.py

@@ -1,27 +1,25 @@
 # -*- coding: UTF-8 -*-
-
+import os
 
 from sympy import symbols
-from sympy import sin,pi
+from sympy import sin, pi
 
 from sympde.calculus import grad, dot
-
 from sympde.topology import ScalarFunctionSpace
-from sympde.topology import elements_of,LogicalExpr
+from sympde.topology import elements_of, LogicalExpr
 from sympde.topology import Square
-from sympde.topology import IdentityMapping,Mapping ,PolarMapping
+from sympde.topology import Mapping, IdentityMapping, PolarMapping
 from sympde.expr     import integral
 from sympde.expr     import LinearForm
 from sympde.expr     import BilinearForm
-from sympde.expr     import Norm
+from sympde.expr     import Norm, SemiNorm
 from sympde.expr.evaluation import TerminalExpr
 
 from psydac.api.ast.fem          import AST
 from psydac.api.ast.parser       import parse
 from psydac.api.discretization   import discretize
 from psydac.api.printing.pycode  import pycode
 
-import os
 # ...
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
@@ -32,6 +30,7 @@
     mesh_dir = os.path.join(base_dir, 'mesh')
     filename = os.path.join(mesh_dir, 'identity_2d.h5')
 
+
 def test_codegen():
     domain  = Square()
     M       = Mapping('M',2)
@@ -51,22 +50,22 @@ def test_codegen():
     Vh = discretize(V, domain_h)
 
     error  = u - sin(pi*x)*sin(pi*y)
-    l2norm = LogicalExpr(M, Norm(error, domain, kind='l2'))
-    h1norm = LogicalExpr(M, Norm(error, domain, kind='h1'))
+    l2norm = LogicalExpr(M,     Norm(error, domain, kind='l2'))
+    h1norm = LogicalExpr(M, SemiNorm(error, domain, kind='h1'))
 
     ast_b = AST(b, TerminalExpr(b)[0],[Vh, Vh])
-    ast_b = parse(ast_b.expr, settings={'dim':2,'nderiv':1,'mapping':M,'target':domain.logical_domain})
+    ast_b = parse(ast_b.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain})
     print(pycode(ast_b))
 
     print('==============================================================================================================')
     ast_l = AST(l, TerminalExpr(l)[0], Vh)
-    ast_l = parse(ast_l.expr, settings={'dim':2,'nderiv':1,'mapping':M,'target':domain.logical_domain})
+    ast_l = parse(ast_l.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain})
     print(pycode(ast_l))
 
 
     print('==============================================================================================================')
     ast_l2norm = AST(l2norm, TerminalExpr(l2norm)[0], Vh)
-    ast_l2norm = parse(ast_l2norm.expr, settings={'dim':2,'nderiv':1,'mapping':M, 'target':domain.logical_domain})
+    ast_l2norm = parse(ast_l2norm.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain})
     print(pycode(ast_l2norm))
-test_codegen()
 
+test_codegen()

psydac/api/ast/tests/system_1.py

@@ -1,4 +1,5 @@
 # -*- coding: UTF-8 -*-
+import os
 
 from sympy import pi, sin, Tuple, Matrix
 
@@ -10,8 +11,7 @@
 from sympde.expr     import integral
 from sympde.expr     import LinearForm
 from sympde.expr     import BilinearForm
-from sympde.expr     import Norm
-
+from sympde.expr     import Norm, SemiNorm
 from sympde.expr.evaluation import TerminalExpr
 
 from psydac.api.ast.fem          import AST
@@ -20,7 +20,6 @@
 from psydac.api.printing.pycode  import pycode
 
 
-import os
 # ...
 try:
     mesh_dir = os.environ['PSYDAC_MESH_DIR']
@@ -39,7 +38,7 @@ def test_codegen():
 
     x,y = domain.coordinates
 
-    f   = Tuple(2*pi**2*sin(pi*x)*sin(pi*y),
+    f = Tuple(2*pi**2*sin(pi*x)*sin(pi*y),
               2*pi**2*sin(pi*x)*sin(pi*y))
 
     Fe = Tuple(sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y))
@@ -50,12 +49,12 @@ def test_codegen():
 
     int_0 = lambda expr: integral(domain , expr)
 
-    b = BilinearForm((u,v), int_0(inner(grad(v), grad(u))))
+    b = BilinearForm((u, v), int_0(inner(grad(v), grad(u))))
     l = LinearForm(v, int_0(dot(f, v)))
 
     error = Matrix([F[0]-Fe[0], F[1]-Fe[1]])
-    l2norm_F = Norm(error, domain, kind='l2')
-    h1norm_F = Norm(error, domain, kind='h1')
+    l2norm_F =     Norm(error, domain, kind='l2')
+    h1norm_F = SemiNorm(error, domain, kind='h1')
 
     # Create computational domain from topological domain
     domain_h = discretize(domain, filename=filename)
@@ -64,17 +63,16 @@ def test_codegen():
     Vh = discretize(V, domain_h)
 
     print('============================================BilinearForm=========================================')
-    ast_b    = AST(b, TerminalExpr(b)[0], [Vh, Vh])
-    stmt_b = parse(ast_b.expr, settings={'dim':2,'nderiv':1, 'mapping':Vh.symbolic_mapping})
+    ast_b  = AST(b, TerminalExpr(b)[0], [Vh, Vh])
+    stmt_b = parse(ast_b.expr, settings={'dim':2, 'nderiv':1, 'mapping':Vh.symbolic_mapping})
     print(pycode(stmt_b))
 
     print('============================================LinearForm===========================================')
-    ast_l    = AST(l, TerminalExpr(l)[0], Vh)
-    stmt_l = parse(ast_l.expr, settings={'dim':2,'nderiv':1, 'mapping':Vh.symbolic_mapping})
+    ast_l  = AST(l, TerminalExpr(l)[0], Vh)
+    stmt_l = parse(ast_l.expr, settings={'dim':2, 'nderiv':1, 'mapping':Vh.symbolic_mapping})
     print(pycode(stmt_l))
 
-    print('============================================Norm===========================================')
+    print('============================================SemiNorm===========================================')
     ast_norm = AST(h1norm_F, TerminalExpr(h1norm_F)[0], Vh)
-    stmt_n = parse(ast_norm.expr, settings={'dim':2,'nderiv':1, 'mapping':Vh.symbolic_mapping})
+    stmt_n = parse(ast_norm.expr, settings={'dim':2, 'nderiv':1, 'mapping':Vh.symbolic_mapping})
     print(pycode(stmt_n))
-

psydac/api/ast/tests/test_nodes.py

@@ -579,11 +579,11 @@ def assembly_bilinear_form_2d_1():
 #==============================================================================
 
 def teardown_module():
-    from sympy import cache
+    from sympy.core import cache
     cache.clear_cache()
 
 def teardown_function():
-    from sympy import cache
+    from sympy.core import cache
     cache.clear_cache()
 
 

psydac/api/ast/utilities.py

@@ -891,11 +891,89 @@ def build_pythran_types_header(name, args, order=None):
 from sympy import preorder_traversal
 from sympy import NumberSymbol
 from sympy import Pow, S
-from sympy.printing.pycode import _known_functions_math
-from sympy.printing.pycode import _known_constants_math
-from sympy.printing.pycode import _known_functions_mpmath
-from sympy.printing.pycode import _known_constants_mpmath
-from sympy.printing.pycode import _known_functions_numpy
+
+_known_functions_math = {
+    'acos': 'acos',
+    'acosh': 'acosh',
+    'asin': 'asin',
+    'asinh': 'asinh',
+    'atan': 'atan',
+    'atan2': 'atan2',
+    'atanh': 'atanh',
+    'ceiling': 'ceil',
+    'cos': 'cos',
+    'cosh': 'cosh',
+    'erf': 'erf',
+    'erfc': 'erfc',
+    'exp': 'exp',
+    'expm1': 'expm1',
+    'factorial': 'factorial',
+    'floor': 'floor',
+    'gamma': 'gamma',
+    'hypot': 'hypot',
+    'loggamma': 'lgamma',
+    'log': 'log',
+    'ln': 'log',
+    'log10': 'log10',
+    'log1p': 'log1p',
+    'log2': 'log2',
+    'sin': 'sin',
+    'sinh': 'sinh',
+    'Sqrt': 'sqrt',
+    'tan': 'tan',
+    'tanh': 'tanh'
+
+}  # Not used from ``math``: [copysign isclose isfinite isinf isnan ldexp frexp pow modf
+# radians trunc fmod fsum gcd degrees fabs]
+_known_constants_math = {
+    'Exp1': 'e',
+    'Pi': 'pi',
+    'E': 'e'
+    # Only in python >= 3.5:
+    # 'Infinity': 'inf',
+    # 'NaN': 'nan'
+}
+
+_not_in_mpmath = 'log1p log2'.split()
+_in_mpmath = [(k, v) for k, v in _known_functions_math.items() if k not in _not_in_mpmath]
+_known_functions_mpmath = dict(_in_mpmath, **{
+    'beta': 'beta',
+    'fresnelc': 'fresnelc',
+    'fresnels': 'fresnels',
+    'sign': 'sign',
+})
+_known_constants_mpmath = {
+    'Exp1': 'e',
+    'Pi': 'pi',
+    'GoldenRatio': 'phi',
+    'EulerGamma': 'euler',
+    'Catalan': 'catalan',
+    'NaN': 'nan',
+    'Infinity': 'inf',
+    'NegativeInfinity': 'ninf'
+}
+
+_not_in_numpy = 'erf erfc factorial gamma loggamma'.split()
+_in_numpy = [(k, v) for k, v in _known_functions_math.items() if k not in _not_in_numpy]
+_known_functions_numpy = dict(_in_numpy, **{
+    'acos': 'arccos',
+    'acosh': 'arccosh',
+    'asin': 'arcsin',
+    'asinh': 'arcsinh',
+    'atan': 'arctan',
+    'atan2': 'arctan2',
+    'atanh': 'arctanh',
+    'exp2': 'exp2',
+    'sign': 'sign',
+})
+_known_constants_numpy = {
+    'Exp1': 'e',
+    'Pi': 'pi',
+    'EulerGamma': 'euler_gamma',
+    'NaN': 'nan',
+    'Infinity': 'PINF',
+    'NegativeInfinity': 'NINF'
+}
 
 
 def math_atoms_as_str(expr, lib='math'):
@@ -929,7 +1007,7 @@ def math_atoms_as_str(expr, lib='math'):
         known_constants = _known_constants_mpmath
     elif lib == 'numpy':
         known_functions = _known_functions_numpy
-        known_constants = _known_constants_math   # numpy version missing
+        known_constants = _known_constants_numpy   # numpy version missing
     else:
         raise ValueError("Library {} not supported.".format(mod))
 

psydac/api/discretization.py

@@ -3,6 +3,7 @@
 # TODO: - init_fem is called whenever we call discretize. we should check that
 #         nderiv has not been changed. shall we add nquads too?
 import os
+
 from sympy import Expr as sym_Expr
 import numpy as np
 
@@ -11,7 +12,7 @@
 from sympde.expr     import LinearForm as sym_LinearForm
 from sympde.expr     import Functional as sym_Functional
 from sympde.expr     import Equation as sym_Equation
-from sympde.expr     import Norm as sym_Norm
+from sympde.expr     import Norm as sym_Norm, SemiNorm as sym_SemiNorm
 from sympde.expr     import TerminalExpr
 
 from sympde.topology import BasicFunctionSpace
@@ -39,9 +40,8 @@
 from psydac.fem.vector       import ProductFemSpace, VectorFemSpace
 from psydac.cad.geometry     import Geometry
 from psydac.mapping.discrete import NurbsMapping
-
-from psydac.linalg.stencil     import StencilVectorSpace
-from psydac.linalg.block       import BlockVectorSpace
+from psydac.linalg.stencil   import StencilVectorSpace
+from psydac.linalg.block     import BlockVectorSpace
 
 __all__ = ('discretize', 'discretize_derham', 'reduce_space_degrees', 'discretize_space', 'discretize_domain')
 
@@ -445,7 +445,7 @@ def discretize(a, *args, **kwargs):
         kwargs['symbolic_mapping'] = mapping
 
     if isinstance(a, sym_BasicForm):
-        if isinstance(a, sym_Norm):
+        if isinstance(a, (sym_Norm, sym_SemiNorm)):
             kernel_expr = TerminalExpr(a, domain)
             if not mapping is None:
                 kernel_expr = tuple(LogicalExpr(i, domain) for i in kernel_expr)