[WIP] Redefine dot and inner binary operators for real and complex fields

sympde/calculus/core.py

@@ -500,25 +500,27 @@ def __new__(cls, arg1, arg2, **options):
         else:
             b = [arg2]
 
-        args_1 = [i for i in a if not i.is_commutative]
+        args_1 = [i for i in a if not i.is_commutative or has(i, types)]
         c1     = [i for i in a if not i in args_1]
-        args_2 = [i for i in b if not i.is_commutative]
+        args_2 = [i for i in b if not i.is_commutative or has(i, types)]
         c2     = [i for i in b if not i in args_2]
 
-        a = reduce(mul, args_1)
-        b = reduce(mul, args_2)
-        c = Mul(*c1)*Mul(*c2)
+        a = reduce(mul, args_1, S.One)
+        b = reduce(mul, args_2, S.One)
+        c = Mul(*c1) * Mul(*c2)
 
+        # TODO [YG 20.05.2025] : apply conjugate if arguments are flipped
         if str(a) > str(b):
-            a,b = b,a
+            a, b = b, a
 
         obj = Basic.__new__(cls, a, b)
 
+        # This is a norm when a == b
         if a == b:
             obj.is_real     = True
             obj.is_positive = True
 
-        return c*obj
+        return c * obj
 
 #==============================================================================
 class Outer(BasicOperator):

sympde/calculus/tests/test_calculus.py

@@ -76,11 +76,14 @@ def test_Cross(dim):
     # TODO: check exceptions
 
 #==============================================================================
+@pytest.mark.parametrize('vector', [True, False])
 @pytest.mark.parametrize('dim', [1, 2, 3])
-def test_Inner(dim):
+def test_Inner(dim, vector):
+
+    FunctionSpace = VectorFunctionSpace if vector else ScalarFunctionSpace
 
     domain  = Domain('Omega', dim=dim)
-    W       = VectorFunctionSpace('W', domain)
+    W       = FunctionSpace('W', domain)
     a, b, c = elements_of(W, names='a, b, c')
     r       = Constant('r')
 

sympde/core/algebra.py

@@ -1,38 +1,17 @@
 # coding: utf-8
 
-import numpy as np
-from itertools import groupby
-
 #from sympy.core.sympify import sympify
 from sympy.simplify.simplify import simplify
-from sympy import Symbol
-from sympy import Lambda
-from sympy import Function
-from sympy import bspline_basis
-from sympy import lambdify
-from sympy import cos
-from sympy import sin
-from sympy import Rational
-from sympy import diff
 from sympy import Matrix, ImmutableDenseMatrix
-from sympy import latex
-from sympy import I as sympy_I
 from sympy.core import Basic
 from sympy.core.singleton import S
-from sympy.simplify.simplify import nsimplify
-from sympy.utilities.lambdify import implemented_function
-from sympy.matrices.dense import MutableDenseMatrix
 from sympy import Mul, Add
-from sympy import postorder_traversal
-from sympy import preorder_traversal
-from sympy.core.expr import Expr
 from sympy.core.containers import Tuple
-from sympy import Integer, Float
 from sympy import Add, Mul
 from sympy import simplify
 from sympy import S
 from sympy import Basic
-from sympy import Indexed, IndexedBase
+from sympy import Indexed
 
 from sympde.old_sympy_utilities import is_sequence
 from .basic import CalculusFunction
@@ -274,48 +253,37 @@ def __new__(cls, *args, **options):
         else:
             return r
 
-class Inner_2d(InnerBasic):
-
     @classmethod
     def eval(cls, *_args):
-        """."""
 
         if not _args:
             return
 
-        if not( len(_args) == 2):
+        if len(_args) != 2:
             raise ValueError('Expecting two arguments')
 
-        u = _args[0]
-        v = _args[1]
+        u, v = _args
+
+        # Matrix constructor requires a sequence
+        if not is_sequence(u): u = [u]
+        if not is_sequence(v): v = [v]
 
+        # Convert arguments to Matrix (nothing happens if already Matrix)
         u = Matrix(u)
         v = Matrix(v)
 
         # TODO add conjugate
-        M = u.transpose()*v
+        M = u.transpose() * v
         return M.trace()
 
+# ...
+class Inner_1d(InnerBasic):
+    pass
 
-class Inner_3d(InnerBasic):
-
-    @classmethod
-    def eval(cls, *_args):
-        """."""
-
-        if not _args:
-            return
-
-        if not( len(_args) == 2):
-            raise ValueError('Expecting two arguments')
-
-        u = _args[0]
-        v = _args[1]
-
-        u = Matrix(u)
-        v = Matrix(v)
+# ...
+class Inner_2d(InnerBasic):
+    pass
 
-        # TODO add conjugate
-        M = u.transpose()*v
-        return M.trace()
 # ...
+class Inner_3d(InnerBasic):
+    pass

sympde/expr/evaluation.py

@@ -1,13 +1,7 @@
 # coding: utf-8
 
-from itertools import product
-
-
-import numpy as np
-
 from sympy import Abs, S, cacheit
 from sympy import Indexed, Matrix, ImmutableDenseMatrix
-from sympy import expand
 from sympy.core import Basic, Symbol
 from sympy.core import Add, Mul, Pow
 from sympy.core.expr import AtomicExpr
@@ -16,22 +10,20 @@
 
 from sympde.core.basic import _coeffs_registery
 from sympde.core.basic import CalculusFunction
-from sympde.core.algebra import (Dot_1d,
+from sympde.core.algebra import (Dot_1d, Inner_1d,
                                  Dot_2d, Inner_2d, Cross_2d,
                                  Dot_3d, Inner_3d, Cross_3d)
-from sympde.core.utils import random_string
 
-from sympde.calculus import jump, avg, minus, plus
+from sympde.calculus import jump, minus, plus
 from sympde.calculus import Jump, is_zero
-from sympde.calculus.core import _generic_ops, _diff_ops
 from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose
-from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
+from sympde.calculus.matrices import MatrixElement, SymbolicTrace
 
 from sympde.topology.basic   import BasicDomain, Union, Interval
 from sympde.topology.basic   import Boundary, Interface
 from sympde.topology.basic   import InteriorDomain
 from sympde.topology.domain  import NormalVector, TangentVector, NCube, NCubeInterior
-from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, MultiPatchMapping, JacobianInverseSymbol
+from sympde.topology.mapping import JacobianSymbol, InterfaceMapping, JacobianInverseSymbol
 from sympde.topology.mapping import LogicalExpr, PullBack
 
 # TODO fix circular dependency between sympde.expr.evaluation and sympde.topology.mapping
@@ -40,7 +32,6 @@
 from sympde.topology.space import VectorFunction
 from sympde.topology.space import IndexedVectorFunction
 from sympde.topology.space import Trace
-from sympde.topology.space import element_of
 from sympde.topology.space import ScalarFunctionSpace
 
 from sympde.topology.derivatives import _partial_derivatives
@@ -49,15 +40,22 @@
 from sympde.topology.derivatives import get_atom_logical_derivatives
 from sympde.topology.derivatives import dx, dy, dz
 from sympde.topology.derivatives import dx1, dx2, dx3
-from sympde.topology.derivatives import (Grad_1d, Div_1d,
-                                         Grad_2d, Curl_2d, Rot_2d, Div_2d,
-                                         Grad_3d, Curl_3d, Div_3d)
+
+from sympde.calculus.core import Grad, Div, Curl, Rot, Bracket, Laplace, Hessian
+from sympde.calculus.core import Dot, Inner, Cross
+
+from sympde.topology.derivatives import Grad_1d, Grad_2d, Grad_3d
+from sympde.topology.derivatives import Div_1d, Div_2d, Div_3d
+from sympde.topology.derivatives import Curl_2d, Curl_3d
+from sympde.topology.derivatives import Rot_2d
 from sympde.topology.derivatives import Bracket_2d
 from sympde.topology.derivatives import Laplace_1d, Laplace_2d, Laplace_3d
 from sympde.topology.derivatives import Hessian_1d, Hessian_2d, Hessian_3d
-from sympde.topology.derivatives import (LogicalGrad_1d, LogicalDiv_1d,
-                                         LogicalGrad_2d, LogicalCurl_2d, LogicalRot_2d, LogicalDiv_2d,
-                                         LogicalGrad_3d, LogicalCurl_3d, LogicalDiv_3d)
+
+from sympde.topology.derivatives import LogicalGrad_1d, LogicalGrad_2d, LogicalGrad_3d
+from sympde.topology.derivatives import LogicalDiv_1d, LogicalDiv_2d, LogicalDiv_3d
+from sympde.topology.derivatives import LogicalCurl_2d, LogicalCurl_3d
+from sympde.topology.derivatives import LogicalRot_2d
 from sympde.topology.derivatives import LogicalBracket_2d
 from sympde.topology.derivatives import LogicalLaplace_1d, LogicalLaplace_2d, LogicalLaplace_3d
 from sympde.topology.derivatives import LogicalHessian_1d, LogicalHessian_2d, LogicalHessian_3d
@@ -91,6 +89,33 @@
     'is_sequence',
 )
 
+#==============================================================================
+differential_operators = {
+    Grad   : (None,    Grad_1d,    Grad_2d,    Grad_3d),
+    Div    : (None,     Div_1d,     Div_2d,     Div_3d),
+    Curl   : (None,       None,    Curl_2d,    Curl_3d),
+    Rot    : (None,       None,     Rot_2d,       None),
+    Bracket: (None,       None, Bracket_2d,       None),
+    Laplace: (None, Laplace_1d, Laplace_2d, Laplace_3d),
+    Hessian: (None, Hessian_1d, Hessian_2d, Hessian_3d),
+}
+
+logical_differential_operators = {
+    Grad   : (None,    LogicalGrad_1d,    LogicalGrad_2d,    LogicalGrad_3d),
+    Div    : (None,     LogicalDiv_1d,     LogicalDiv_2d,     LogicalDiv_3d),
+    Curl   : (None,              None,    LogicalCurl_2d,    LogicalCurl_3d),
+    Rot    : (None,              None,     LogicalRot_2d,              None),
+    Bracket: (None,              None, LogicalBracket_2d,              None),
+    Laplace: (None, LogicalLaplace_1d, LogicalLaplace_2d, LogicalLaplace_3d),
+    Hessian: (None, LogicalHessian_1d, LogicalHessian_2d, LogicalHessian_3d),
+}
+
+generic_operators = {
+    Dot  : (None,   Dot_1d,   Dot_2d,   Dot_3d),
+    Inner: (None, Inner_1d, Inner_2d, Inner_3d),
+    Cross: (None,     None, Cross_2d, Cross_3d),
+}
+
 #==============================================================================
 def is_sequence(a):
     return isinstance(a, (list,tuple,Tuple))
@@ -769,23 +794,35 @@ def eval(cls, expr, domain):
         elif isinstance(expr, BasicExpr):
             return cls.eval(expr.expr, domain=domain)
 
-        elif isinstance(expr, _diff_ops):
-            op   = type(expr)
+        elif isinstance(expr, (*differential_operators,)):
             if domain.mapping is None:
-                new  = eval('Logical{0}_{1}d'.format(op, dim))
+                new = logical_differential_operators[type(expr)][dim]
             else:
-                new  = eval('{0}_{1}d'.format(op, dim))
-
+                new = differential_operators[type(expr)][dim]
             args = [cls.eval(i, domain=domain) for i in expr.args]
             return new(*args)
 
-        elif isinstance(expr, _generic_ops):
-            # if i = Dot(...) then type(i) is Grad
-            op = type(expr)
-            new  = eval('{0}_{1}d'.format(op, dim))
+        elif isinstance(expr, (*generic_operators,)):
+            new = generic_operators[type(expr)][dim]
             args = [cls.eval(i, domain=domain) for i in expr.args]
             return new(*args)
 
+#        elif isinstance(expr, _diff_ops):
+#            op   = type(expr)
+#            if domain.mapping is None:
+#                new  = eval('Logical{0}_{1}d'.format(op, dim))
+#            else:
+#                new  = eval('{0}_{1}d'.format(op, dim))
+#            args = [cls.eval(i, domain=domain) for i in expr.args]
+#            return new(*args)
+#
+#        elif isinstance(expr, _generic_ops):
+#            # if i = Dot(...) then type(i) is Grad
+#            op = type(expr)
+#            new  = eval('{0}_{1}d'.format(op, dim))
+#            args = [cls.eval(i, domain=domain) for i in expr.args]
+#            return new(*args)
+
         elif isinstance(expr, Trace):
             # TODO treate different spaces
             if expr.order == 0: