From 0139092b5b601e1467a99adda2310634d1ed10ad Mon Sep 17 00:00:00 2001 From: saidctb Date: Sun, 15 May 2022 17:26:48 +0200 Subject: [PATCH 1/8] fix mp file bug --- sympde/expr/evaluation.py | 7 ++++- sympde/topology/domain.py | 55 +++++++++++++++++++++++++++++++++++--- sympde/topology/mapping.py | 2 ++ 3 files changed, 59 insertions(+), 5 deletions(-) diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index 58cf4f16..da06d9cb 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -594,7 +594,12 @@ def eval(cls, expr, domain): newcoords = [Symbol(c.name+"_plus") for c in coordinates] subs = list(zip(coordinates, newcoords)) J = J.subs(subs) - + elif mapping.is_plus: + M = list(J.atoms(type(mapping)))[0] + coordinates = [M[i] for i in range(domain.dim)] + newcoords = [mapping[i] for i in range(domain.dim)] + subs = list(zip(coordinates, newcoords)) + J = J.subs(subs) return J elif isinstance(expr, SymbolicDeterminant): diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index f1ace122..10048f89 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -104,10 +104,16 @@ def __new__(cls, name, *, interiors=None, boundaries=None, dim=None, if len(interiors) == 0 and dim is None: raise TypeError('No interior domain found') - else: + elif len(interiors) == 1: dtype = interiors[0].dtype dim = interiors[0].dim interiors = Union(*interiors) + else: + dim = interiors[0].dim + interiors = Union(*interiors) + dtype = [i.dtype for i in interiors] + + assert mapping is None and logical_domain is None or \ mapping is not None and logical_domain is not None @@ -297,13 +303,54 @@ def from_file( cls, filename ): d_interior = yml['interior'] d_boundary = yml['boundary'] d_connectivity = yml['connectivity'] - mapping = Mapping('{}_mapping'.format(domain_name), dim=int(dim)) if dtype == 'None': dtype = None if dtype is not None: - constructor = globals()[dtype['type']] - return mapping(constructor(domain_name, **dtype['parameters'])) + if len(d_interior)==1: + constructor = globals()[dtype['type']] + mapping = Mapping('mapping_0', dim=int(dim)) + domain = mapping(constructor(domain_name, **dtype['parameters'])) + else: + constructors = [globals()[dt['type']] for dt in dtype] + interiors = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)] + patch_index = {i['name']:ind for ind,i in enumerate(d_interior)} + mappings = [Mapping('mapping_{}'.format(i), dim=int(dim)) for i in range(len(interiors))] + domains = [mapping(i) for i,mapping in zip(interiors, mappings)] + + boundaries = [] + for bd in d_boundary: + name = bd['patch'] + axis = bd['axis'] + ext = bd['ext'] + i = patch_index[name] + bd = domains[i].get_boundary(axis=int(axis), ext=int(ext)) + boundaries.append(bd) + + interfaces = [] + for edge,(minus, plus) in d_connectivity.items(): + minus_name = minus['patch'] + minus_axis = minus['axis'] + minus_ext = minus['ext'] + minus_patch_i = patch_index[minus_name] + minus_bd = domains[minus_patch_i].get_boundary(axis=int(minus_axis), ext=int(minus_ext)) + + plus_name = plus['patch'] + plus_axis = plus['axis'] + plus_ext = plus['ext'] + plus_patch_i = patch_index[plus_name] + plus_bd = domains[plus_patch_i].get_boundary(axis=int(plus_axis), ext=int(plus_ext)) + + interfaces.append([minus_bd, plus_bd, 1]) + + domain = domains[0] + for p in domains[1:]: + domain = domain.join(p, name=domain_name) + + for I in interfaces: + domain = domain.join(domain, domain.name, bnd_minus=I[0], bnd_plus=I[1], direction=I[2]) + + return domain # ... create sympde InteriorDomain (s) interior = [InteriorDomain(i['name'], dim=dim) for i in d_interior] diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py index eadd4fba..83046381 100644 --- a/sympde/topology/mapping.py +++ b/sympde/topology/mapping.py @@ -1288,6 +1288,8 @@ def eval(cls, *_args, **kwargs): else: raise ValueError('Wrong index') + if base.is_plus: + name = name + '_plus' else: name = '{base}_{i}'.format(base=base.name, i=expr.indices[0]) From 4c5b429dacb81c457f930941b11e02983325a8ec Mon Sep 17 00:00:00 2001 From: saidctb Date: Mon, 16 May 2022 10:34:45 +0200 Subject: [PATCH 2/8] clean domain.from_file --- sympde/topology/basic.py | 12 +- sympde/topology/domain.py | 159 ++++++++++++------------- sympde/topology/tests/test_gallery.py | 25 ++-- sympde/topology/tests/test_topology.py | 56 ++++++--- 4 files changed, 135 insertions(+), 117 deletions(-) diff --git a/sympde/topology/basic.py b/sympde/topology/basic.py index 6ef66a0c..f8da7b67 100644 --- a/sympde/topology/basic.py +++ b/sympde/topology/basic.py @@ -103,7 +103,8 @@ def _sympystr(self, printer): return '{}'.format(sstr(self.name)) def todict(self): - return {'name': str(self.name)} + return {'name': str(self.logical_domain.name if self.logical_domain else self.name ), + 'mapping':str(self.mapping.name if self.mapping else None)} #============================================================================== @@ -351,11 +352,14 @@ def __add__(self, other): return Union(self, other) def todict(self): - return {'axis' : str(self.axis), + name = self.domain.logical_domain.name if self.logical_domain else self.domain.name + mapping = self.domain.mapping.name if self.domain.mapping else 'None' + d = {'axis' : str(self.axis), 'ext' : str(self.ext), 'name' : str(self.name), - 'patch': str(self.domain.name)} - + 'patch': str(name), + 'mapping':str(mapping)} + return d #============================================================================== class CornerBoundary(BasicDomain): """ diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index 10048f89..d1e65e6a 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -298,7 +298,7 @@ def from_file( cls, filename ): yml = yaml.load( h5['topology.yml'][()], Loader=yaml.SafeLoader ) domain_name = yml['name'] - dim = yml['dim'] + dim = int(yml['dim']) dtype = yml['dtype'] d_interior = yml['interior'] d_boundary = yml['boundary'] @@ -307,105 +307,110 @@ def from_file( cls, filename ): if dtype == 'None': dtype = None if dtype is not None: - if len(d_interior)==1: - constructor = globals()[dtype['type']] - mapping = Mapping('mapping_0', dim=int(dim)) - domain = mapping(constructor(domain_name, **dtype['parameters'])) - else: - constructors = [globals()[dt['type']] for dt in dtype] - interiors = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)] - patch_index = {i['name']:ind for ind,i in enumerate(d_interior)} - mappings = [Mapping('mapping_{}'.format(i), dim=int(dim)) for i in range(len(interiors))] - domains = [mapping(i) for i,mapping in zip(interiors, mappings)] - - boundaries = [] - for bd in d_boundary: - name = bd['patch'] - axis = bd['axis'] - ext = bd['ext'] - i = patch_index[name] - bd = domains[i].get_boundary(axis=int(axis), ext=int(ext)) - boundaries.append(bd) - - interfaces = [] - for edge,(minus, plus) in d_connectivity.items(): - minus_name = minus['patch'] - minus_axis = minus['axis'] - minus_ext = minus['ext'] - minus_patch_i = patch_index[minus_name] - minus_bd = domains[minus_patch_i].get_boundary(axis=int(minus_axis), ext=int(minus_ext)) - - plus_name = plus['patch'] - plus_axis = plus['axis'] - plus_ext = plus['ext'] - plus_patch_i = patch_index[plus_name] - plus_bd = domains[plus_patch_i].get_boundary(axis=int(plus_axis), ext=int(plus_ext)) - - interfaces.append([minus_bd, plus_bd, 1]) - - domain = domains[0] - for p in domains[1:]: - domain = domain.join(p, name=domain_name) - - for I in interfaces: - domain = domain.join(domain, domain.name, bnd_minus=I[0], bnd_plus=I[1], direction=I[2]) + if isinstance(d_interior, dict): + d_interior = [d_interior] + dtype = [dtype] + + if dtype is not None and all(dtype): + constructors = [globals()[dt['type']] for dt in dtype] + interiors = [cs(i['name'], **dt['parameters']) for cs,i,dt in zip(constructors, d_interior, dtype)] + mappings = [Mapping(I['mapping'], dim=dim) if I.get('mapping', "None") != "None" else None for I in d_interior] + domains = [mapping(i) if mapping else i for i,mapping in zip(interiors, mappings)] + patch_index = {I.name:ind for ind,I in enumerate(interiors)} + + boundaries = [] + for bd in d_boundary: + name = bd['patch'] + axis = bd['axis'] + ext = bd['ext'] + i = patch_index[name] + bd = domains[i].get_boundary(axis=int(axis), ext=int(ext)) + boundaries.append(bd) + + interfaces = [] + for edge,(minus, plus) in d_connectivity.items(): + minus_name = minus['patch'] + minus_axis = minus['axis'] + minus_ext = minus['ext'] + minus_patch_i = patch_index[minus_name] + minus_bd = domains[minus_patch_i].get_boundary(axis=int(minus_axis), ext=int(minus_ext)) + + plus_name = plus['patch'] + plus_axis = plus['axis'] + plus_ext = plus['ext'] + plus_patch_i = patch_index[plus_name] + plus_bd = domains[plus_patch_i].get_boundary(axis=int(plus_axis), ext=int(plus_ext)) + + interfaces.append([minus_bd, plus_bd, 1]) + + domain = domains[0] + for p in domains[1:]: + domain = domain.join(p, name=domain_name) + + for I in interfaces: + domain = domain.join(domain, domain.name, bnd_minus=I[0], bnd_plus=I[1], direction=I[2]) return domain # ... create sympde InteriorDomain (s) - interior = [InteriorDomain(i['name'], dim=dim) for i in d_interior] + l_interiors = [Domain(i['name'], dim=dim).interior for i in d_interior] + mappings = [Mapping(i['mapping'], dim=dim) if i['mapping'] != 'None' else None for i in d_interior] + l_interiors_index = {I.name:ind for ind,I in enumerate(l_interiors)} # create a dict of interiors accessed by name => needed for boundaries - d_interior = {} - for i in interior: - d_interior[i.name] = i - - if len(interior) == 1: - interior = interior[0] - # ... + d_l_interiors = {} + d_l_boundaries = {} + for i in l_interiors: + d_l_interiors[i.name] = i + d_l_boundaries[i.name] = [] # ... create sympde Boundary (s) - boundary = [] + for desc in d_boundary: name = desc['name'] - patch = d_interior[desc['patch']] + patch = d_l_interiors[desc['patch']] axis = desc['axis'] ext = desc['ext'] if axis == 'None': axis = None if ext == 'None': ext = None - bnd = Boundary(name, patch, axis=axis, ext=ext) - boundary.append(bnd) - - if len(boundary) == 1: - boundary = boundary[0] + lbnd = Boundary(name, patch, axis=axis, ext=ext) + d_l_boundaries[patch.name].append(lbnd) # ... - # ... create connectivity - connectivity = Connectivity() + subdomains = [] + for name,M in zip(d_l_interiors, mappings): + I = d_l_interiors[name] + B = d_l_boundaries[name] + D = Domain(I.name, dim=dim, interiors=I, boundaries=B) + D = M(D) if M else D + subdomains.append(D) + +# ... create connectivity + interfaces = [] for edge,pair in d_connectivity.items(): bnds = [] for desc in pair: - patch = d_interior[desc['patch']] + patch = d_l_interiors[desc['patch']] + patch_index = l_interiors_index[patch.name] name = desc['name'] bnd = Boundary(name, patch) + M = mappings[patch_index] + bnd = M(bnd) if M else bnd bnds.append(bnd) - connectivity[edge] = Interface(edge, bnds[0], bnds[1]) + interfaces.append(bnds) # ... - logical_domain = Domain(domain_name, dim=dim, - interiors=interior, - boundaries=boundary, - connectivity=connectivity) + domain = subdomains[0] + for p in subdomains[1:]: + domain = domain.join(p, name=domain_name) - obj = Domain.__new__(cls, domain_name, - interiors=interior, - boundaries=boundary, - connectivity=connectivity, - mapping = mapping, logical_domain=logical_domain) - return obj + for I in interfaces: + domain = domain.join(domain, domain.name, bnd_minus=I[0], bnd_plus=I[1]) + + return domain def join(self, other, name, bnd_minus=None, bnd_plus=None, direction=None): @@ -632,10 +637,7 @@ def __new__(cls, name, dim, min_coords, max_coords): raise ValueError("Min coordinates must be smaller than max") coord_names = 'x1:{}'.format(dim + 1) - coordinates = symbols(coord_names) - intervals = [Interval('{}_{}'.format(name, c.name), coordinate=c, bounds=(xmin, xmax)) - for c, xmin, xmax in zip(coordinates, min_coords, max_coords)] # Choose which type to use: # a) if dim <= 3, use Line, Square or Cube; @@ -673,12 +675,7 @@ def __new__(cls, name, dim, min_coords, max_coords): 'min_coords': [*min_coords], 'max_coords': [*max_coords]}} - if len(intervals) == 1: - interior = intervals[0] - else: - interior = ProductDomain(*intervals, name=name) - - interior = NCubeInterior(interior, dim=dim, dtype=dtype, min_coords=tuple(min_coords), max_coords=tuple(max_coords)) + interior = NCubeInterior(name, dim=dim, dtype=dtype, min_coords=tuple(min_coords), max_coords=tuple(max_coords)) # Create instance of given type obj = super().__new__(cls, name, interiors=[interior], boundaries=interior.boundary) diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py index dc7fa93a..9fd5c7d1 100644 --- a/sympde/topology/tests/test_gallery.py +++ b/sympde/topology/tests/test_gallery.py @@ -50,7 +50,7 @@ def test_unit_line(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== def test_generic_line(): @@ -81,7 +81,7 @@ def test_generic_line(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== def test_unit_square(): @@ -98,9 +98,7 @@ def test_unit_square(): # Domain's attributes assert isinstance(domain.interior, InteriorDomain) - assert isinstance(domain.interior.target, ProductDomain) - assert all(isinstance(i, Interval) for i in domain.interior.target.domains) - assert len(domain.interior.target.domains) == 2 + assert len(domain.boundary) == 4 assert domain.dtype == {'type': 'Square', 'parameters': {'bounds1': [0, 1], @@ -117,7 +115,7 @@ def test_unit_square(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== def test_rectangle(): @@ -134,9 +132,6 @@ def test_rectangle(): # Domain's attributes assert isinstance(domain.interior, InteriorDomain) - assert isinstance(domain.interior.target, ProductDomain) - assert all(isinstance(i, Interval) for i in domain.interior.target.domains) - assert len(domain.interior.target.domains) == 2 assert len(domain.boundary) == 4 assert domain.dtype == {'type': 'Square', 'parameters': {'bounds1': [1, 5], @@ -153,7 +148,7 @@ def test_rectangle(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== def test_unit_cube(): @@ -170,9 +165,6 @@ def test_unit_cube(): # Domain's attributes assert isinstance(domain.interior, InteriorDomain) - assert isinstance(domain.interior.target, ProductDomain) - assert all(isinstance(i, Interval) for i in domain.interior.target.domains) - assert len(domain.interior.target.domains) == 3 assert len(domain.boundary) == 6 assert domain.dtype == {'type': 'Cube', 'parameters': {'bounds1': [0, 1], @@ -191,7 +183,7 @@ def test_unit_cube(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== def test_orthogonal_hexahedron(): @@ -208,9 +200,6 @@ def test_orthogonal_hexahedron(): # Domain's attributes assert isinstance(domain.interior, InteriorDomain) - assert isinstance(domain.interior.target, ProductDomain) - assert all(isinstance(i, Interval) for i in domain.interior.target.domains) - assert len(domain.interior.target.domains) == 3 assert len(domain.boundary) == 6 assert domain.dtype == {'type': 'Cube', 'parameters': {'bounds1': [-1, 1], @@ -229,7 +218,7 @@ def test_orthogonal_hexahedron(): # Export to file, read it again and compare with original domain domain.export('domain.h5') D = Domain.from_file('domain.h5') - assert D.logical_domain == domain + assert D == domain #============================================================================== # CLEAN UP SYMPY NAMESPACE diff --git a/sympde/topology/tests/test_topology.py b/sympde/topology/tests/test_topology.py index 11da9f07..e0e8053e 100644 --- a/sympde/topology/tests/test_topology.py +++ b/sympde/topology/tests/test_topology.py @@ -21,8 +21,8 @@ def test_interior_domain(): D1 = InteriorDomain('D1', dim=2) D2 = InteriorDomain('D2', dim=2) - assert( D1.todict() == {'name': 'D1'} ) - assert( D2.todict() == {'name': 'D2'} ) + assert( D1.todict() == {'name': 'D1', 'mapping':'None'} ) + assert( D2.todict() == {'name': 'D2', 'mapping':'None'} ) assert( Union(D2, D1) == Union(D1, D2) ) @@ -30,8 +30,8 @@ def test_interior_domain(): assert(D.dim == 2) assert(len(D) == 2) - assert( D.todict() == [{'name': 'D1'}, - {'name': 'D2'}] ) + assert( D.todict() == [{'name': 'D1', 'mapping':'None'}, + {'name': 'D2', 'mapping':'None'}] ) #============================================================================== def test_topology_1(): @@ -49,20 +49,47 @@ def test_topology_1(): bnd_B_2 = Boundary('Gamma_2', B) bnd_B_3 = Boundary('Gamma_3', B) - connectivity['I'] = Interface('I', bnd_A_1, bnd_B_2) - - mapping = Mapping('M', dim=2) - logical_domain = Domain('D', dim=2) + mapping = Mapping('M', dim=2) interiors = [A, B] boundaries = [bnd_A_2, bnd_A_3, bnd_B_1, bnd_B_3] - Omega = Domain('Omega', + lOmega = Domain('Omega', interiors=interiors, boundaries=boundaries, - connectivity=connectivity, - mapping=mapping, - logical_domain=logical_domain) + connectivity=connectivity) + + Omega = mapping(lOmega) + Omega = Omega.join(Omega, name=Omega.name, bnd_minus=mapping(bnd_A_1), bnd_plus=mapping(bnd_B_2)) + + interfaces = Omega.interfaces + assert(isinstance(interfaces, Interface)) + + # export + Omega.export('omega.h5') + # ... + + # read it again and check that it has the same description as Omega + D = Domain.from_file('omega.h5') + assert( D.todict() == Omega.todict() ) + +#============================================================================== +def test_topology_2(): + + # ... create a domain with 2 subdomains D1 and D2 + A = Square('A') + B = Square('B') + # ... + + M1 = Mapping('M1', dim=2) + M2 = Mapping('M2', dim=2) + + D1 = M1(A) + D2 = M2(B) + + Omega = D1.join(D2, name = 'Omega', + bnd_minus = D1.get_boundary(axis=0, ext=1), + bnd_plus = D2.get_boundary(axis=0, ext=-1)) interfaces = Omega.interfaces assert(isinstance(interfaces, Interface)) @@ -73,6 +100,7 @@ def test_topology_1(): # read it again and check that it has the same description as Omega D = Domain.from_file('omega.h5') + assert( D.todict() == Omega.todict() ) #============================================================================== @@ -173,7 +201,7 @@ def test_domain_join_line(): print(AB) assert AB.interior == Union(A.interior, B.interior) - assert AB.interfaces == Interface('A_x1|B_x1', AB_bnd_minus, AB_bnd_plus) + assert AB.interfaces == Interface('A|B', AB_bnd_minus, AB_bnd_plus) print(AB.connectivity) print('') # ... @@ -189,7 +217,7 @@ def test_domain_join_line(): print(ABC) assert ABC.interior == Union(A.interior, B.interior, C.interior) - assert ABC.interfaces == Union(Interface('A_x1|B_x1', AB_bnd_minus, AB_bnd_plus),Interface('B_x1|C_x1', BC_bnd_minus, BC_bnd_plus)) + assert ABC.interfaces == Union(Interface('A|B', AB_bnd_minus, AB_bnd_plus),Interface('B|C', BC_bnd_minus, BC_bnd_plus)) print(list(ABC.connectivity.items())) print('') # ... From d741c7f9a53da8ac4b154aeeec834c42e0dbc343 Mon Sep 17 00:00:00 2001 From: saidctb Date: Mon, 16 May 2022 11:17:05 +0200 Subject: [PATCH 3/8] remove unused variables --- sympde/topology/domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/topology/domain.py b/sympde/topology/domain.py index d1e65e6a..0c237b2a 100644 --- a/sympde/topology/domain.py +++ b/sympde/topology/domain.py @@ -328,7 +328,7 @@ def from_file( cls, filename ): boundaries.append(bd) interfaces = [] - for edge,(minus, plus) in d_connectivity.items(): + for _,(minus, plus) in d_connectivity.items(): minus_name = minus['patch'] minus_axis = minus['axis'] minus_ext = minus['ext'] From 5189605ad907a081e94201766c1a272de6ad9944 Mon Sep 17 00:00:00 2001 From: saidctb Date: Wed, 1 Jun 2022 00:23:11 +0200 Subject: [PATCH 4/8] generalize the div operator to take a matrix as an argument --- sympde/calculus/core.py | 1 - sympde/expr/evaluation.py | 4 ++-- sympde/expr/expr.py | 2 +- sympde/topology/derivatives.py | 23 +++++++++++++++++++++++ sympde/topology/space.py | 10 ++++++++-- 5 files changed, 34 insertions(+), 6 deletions(-) diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py index 82919187..ac831aaa 100644 --- a/sympde/calculus/core.py +++ b/sympde/calculus/core.py @@ -366,7 +366,6 @@ def __new__(cls, arg1, arg2, **options): return c*obj - #============================================================================== #TODO add cross(u,v) + cross(v,u) = 0 class Cross(BasicOperator): diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py index da06d9cb..7cd02bbe 100644 --- a/sympde/expr/evaluation.py +++ b/sympde/expr/evaluation.py @@ -344,7 +344,7 @@ def _nullify(expr, u, us): newexpr = newexpr.subs(mapping, mapping.plus) for nn in newexpr.atoms(NormalVector): - newexpr = newexpr.subs(nn, -nn) + newexpr = newexpr.replace(nn, -nn) if not is_zero(newexpr): if interface.plus in bnd_expressions: @@ -760,7 +760,7 @@ def eval(cls, expr, domain): elif isinstance(expr, NormalVector): lines = [[expr[i] for i in range(dim)]] - return ImmutableDenseMatrix(lines) + return ImmutableDenseMatrix(lines).T elif isinstance(expr, TangentVector): lines = [[expr[i] for i in range(dim)]] diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 7428df7f..7e4d00da 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -670,7 +670,7 @@ def is_linear_expression(expr, args, integral=True, debug=True): # ... left_args = [] right_args = [] - + return True for arg in args: tag = random_string( 4 ) diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index 094e9432..b1b6ed0a 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -732,6 +732,11 @@ def eval(cls, *_args): return u = _args[0] + if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: + div = [] + for j in range(u.shape[1]): + div.append(dx(u[0,j])+dy(u[1,j])) + return ImmutableDenseMatrix(div) return dx(u[0]) + dy(u[1]) @@ -745,6 +750,12 @@ def eval(cls, *_args): u = _args[0] + if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: + div = [] + for j in range(u.shape[1]): + div.append(dx(u[0,j])+dy(u[1,j])+dz(u[2,j])) + return ImmutableDenseMatrix(div) + return dx(u[0]) + dy(u[1]) + dz(u[2]) #============================================================================== @@ -1056,6 +1067,12 @@ def eval(cls, *_args): u = _args[0] + if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: + div = [] + for j in range(u.shape[1]): + div.append(dx1(u[0,j])+dx2(u[1,j])) + return ImmutableDenseMatrix(div) + return dx1(u[0]) + dx2(u[1]) class LogicalDiv_3d(DivBasic): @@ -1068,6 +1085,12 @@ def eval(cls, *_args): u = _args[0] + if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: + div = [] + for j in range(u.shape[1]): + div.append(dx1(u[0,j])+dx2(u[1,j])+dx3(u[2,j])) + return ImmutableDenseMatrix(div) + return dx1(u[0]) + dx2(u[1]) + dx3(u[2]) #============================================================================== diff --git a/sympde/topology/space.py b/sympde/topology/space.py index a3b40e80..cce1f556 100644 --- a/sympde/topology/space.py +++ b/sympde/topology/space.py @@ -363,6 +363,7 @@ def __new__(cls, *spaces): name = ''.join(i.name for i in spaces) # ... + assert all(i.is_broken for i in spaces) or all(not i.is_broken for i in spaces) # ... def _get_name(V): if V.ldim == 1: @@ -383,9 +384,10 @@ def _get_name(V): # ... obj = Basic.__new__(cls, spaces) - obj._shape = shape + obj._shape = shape obj._coordinates = coordinates - obj._name = name + obj._name = name + obj._is_broken = spaces[0].is_broken # ... return obj @@ -523,6 +525,10 @@ def free_symbols(self): def ldim(self): return self.base.space.ldim + @property + def space(self): + return self.base.space + def __hash__(self): return hash(self._args) From bf0c473991479158484e2915db19f4f59a9080c7 Mon Sep 17 00:00:00 2001 From: Said Hadjout Date: Thu, 30 Jun 2022 21:41:32 +0200 Subject: [PATCH 5/8] Update expr.py --- sympde/expr/expr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py index 7e4d00da..7428df7f 100644 --- a/sympde/expr/expr.py +++ b/sympde/expr/expr.py @@ -670,7 +670,7 @@ def is_linear_expression(expr, args, integral=True, debug=True): # ... left_args = [] right_args = [] - return True + for arg in args: tag = random_string( 4 ) From 0f15fde0582a0488d3b6d414565ee82cbf04c749 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 21 Jul 2022 17:37:07 +0200 Subject: [PATCH 6/8] Make Grad code simpler --- sympde/topology/derivatives.py | 75 ++++++++++------------------------ 1 file changed, 21 insertions(+), 54 deletions(-) diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index b1b6ed0a..c79f1812 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -520,22 +520,31 @@ def eval(cls, *_args): #============================================================================== class GradBasic(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Grad' @cacheit - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) + return cls.eval(u) else: - r = None + return Basic.__new__(cls, u, **options) - if r is None: - return Basic.__new__(cls, *args, **options) + @classmethod + @cacheit + def eval(cls, u): + du = [dxi(u) for dxi in cls._derivatives] + + # 1D case: scalar output + if len(du) == 1: return du[0] + + # 2D and 3D cases + if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)): + lines = [list(d[:]) for d in du] else: - return r + lines = [[d] for d in du] + return ImmutableDenseMatrix(lines) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -544,59 +553,17 @@ def __getitem__(self, indices, **kw_args): else: return Indexed(self, indices, **kw_args) -class Grad_1d(GradBasic): - - @classmethod - @cacheit - def eval(cls, *_args): - - if not _args: - return - u = _args[0] +class Grad_1d(GradBasic): + _derivatives = (dx,) - return dx(u) class Grad_2d(GradBasic): + _derivatives = (dx, dy) - @classmethod - @cacheit - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - du = (dx(u), dy(u)) - - if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)): - lines = [list(d[:]) for d in du] - else: - lines = [[d] for d in du] - - v = ImmutableDenseMatrix(lines) - return v class Grad_3d(GradBasic): - - @classmethod - @cacheit - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - du = (dx(u), dy(u), dz(u)) - - if isinstance(du[0], (Tuple, Matrix, ImmutableDenseMatrix)): - lines = [list(d[:]) for d in du] - else: - lines = [[d] for d in du] - - v = ImmutableDenseMatrix(lines) - - return v + _derivatives = (dx, dy, dz) #============================================================================== class CurlBasic(CalculusFunction): From 7b234fa60ca6038c6cb5c1ec6aab20397ee2e40e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 21 Jul 2022 17:42:39 +0200 Subject: [PATCH 7/8] Simplify code for Curl and Rot --- sympde/topology/derivatives.py | 54 ++++++++-------------------------- 1 file changed, 13 insertions(+), 41 deletions(-) diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index c79f1812..99540100 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -568,21 +568,15 @@ class Grad_3d(GradBasic): #============================================================================== class CurlBasic(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Curl' - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) - else: - r = None - - if r is None: - return Basic.__new__(cls, *args, **options) + return cls.eval(u) else: - return r + return Basic.__new__(cls, u, **options) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -591,28 +585,18 @@ def __getitem__(self, indices, **kw_args): else: return Indexed(self, indices, **kw_args) + class Curl_2d(CurlBasic): @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] + def eval(cls, u): + return dx(u[1]) - dy(u[0]) - return dx(u[1])-dy(u[0]) class Curl_3d(CurlBasic): @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - + def eval(cls, u): return ImmutableDenseMatrix([[dy(u[2]) - dz(u[1])], [dz(u[0]) - dx(u[2])], [dx(u[1]) - dy(u[0])]]) @@ -623,18 +607,12 @@ class Rot_2d(CalculusFunction): nargs = None name = 'Grad' - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) - else: - r = None - - if r is None: - return Basic.__new__(cls, *args, **options) + return cls.eval(u) else: - return r + return Basic.__new__(cls, u, **options) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -644,14 +622,8 @@ def __getitem__(self, indices, **kw_args): return Indexed(self, indices, **kw_args) @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - - return ImmutableDenseMatrix([[dy(u)],[-dx(u)]]) + def eval(cls, u): + return ImmutableDenseMatrix([[dy(u)], [-dx(u)]]) #============================================================================== class DivBasic(CalculusFunction): From 69f5b91d41059f4c0dd87abf6fd3522ce9557632 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 21 Jul 2022 19:23:29 +0200 Subject: [PATCH 8/8] Simplify code for Div, Laplace and Hessian --- sympde/topology/derivatives.py | 210 ++++++++++++++------------------- 1 file changed, 88 insertions(+), 122 deletions(-) diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py index 99540100..57c67fe2 100644 --- a/sympde/topology/derivatives.py +++ b/sympde/topology/derivatives.py @@ -604,7 +604,7 @@ def eval(cls, u): #============================================================================== class Rot_2d(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Grad' def __new__(cls, u, **options): @@ -628,21 +628,15 @@ def eval(cls, u): #============================================================================== class DivBasic(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Div' - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) - else: - r = None - - if r is None: - return Basic.__new__(cls, *args, **options) + return cls.eval(u) else: - return r + return Basic.__new__(cls, u, **options) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -651,70 +645,78 @@ def __getitem__(self, indices, **kw_args): else: return Indexed(self, indices, **kw_args) -class Div_1d(DivBasic): - - @classmethod - def eval(cls, *_args): +# @classmethod +# def scalar_eval(cls, u): +# return sum(dxi(u[i]) for i, dxi in enumerate(cls._derivatives)) +# +# @classmethod +# def eval(cls, u): +# if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1: +# div = [cls.scalar_eval(u[:, j]) for j in range(u.shape[1])] +# return ImmutableDenseMatrix(div) +# else: +# return cls.scalar_eval(u) - if not _args: - return - - u = _args[0] - return dx(u[0]) - -class Div_2d(DivBasic): @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: - div = [] - for j in range(u.shape[1]): - div.append(dx(u[0,j])+dy(u[1,j])) + def eval(cls, u): + if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1: + div = [cls.eval(u[:, j]) for j in range(u.shape[1])] return ImmutableDenseMatrix(div) + else: + return sum(dxi(u[i]) for i, dxi in enumerate(cls._derivatives)) - return dx(u[0]) + dy(u[1]) -class Div_3d(DivBasic): +class Div_1d(DivBasic): + _derivatives = (dx, ) - @classmethod - def eval(cls, *_args): - if not _args: - return +class Div_2d(DivBasic): + _derivatives = (dx, dy) - u = _args[0] - if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1]>1: - div = [] - for j in range(u.shape[1]): - div.append(dx(u[0,j])+dy(u[1,j])+dz(u[2,j])) - return ImmutableDenseMatrix(div) +class Div_3d(DivBasic): + _derivatives = (dx, dy, dz) - return dx(u[0]) + dy(u[1]) + dz(u[2]) + +#class Div_1d(DivBasic): +# @classmethod +# def eval(cls, u): +# return dx(u[0]) +# +# +#class Div_2d(DivBasic): +# @classmethod +# def eval(cls, u): +# if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1: +# div = [dx(u[0, j]) + dy(u[1, j]) for j in range(u.shape[1])] +# return ImmutableDenseMatrix(div) +# else: +# return dx(u[0]) + dy(u[1]) +# +# +#class Div_3d(DivBasic): +# @classmethod +# def eval(cls, u): +# if isinstance(u, (Matrix, ImmutableDenseMatrix)) and u.shape[1] > 1: +# div = [dx(u[0, j]) + dy(u[1, j]) + dz(u[2, j]) +# for j in range(u.shape[1])] +# return ImmutableDenseMatrix(div) +# else: +# return dx(u[0]) + dy(u[1]) + dz(u[2]) #============================================================================== class LaplaceBasic(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Laplace' - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) - else: - r = None - - if r is None: - return Basic.__new__(cls, *args, **options) + return cls.eval(u) else: - return r + return Basic.__new__(cls, u, **options) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -723,66 +725,37 @@ def __getitem__(self, indices, **kw_args): else: return Indexed(self, indices, **kw_args) -class Laplace_1d(LaplaceBasic): - @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] + def eval(cls, u): if isinstance(u, VectorFunction): raise NotImplementedError('TODO') - return dx(dx(u)) + return sum(dxi(dxi(u)) for dxi in cls._derivatives) -class Laplace_2d(LaplaceBasic): - @classmethod - def eval(cls, *_args): +class Laplace_1d(LaplaceBasic): + _derivatives = (dx,) - if not _args: - return - u = _args[0] - if isinstance(u, VectorFunction): - raise NotImplementedError('TODO') +class Laplace_2d(LaplaceBasic): + _derivatives = (dx, dy) - return dx(dx(u)) + dy(dy(u)) class Laplace_3d(LaplaceBasic): - - @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] - if isinstance(u, VectorFunction): - raise NotImplementedError('TODO') - - return dx(dx(u)) + dy(dy(u)) + dz(dz(u)) + _derivatives = (dx, dy, dz) #============================================================================== class HessianBasic(CalculusFunction): - nargs = None + nargs = None # TODO [YG, 21.07.2022]: Check why we need this name = 'Hessian' - def __new__(cls, *args, **options): + def __new__(cls, u, **options): # (Try to) sympify args first - if options.pop('evaluate', True): - r = cls.eval(*args) - else: - r = None - - if r is None: - return Basic.__new__(cls, *args, **options) + return cls.eval(u) else: - return r + return Basic.__new__(cls, u, **options) def __getitem__(self, indices, **kw_args): if is_sequence(indices): @@ -791,50 +764,43 @@ def __getitem__(self, indices, **kw_args): else: return Indexed(self, indices, **kw_args) -class Hessian_1d(HessianBasic): +# @classmethod +# def eval(cls, u): +# if isinstance(u, VectorFunction): +# raise NotImplementedError('TODO') +# +# hessian = ImmutableDenseMatrix( +# [[di(dj(u)) for dj in cls._derivatives] +# for di in cls._derivatives] +# ) +# +# return hessian[0, 0] if hessian.shape == (1, 1) else hessian +class Hessian_1d(HessianBasic): @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] + def eval(cls, u): if isinstance(u, VectorFunction): raise NotImplementedError('TODO') + return dx(u) - return dx(dx(u)) class Hessian_2d(HessianBasic): - @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] + def eval(cls, u): if isinstance(u, VectorFunction): raise NotImplementedError('TODO') - return ImmutableDenseMatrix([[dx(dx(u)), dx(dy(u))], - [dx(dy(u)), dy(dy(u))]]) + [dx(dy(u)), dy(dy(u))]]) -class Hessian_3d(HessianBasic): +class Hessian_3d(HessianBasic): @classmethod - def eval(cls, *_args): - - if not _args: - return - - u = _args[0] + def eval(cls, u): if isinstance(u, VectorFunction): raise NotImplementedError('TODO') - return ImmutableDenseMatrix([[dx(dx(u)), dx(dy(u)), dx(dz(u))], - [dx(dy(u)), dy(dy(u)), dy(dz(u))], - [dx(dz(u)), dy(dz(u)), dz(dz(u))]]) + [dx(dy(u)), dy(dy(u)), dy(dz(u))], + [dx(dz(u)), dy(dz(u)), dz(dz(u))]]) #============================================================================== class BracketBasic(CalculusFunction):