Add new LazyBinaryOperator pow

src/LazyOperators/LazyOperators.jl

@@ -1,7 +1,27 @@
 module LazyOperators
 
 import Base:
-    *, /, +, -, max, min, sqrt, abs, tan, sin, cos, tanh, sinh, cosh, atan, asin, acos, zero
+    *,
+    /,
+    +,
+    -,
+    ^,
+    max,
+    min,
+    sqrt,
+    abs,
+    tan,
+    sin,
+    cos,
+    tanh,
+    sinh,
+    cosh,
+    atan,
+    asin,
+    acos,
+    zero,
+    one,
+    materialize
 
 using LinearAlgebra
 import LinearAlgebra: dot, transpose, tr

src/LazyOperators/algebra.jl

@@ -2,7 +2,7 @@
 # Define rules for binary operators
 ###############################################################
 
-const LazyBinaryOp = (:*, :/, :+, :-, :max, :min, :dot)
+const LazyBinaryOp = (:*, :/, :+, :-, :^, :max, :min, :dot)
 
 for f in LazyBinaryOp
     @eval ($f)(a::AbstractLazy, b::AbstractLazy) = LazyOperator($f, a, b)
@@ -31,6 +31,7 @@
     :asin,
     :acos,
     :zero,
+    :one,
 )
 
 for f in LazyUnaryOp
@@ -65,6 +66,13 @@
 Base.:/(a::NullOperator, b) = a
 Base.:/(::NullOperator, ::NullOperator) = NullOperator()
 
+# For binary `^`:
+Base.:^(a, ::NullOperator) = error("Undefined")
+Base.:^(a::NullOperator, b) = a
+Base.:^(::NullOperator, ::NullOperator) = NullOperator() # or "I" ?
+Base.:^(a::AbstractLazy, ::NullOperator) = error("Undefined")
+Base.:^(a::NullOperator, b::AbstractLazy) = a
+
 ###############################################################
 # Define rules with `broadcasted`
 ###############################################################

src/LazyOperators/mapover.jl

@@ -102,6 +102,8 @@
 Base.:+(::NullOperator, b::AbstractMapOver) = b
 Base.:-(a::AbstractMapOver, ::NullOperator) = a
 Base.:-(::NullOperator, b::AbstractMapOver) = -b
+Base.:^(::AbstractMapOver, ::NullOperator) = error("Undefined")
+Base.:^(::NullOperator, ::AbstractMapOver) = NullOperator()
 Base.max(a::AbstractMapOver, ::NullOperator) = a
 Base.max(::NullOperator, b::AbstractMapOver) = b
 Base.min(a::AbstractMapOver, ::NullOperator) = a

test/operator/test_algebra.jl

@@ -195,4 +195,19 @@
         a = assemble_linear(l, V)
         @test all(a .≈ [-3.0, 3.0, -3.0, 3.0, -3.0, -3.0, 3.0, 3.0])
     end
+
+    @testset "Pow" begin
+        mesh = one_cell_mesh(:quad)
+        u = PhysicalFunction(x -> 3.0)
+        v = PhysicalFunction(x -> 2.0)
+        cinfo = CellInfo(mesh, 1)
+        cpoint = Bcube.CellPoint(SA[0.1, 0.3], cinfo, Bcube.ReferenceDomain())
+
+        @test Bcube.materialize(u * u, cinfo)(cpoint) ≈ 9
+        @test Bcube.materialize(u^2, cinfo)(cpoint) ≈ 9
+        @test Bcube.materialize(u^v, cinfo)(cpoint) ≈ 9
+
+        a = Bcube.NullOperator()
+        @test a^u == a
+    end
 end