From dcb5b55449e8ba39252c5a719cfc729069467954 Mon Sep 17 00:00:00 2001 From: e-moral-sanchez Date: Tue, 27 Jun 2023 16:11:59 +0200 Subject: [PATCH] eliminate boundary integrals before computing symbol --- gelato/expr.py | 20 ++++++++++++++-- gelato/tests/test_gelatize_1d.py | 38 +++++++++++++++++++++++++++++- gelato/tests/test_gelatize_2d.py | 40 ++++++++++++++++++++++++++++++-- gelato/tests/test_gelatize_3d.py | 39 ++++++++++++++++++++++++++++++- 4 files changed, 131 insertions(+), 6 deletions(-) diff --git a/gelato/expr.py b/gelato/expr.py index 299ab0f..e3a73d4 100644 --- a/gelato/expr.py +++ b/gelato/expr.py @@ -8,7 +8,7 @@ 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 @@ -16,7 +16,7 @@ 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) # ... diff --git a/gelato/tests/test_gelatize_1d.py b/gelato/tests/test_gelatize_1d.py index 59f68fa..08680e8 100644 --- a/gelato/tests/test_gelatize_1d.py +++ b/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 #============================================================================== diff --git a/gelato/tests/test_gelatize_2d.py b/gelato/tests/test_gelatize_2d.py index 13a4b32..d559cde 100644 --- a/gelato/tests/test_gelatize_2d.py +++ b/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) diff --git a/gelato/tests/test_gelatize_3d.py b/gelato/tests/test_gelatize_3d.py index f180670..1c9c79a 100644 --- a/gelato/tests/test_gelatize_3d.py +++ b/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