Skip to content

[not meant to be merged] exploratory changes to improve single and multi-patch compatibility #467

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
campospinto wants to merge 2 commits into
base: devel
Choose a base branch
from
Draft

[not meant to be merged] exploratory changes to improve single and multi-patch compatibility #467

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8f71751
exploratory changes to make multipatch solver working with single pat…
campospinto Feb 21, 2025
4ed86dd
Merge branch 'devel' into multipatch_compatible
campospinto Mar 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion psydac/api/feec.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
from psydac.fem.basic import FemSpace
from psydac.fem.vector import VectorFemSpace


__all__ = ('DiscreteDerham',)

#==============================================================================
Expand Down Expand Up @@ -263,4 +264,3 @@ def projectors(self, *, kind='global', nquads=None):
return P0, P1, P2, P3, Pvec
else :
return P0, P1, P2, P3

146 changes: 98 additions & 48 deletions psydac/feec/multipatch/api.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,20 @@

from sympde.topology import Derham
from sympde.topology import element_of, elements_of
from sympde.topology.mapping import MultiPatchMapping
from sympde.topology.space import ScalarFunction
from sympde.calculus import grad, dot, inner, rot, div
from sympde.calculus import laplace, bracket, convect
from sympde.calculus import jump, avg, Dn, minus, plus
from sympde.expr.expr import LinearForm, BilinearForm, integral


from psydac.api.settings import PSYDAC_BACKENDS

from psydac.api.discretization import discretize as discretize_single_patch
from psydac.api.discretization import discretize_space
from psydac.api.discretization import DiscreteDerham
from psydac.feec.pull_push import pull_2d_h1, pull_2d_hcurl, pull_2d_l2
from psydac.feec.multipatch.operators import BrokenGradient_2D
from psydac.feec.multipatch.operators import BrokenScalarCurl_2D
from psydac.feec.multipatch.operators import Multipatch_Projector_H1
Expand Down Expand Up @@ -89,6 +92,13 @@ def __init__(self, *, mapping, domain_h, spaces, sequence=None):
else:
raise ValueError('Dimension {} is not available'.format(dim))

if mapping:
assert isinstance(mapping, MultiPatchMapping)
mappings = mapping.args[0]
self._mappings_list = list(mappings.values())
else:
self._mappings_list = list([])

#--------------------------------------------------------------------------
@property
def sequence(self):
Expand All @@ -104,6 +114,20 @@ def broken_derivatives_as_operators(self):
def broken_derivatives_as_matrices(self):
return tuple(b_diff.matrix for b_diff in self._broken_diff_ops)


# for compatibility with DiscreteDerham sequence:
@property
def derivatives(self):
return self.broken_derivatives_as_operators

@property
def derivatives_as_matrices(self):
return self.broken_derivatives_as_matrices

@property
def mappings_list(self):
return self._mappings_list

#--------------------------------------------------------------------------
def projectors(self, *, kind='global', nquads=None):
"""
Expand Down Expand Up @@ -150,6 +174,25 @@ def projectors(self, *, kind='global', nquads=None):
raise NotImplementedError('2D sequence with H-div not available yet')

P2 = Multipatch_Projector_L2(self.V2, nquads=nquads)

if self.mappings_list:
P0_m = lambda f: P0([pull_2d_h1(f, m.get_callable_mapping())
for m in self.mappings_list])

P2_m = lambda f: P2([pull_2d_l2(f, m.get_callable_mapping())
for m in self.mappings_list])
if self.sequence[1] == 'hcurl':
# f_x = lambdify(domain.coordinates, f_phys[0])
# f_y = lambdify(domain.coordinates, f_phys[1])
# f_log = [pull_2d_hcurl([f_x, f_y], m.get_callable_mapping())
# for m in mappings_list]
P1_m = lambda f: P1([pull_2d_hcurl(f, m.get_callable_mapping())
for m in self.mappings_list])
else:
raise NotImplementedError('2D sequence with H-div not available yet')

return P0_m, P1_m, P2_m

return P0, P1, P2

elif self.dim == 3:
Expand Down Expand Up @@ -225,61 +268,68 @@ def conforming_projection(self, space, hom_bc=False, backend_language="python",

return cP

def get_dual_dofs(self, space, f, backend_language="python", return_format='stencil_array'):
"""
return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(V^k)) of a given function f, as a stencil array or numpy array

Parameters
----------
space : <str>
The space of the dual dofs

f : <sympy.Expr>
The function used for evaluation

backend_language: <str>
The backend used to accelerate the code

return_format: <str>
The format of the dofs, can be 'stencil_array' or 'numpy_array'

Returns
-------
tilde_f:<Vector|ndarray>
The dual dofs
"""
if space == 'V0':
Vh = self.V0
elif space == 'V1':
Vh = self.V1
elif space == 'V2':
Vh = self.V2
else:
raise NotImplementedError("The space of kind {} is not available".format(space))

V = Vh.symbolic_space
v = element_of(V, name='v')

if isinstance(v, ScalarFunction):
expr = f*v
else:
expr = dot(f,v)

l = LinearForm(v, integral( V.domain, expr))
lh = discretize(l, self._domain_h, Vh, backend=PSYDAC_BACKENDS[backend_language])
tilde_f = lh.assemble()

if return_format == 'numpy_array':
return tilde_f.toarray()
else:
return tilde_f
# def get_dual_dofs(self, space, f, backend_language="python", return_format='stencil_array'):
# """
# return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(V^k)) of a given function f, as a stencil array or numpy array

# Parameters
# ----------
# space : <str>
# The space of the dual dofs

# f : <sympy.Expr>
# The function used for evaluation

# backend_language: <str>
# The backend used to accelerate the code

# return_format: <str>
# The format of the dofs, can be 'stencil_array' or 'numpy_array'

# Returns
# -------
# tilde_f:<Vector|ndarray>
# The dual dofs
# """
# if space == 'V0':
# Vh = self.V0
# elif space == 'V1':
# Vh = self.V1
# elif space == 'V2':
# Vh = self.V2
# else:
# raise NotImplementedError("The space of kind {} is not available".format(space))

# V = Vh.symbolic_space
# v = element_of(V, name='v')

# if isinstance(v, ScalarFunction):
# expr = f*v
# else:
# expr = dot(f,v)

# l = LinearForm(v, integral( V.domain, expr))
# lh = discretize(l, self._domain_h, Vh, backend=PSYDAC_BACKENDS[backend_language])
# tilde_f = lh.assemble()

# if return_format == 'numpy_array':
# return tilde_f.toarray()
# else:
# return tilde_f

#==============================================================================
def discretize_derham_multipatch(derham, domain_h, *args, **kwargs):

ldim = derham.shape
mapping = derham.spaces[0].domain.mapping

print(f'[o] mapping: ', mapping)
print(f'[o] type(mapping): ', type(mapping))
print(f'[o] mapping.args: ', mapping.args)
# print(f'[o] type(mapping[0]): ', type(mapping[0]))
print(f'[o] type(mapping.args[0]): ', type(mapping.args[0]))
# print(f'[o] type(mapping): ', type(mapping))

bases = ['B'] + ldim * ['M']
spaces = [discretize_space(V, domain_h, *args, basis=basis, **kwargs) \
for V, basis in zip(derham.spaces, bases)]
Expand Down
9 changes: 8 additions & 1 deletion psydac/feec/multipatch/examples/ppc_test_cases.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
import os
import numpy as np

from sympy import lambdify
from sympy import pi, cos, sin, Tuple, exp, atan, atan2

from sympde.topology import Derham
Expand Down Expand Up @@ -336,7 +337,13 @@ def get_source_and_solution_hcurl(
else:
raise ValueError(source_type)

# u_ex = Tuple(0, 1) # DEBUG
u_bc = [lambdify(domain.coordinates, u_bc[i]) for i in [0,1]]
u_ex = [lambdify(domain.coordinates, u_ex[i]) for i in [0,1]]

# f_y = lambdify(domain.coordinates, f_phys[1])

# u_bc = lambdify(domain.coordinates, u_bc)
# u_ex = lambdify(domain.coordinates, u_ex)
return f_vect, u_bc, u_ex, curl_u_ex, div_u_ex # , phi, grad_phi


Expand Down
2 changes: 1 addition & 1 deletion psydac/feec/multipatch/fem_linear_operators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ def to_sparse_matrix( self , **kwargs):
elif self._matrix is not None:
return self._matrix.tosparse()
else:
raise NotImplementedError('Class does not provide a get_sparse_matrix() method without a matrix')
raise NotImplementedError('Class does not provide a to_sparse_matrix() method without a matrix')

# ...
def __call__( self, f ):
Expand Down
11 changes: 9 additions & 2 deletions psydac/feec/multipatch/multipatch_domain_utilities.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -954,8 +954,12 @@ def build_cartesian_multipatch_domain(ncells, log_interval_x, log_interval_y, ma
"""
ax, bx = log_interval_x
ay, by = log_interval_y

print(f'type(ncells) = {type(ncells)}')
ncells = np.array(ncells)
nb_patchx, nb_patchy = np.shape(ncells)

print(f'type(ncells) = {type(ncells)}')

# equidistant logical patches
# ensure the following lists have the same shape as ncells
list_log_patches = [[Square('Log_' + str(j) + '_' + str(i),
Expand All @@ -976,10 +980,13 @@ def build_cartesian_multipatch_domain(ncells, log_interval_x, log_interval_y, ma
for i in range(nb_patchy)] for j in range(nb_patchx)]

# flatten for the join function
# print(f'ncells = {ncells}, np.shape(ncells) = {np.shape(ncells)}')
patches = []
for i in range(nb_patchx):
for j in range(nb_patchy):
if ncells[i, j] is not None:
# print(f'nb_patchx, nb_patchy, i, j = {nb_patchx, nb_patchy, i, j}')
if ncells[i, j] is not None: # raises error (TypeError: list indices must be integers or slices, not tuple) but why now?
# if ncells[i][j] is not None:
patches.append(list_patches[j][i])

axis_0 = 0
Expand Down
Loading