GLT symbol of bilinear forms containing boundary integrals

gelato/expr.py

@@ -8,15 +8,15 @@
 from sympy import simplify, expand
 
 from sympde.expr import BilinearForm
-from sympde.expr import TensorExpr
+from sympde.expr import TensorExpr, IntAdd
 
 from sympde.expr import Mass as MassForm
 from sympde.expr import Stiffness as StiffnessForm
 from sympde.expr import Advection as AdvectionForm
 from sympde.expr import AdvectionT as AdvectionTForm
 from sympde.expr import Bilaplacian as BilaplacianForm
 from sympde.expr import Basic1dForm
-from sympde.topology import SymbolicExpr
+from sympde.topology import SymbolicExpr, Union
 from sympde.calculus.matrices import SymbolicDeterminant
 
 from .glt import (BasicGlt, Mass, Stiffness, Advection, Bilaplacian)
@@ -33,6 +33,22 @@ def gelatize(a, degrees=None, n_elements=None, evaluate=False, mapping=None,
     dim = a.ldim
     domain = a.domain
 
+    # ... set to zero if a only contains a boundary integral
+    if hasattr(domain, 'ext'):
+        return TensorExpr(0.)
+
+    # ... ignore the boundary integrals within the bilinear form
+    if isinstance(a.expr, IntAdd):
+        expr_list = []
+        for integ in a.expr.args:
+            if not hasattr(integ.domain, 'ext'):
+                expr_list.append(integ)
+
+        if len(expr_list) == 0:
+            return TensorExpr(0.)
+
+        a = BilinearForm(a.args[0], IntAdd(*expr_list))
+
     # ... compute tensor form
     expr = TensorExpr(a, domain=domain, mapping=mapping, expand=expand)
     # ...

gelato/tests/test_gelatize_1d.py

@@ -12,7 +12,7 @@
 from sympde.calculus import laplace, hessian, bracket, convect
 from sympde.topology import dx, dy, dz, dx1, dx2, dx3
 from sympde.topology import ScalarFunctionSpace
-from sympde.topology import Domain
+from sympde.topology import Domain, Line
 from sympde.topology import elements_of
 from sympde.expr import BilinearForm
 from sympde.expr import integral
@@ -102,6 +102,42 @@ def test_gelatize_1d_4():
     expr = BilinearForm((u,v), integral(domain, expr))
     assert(gelatize(expr) == expected)
 
+#==============================================================================
+def test_gelatize_1d_5():
+    domain = Line('Omega', bounds=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    expected = 0.0
+
+    expr = u*v
+    expr = BilinearForm((u,v), integral(domain.boundary, expr))
+
+    assert(gelatize(expr) == expected)
+
+#==============================================================================
+def test_gelatize_1d_6():
+    domain = Line('Omega', bounds=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    nx = symbols('nx', integer=True)
+    px = symbols('px', integer=True)
+    tx = symbols('tx')
+
+    expected = nx*Stiffness(px,tx)
+
+    expr = dot(grad(v), grad(u))
+
+    expr = BilinearForm((u,v), integral(domain.boundary, u*v) +
+                        integral(domain, expr))
+
+    assert(gelatize(expr) == expected)
+
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
 #==============================================================================

gelato/tests/test_gelatize_2d.py

@@ -12,7 +12,7 @@
 from sympde.calculus import laplace, hessian, bracket, convect
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import Domain
+from sympde.topology import Domain, Square
 from sympde.topology import Mapping
 from sympde.topology import elements_of
 from sympde.expr import BilinearForm
@@ -160,9 +160,45 @@ def test_gelatize_2d_7():
     expr = BilinearForm((u,v), integral(domain, expr))
     assert(gelatize(expr) == expected)
 
+#==============================================================================
+def test_gelatize_2d_8():
+    domain = Square('Omega', bounds1=(0,1), bounds2=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    expected = 0.0
+
+    expr = u*v
+    expr = BilinearForm((u,v), integral(domain.boundary, expr))
+
+    assert(gelatize(expr) == expected)
+
+#==============================================================================
+def test_gelatize_2d_9():
+    domain = Square('Omega', bounds1=(0,1), bounds2=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    nx, ny = symbols('nx ny', integer=True)
+    px, py = symbols('px py', integer=True)
+    tx, ty = symbols('tx ty')
+
+    expected = Mass(px,tx)*ny*Stiffness(py,ty)/nx + Mass(py,ty)*nx*Stiffness(px,tx)/ny
+
+    expr = dot(grad(v), grad(u))
+
+    expr = BilinearForm((u,v), integral(domain.boundary, u*v) +
+                        integral(domain, expr))
+
+    assert(gelatize(expr) == expected)
+
 #==============================================================================
 ## TODO
-#def test_gelatize_2d_8():
+#def test_gelatize_2d_10():
 #    domain = Domain('Omega', dim=DIM)
 #
 #    V = ScalarFunctionSpace('V', domain)

gelato/tests/test_gelatize_3d.py

@@ -12,7 +12,7 @@
 from sympde.calculus import laplace, hessian, bracket, convect
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
-from sympde.topology import Domain
+from sympde.topology import Domain, Cube
 from sympde.topology import Mapping
 from sympde.topology import elements_of
 from sympde.expr import BilinearForm
@@ -165,6 +165,43 @@ def test_gelatize_3d_7():
     expr = BilinearForm((u,v), integral(domain, expr))
     assert(gelatize(expr) == expected)
 
+#==============================================================================
+def test_gelatize_3d_8():
+    domain = Cube('Omega', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    expected = 0.0
+
+    expr = u*v
+    expr = BilinearForm((u,v), integral(domain.boundary, expr))
+
+    assert(gelatize(expr) == expected)
+
+#==============================================================================
+def test_gelatize_3d_9():
+    domain = Cube('Omega', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
+
+    V = ScalarFunctionSpace('V', domain)
+
+    u,v = elements_of(V, names='u,v')
+
+    nx, ny, nz = symbols('nx ny nz', integer=True)
+    px, py, pz = symbols('px py pz', integer=True)
+    tx, ty, tz = symbols('tx ty tz')
+
+    expected = ( nx*Mass(py,ty)*Mass(pz,tz)*Stiffness(px,tx)/(ny*nz) +
+                ny*Mass(px,tx)*Mass(pz,tz)*Stiffness(py,ty)/(nx*nz) +
+                nz*Mass(px,tx)*Mass(py,ty)*Stiffness(pz,tz)/(nx*ny))
+
+    expr = dot(grad(v), grad(u))
+
+    expr = BilinearForm((u,v), integral(domain.boundary, u*v) +
+                        integral(domain, expr))
+
+    assert(gelatize(expr) == expected)
 
 #==============================================================================
 ## TODO