From 8f7175120713f1fbf5d8d1a96735697e3d14181d Mon Sep 17 00:00:00 2001 From: Martin Campos Pinto Date: Fri, 21 Feb 2025 10:18:44 +0100 Subject: [PATCH] exploratory changes to make multipatch solver working with single patch domain --- psydac/api/feec.py | 2 +- psydac/feec/multipatch/api.py | 146 ++++++++++++------ .../multipatch/examples/ppc_test_cases.py | 9 +- .../feec/multipatch/fem_linear_operators.py | 2 +- .../multipatch/multipatch_domain_utilities.py | 11 +- .../feec/multipatch/non_matching_operators.py | 128 ++++++++------- psydac/feec/multipatch/operators.py | 38 +++-- psydac/feec/multipatch/plotting_utilities.py | 66 +++++--- psydac/fem/basic.py | 36 ++++- psydac/fem/projectors.py | 75 ++++++++- psydac/fem/splines.py | 12 ++ psydac/fem/tensor.py | 12 ++ psydac/fem/vector.py | 23 +++ 13 files changed, 410 insertions(+), 150 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index cfec097eb..2c154497a 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -12,6 +12,7 @@ from psydac.fem.basic import FemSpace from psydac.fem.vector import VectorFemSpace + __all__ = ('DiscreteDerham',) #============================================================================== @@ -263,4 +264,3 @@ def projectors(self, *, kind='global', nquads=None): return P0, P1, P2, P3, Pvec else : return P0, P1, P2, P3 - diff --git a/psydac/feec/multipatch/api.py b/psydac/feec/multipatch/api.py index 566815421..0bb45dd40 100644 --- a/psydac/feec/multipatch/api.py +++ b/psydac/feec/multipatch/api.py @@ -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 @@ -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): @@ -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): """ @@ -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: @@ -225,54 +268,54 @@ 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 : - The space of the dual dofs - - f : - The function used for evaluation - - backend_language: - The backend used to accelerate the code - - return_format: - The format of the dofs, can be 'stencil_array' or 'numpy_array' - - Returns - ------- - tilde_f: - 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 : + # The space of the dual dofs + + # f : + # The function used for evaluation + + # backend_language: + # The backend used to accelerate the code + + # return_format: + # The format of the dofs, can be 'stencil_array' or 'numpy_array' + + # Returns + # ------- + # tilde_f: + # 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): @@ -280,6 +323,13 @@ 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)] diff --git a/psydac/feec/multipatch/examples/ppc_test_cases.py b/psydac/feec/multipatch/examples/ppc_test_cases.py index 122a68dca..d32485cce 100644 --- a/psydac/feec/multipatch/examples/ppc_test_cases.py +++ b/psydac/feec/multipatch/examples/ppc_test_cases.py @@ -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 @@ -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 diff --git a/psydac/feec/multipatch/fem_linear_operators.py b/psydac/feec/multipatch/fem_linear_operators.py index 4b8da2553..3982b9aa3 100644 --- a/psydac/feec/multipatch/fem_linear_operators.py +++ b/psydac/feec/multipatch/fem_linear_operators.py @@ -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 ): diff --git a/psydac/feec/multipatch/multipatch_domain_utilities.py b/psydac/feec/multipatch/multipatch_domain_utilities.py index 8ded13651..39b7286b8 100644 --- a/psydac/feec/multipatch/multipatch_domain_utilities.py +++ b/psydac/feec/multipatch/multipatch_domain_utilities.py @@ -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), @@ -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 diff --git a/psydac/feec/multipatch/non_matching_operators.py b/psydac/feec/multipatch/non_matching_operators.py index 461509646..909f1d10b 100644 --- a/psydac/feec/multipatch/non_matching_operators.py +++ b/psydac/feec/multipatch/non_matching_operators.py @@ -15,40 +15,40 @@ from psydac.fem.splines import SplineSpace from psydac.utilities.quadratures import gauss_legendre from psydac.core.bsplines import quadrature_grid, basis_ders_on_quad_grid, find_spans, elements_spans, cell_index, basis_ders_on_irregular_grid - - -def get_patch_index_from_face(domain, face): - """ - Return the patch index of subdomain/boundary - - Parameters - ---------- - domain : - The Symbolic domain - - face : - A patch or a boundary of a patch - - Returns - ------- - i : - The index of a subdomain/boundary in the multipatch domain - """ - - if domain.mapping: - domain = domain.logical_domain - if face.mapping: - face = face.logical_domain - - domains = domain.interior.args - if isinstance(face, Interface): - raise NotImplementedError( - "This face is an interface, it has several indices -- I am a machine, I cannot choose. Help.") - elif isinstance(face, Boundary): - i = domains.index(face.domain) - else: - i = domains.index(face) - return i +from psydac.feec.multipatch.operators import get_patch_index_from_face + +# def get_patch_index_from_face(domain, face): +# """ +# Return the patch index of subdomain/boundary + +# Parameters +# ---------- +# domain : +# The Symbolic domain + +# face : +# A patch or a boundary of a patch + +# Returns +# ------- +# i : +# The index of a subdomain/boundary in the multipatch domain +# """ + +# if domain.mapping: +# domain = domain.logical_domain +# if face.mapping: +# face = face.logical_domain + +# domains = domain.interior.args +# if isinstance(face, Interface): +# raise NotImplementedError( +# "This face is an interface, it has several indices -- I am a machine, I cannot choose. Help.") +# elif isinstance(face, Boundary): +# i = domains.index(face.domain) +# else: +# i = domains.index(face) +# return i class Local2GlobalIndexMap: @@ -583,7 +583,7 @@ def construct_h1_conforming_projection( # moment corrections perpendicular to interfaces # assume same moments everywhere gamma = get_1d_moment_correction( - Vh.spaces[0].spaces[0], p_moments=p_moments) + Vh.patch_space(0).spaces[0], p_moments=p_moments) domain = Vh.symbolic_space.domain ndim = 2 @@ -592,7 +592,7 @@ def construct_h1_conforming_projection( l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) for k in range(n_patches): - Vk = Vh.spaces[k] + Vk = Vh.patch_space(k) # T is a TensorFemSpace and S is a 1D SplineSpace shapes = [S.nbasis for S in Vk.spaces] l2g.set_patch_shapes(k, shapes) @@ -606,8 +606,8 @@ def construct_h1_conforming_projection( def get_vertex_index_from_patch(patch, coords): # coords = co[patch] - nbasis0 = Vh.spaces[patch].spaces[coords[0]].nbasis - 1 - nbasis1 = Vh.spaces[patch].spaces[coords[1]].nbasis - 1 + nbasis0 = Vh.patch_space(patch).spaces[coords[0]].nbasis - 1 + nbasis1 = Vh.patch_space(patch).spaces[coords[1]].nbasis - 1 # patch local index multi_index = [None] * ndim @@ -621,8 +621,8 @@ def vertex_moment_indices(axis, coords, patch, p_moments): if coords[axis] == 0: return range(1, p_moments + 2) else: - return range(Vh.spaces[patch].spaces[coords[axis]].nbasis - 1 - 1, - Vh.spaces[patch].spaces[coords[axis]].nbasis - 1 - p_moments - 2, -1) + return range(Vh.patch_space(patch).spaces[coords[axis]].nbasis - 1 - 1, + Vh.patch_space(patch).spaces[coords[axis]].nbasis - 1 - p_moments - 2, -1) # loop over all vertices for (bd, co) in corners.items(): @@ -772,8 +772,13 @@ def vertex_moment_indices(axis, coords, patch, p_moments): Proj_edge = sparse_eye(dim_tot, format="lil") Interfaces = domain.interfaces - if isinstance(Interfaces, Interface): - Interfaces = (Interfaces, ) + # why is this necessary -> improve Union in sympde ?? + if Interfaces is None: + iterable_Interfaces = () + elif isinstance(Interfaces, Interface): + iterable_Interfaces = (Interfaces, ) + else: + iterable_Interfaces = Interfaces def get_edge_index(j, axis, ext, space, k): multi_index = [None] * ndim @@ -822,7 +827,7 @@ def get_mu_minus(j, coarse_space, fine_space, R): return mu_minus # loop over all interfaces - for I in Interfaces: + for I in iterable_Interfaces: axis = I.axis direction = I.ornt # for now assume the interfaces are along the same direction @@ -830,8 +835,8 @@ def get_mu_minus(j, coarse_space, fine_space, R): k_minus = get_patch_index_from_face(domain, I.minus) k_plus = get_patch_index_from_face(domain, I.plus) - I_minus_ncells = Vh.spaces[k_minus].ncells - I_plus_ncells = Vh.spaces[k_plus].ncells + I_minus_ncells = Vh.patch_space(k_minus).ncells + I_plus_ncells = Vh.patch_space(k_plus).ncells # logical directions normal to interface if I_minus_ncells <= I_plus_ncells: @@ -848,8 +853,8 @@ def get_mu_minus(j, coarse_space, fine_space, R): d_fine = 1 - fine_axis d_coarse = 1 - coarse_axis - space_fine = Vh.spaces[k_fine] - space_coarse = Vh.spaces[k_coarse] + space_fine = Vh.patch_space(k_fine) + space_coarse = Vh.patch_space(k_coarse) coarse_space_1d = space_coarse.spaces[d_coarse] fine_space_1d = space_fine.spaces[d_fine] @@ -962,7 +967,7 @@ def get_mu_minus(j, coarse_space, fine_space, R): if hom_bc: for bn in domain.boundary: k = get_patch_index_from_face(domain, bn) - space_k = Vh.spaces[k] + space_k = Vh.patch_space(k) axis = bn.axis d = 1 - axis @@ -1037,8 +1042,10 @@ def construct_hcurl_conforming_projection( return sparse_eye(dim_tot, format="lil") # moment corrections perpendicular to interfaces + # gamma = [get_1d_moment_correction( + # Vh.spaces[0].spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] gamma = [get_1d_moment_correction( - Vh.spaces[0].spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] + Vh.patch_space(0).spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] domain = Vh.symbolic_space.domain ndim = 2 @@ -1047,7 +1054,7 @@ def construct_hcurl_conforming_projection( l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) for k in range(n_patches): - Vk = Vh.spaces[k] + Vk = Vh.patch_space(k) # T is a TensorFemSpace and S is a 1D SplineSpace shapes = [[S.nbasis for S in T.spaces] for T in Vk.spaces] l2g.set_patch_shapes(k, *shapes) @@ -1057,8 +1064,13 @@ def construct_hcurl_conforming_projection( Proj_edge = sparse_eye(dim_tot, format="lil") Interfaces = domain.interfaces - if isinstance(Interfaces, Interface): - Interfaces = (Interfaces, ) + # why is this necessary -> improve Union in sympde ?? + if Interfaces is None: + iterable_Interfaces = () + elif isinstance(Interfaces, Interface): + iterable_Interfaces = (Interfaces, ) + else: + iterable_Interfaces = Interfaces def get_edge_index(j, axis, ext, space, k): multi_index = [None] * ndim @@ -1075,7 +1087,7 @@ def edge_moment_index(p, i, axis, ext, space, k): return l2g.get_index(k, 1 - axis, multi_index) # loop over all interfaces - for I in Interfaces: + for I in iterable_Interfaces: direction = I.ornt # for now assume the interfaces are along the same direction assert direction == 1 @@ -1086,8 +1098,8 @@ def edge_moment_index(p, i, axis, ext, space, k): minus_axis, plus_axis = I.minus.axis, I.plus.axis # logical directions along the interface d_minus, d_plus = 1 - minus_axis, 1 - plus_axis - I_minus_ncells = Vh.spaces[k_minus].spaces[d_minus].ncells[d_minus] - I_plus_ncells = Vh.spaces[k_plus].spaces[d_plus].ncells[d_plus] + I_minus_ncells = Vh.patch_space(k_minus).spaces[d_minus].ncells[d_minus] + I_plus_ncells = Vh.patch_space(k_plus).spaces[d_plus].ncells[d_plus] # logical directions normal to interface if I_minus_ncells <= I_plus_ncells: @@ -1104,8 +1116,8 @@ def edge_moment_index(p, i, axis, ext, space, k): d_fine = 1 - fine_axis d_coarse = 1 - coarse_axis - space_fine = Vh.spaces[k_fine] - space_coarse = Vh.spaces[k_coarse] + space_fine = Vh.patch_space(k_fine) + space_coarse = Vh.patch_space(k_coarse) coarse_space_1d = space_coarse.spaces[d_coarse].spaces[d_coarse] fine_space_1d = space_fine.spaces[d_fine].spaces[d_fine] @@ -1164,7 +1176,7 @@ def edge_moment_index(p, i, axis, ext, space, k): # boundary condition for bn in domain.boundary: k = get_patch_index_from_face(domain, bn) - space_k = Vh.spaces[k] + space_k = Vh.patch_space(k) axis = bn.axis if not hom_bc: diff --git a/psydac/feec/multipatch/operators.py b/psydac/feec/multipatch/operators.py index 87368a1eb..b8e9fe6d4 100644 --- a/psydac/feec/multipatch/operators.py +++ b/psydac/feec/multipatch/operators.py @@ -37,7 +37,8 @@ def get_patch_index_from_face(domain, face): - """ Return the patch index of subdomain/boundary + """ + Return the patch index of subdomain/boundary Parameters ---------- @@ -58,7 +59,8 @@ def get_patch_index_from_face(domain, face): if face.mapping: face = face.logical_domain - domains = domain.interior.args + domains = domain.subdomains + if isinstance(face, Interface): raise NotImplementedError( "This face is an interface, it has several indices -- I am a machine, I cannot choose. Help.") @@ -70,7 +72,8 @@ def get_patch_index_from_face(domain, face): def get_interface_from_corners(corner1, corner2, domain): - """ Return the interface between two corners from two different patches that correspond to a single (physical) vertex. + """ + Return the interface between two corners from two different patches that correspond to a single (physical) vertex. Parameters ---------- @@ -968,6 +971,10 @@ def __init__( self._backend_language = backend_language self._dual_Hodge_sparse_matrix = None + domain = Vh.symbolic_space.domain + self._singlepatch = (len(domain) == 1) + print("self._singlepatch = ", self._singlepatch) + assert metric == 'identity' self._metric = metric @@ -1035,7 +1042,7 @@ def assemble_primal_Hodge_matrix(self): assert Vh == self.fem_codomain V = Vh.symbolic_space - domain = V.domain + domain = V.domain # domain_h = V0h.domain: would be nice... u, v = elements_of(V, names='u, v') @@ -1068,17 +1075,22 @@ def assemble_dual_Hodge_matrix(self): self.assemble_primal_Hodge_matrix() M = self._matrix # mass matrix of the (primal) basis - nrows = M.n_block_rows - ncols = M.n_block_cols + + if self._singlepatch: + inv_M = inv(self._sparse_matrix.tocsc()) + inv_M.eliminate_zeros() - inv_M_blocks = [] - for i in range(nrows): - Mii = M[i, i].tosparse() - inv_Mii = inv(Mii.tocsc()) - inv_Mii.eliminate_zeros() - inv_M_blocks.append(inv_Mii) + else: + nrows = M.n_block_rows + ncols = M.n_block_cols + inv_M_blocks = [] + for i in range(nrows): + Mii = M[i, i].tosparse() + inv_Mii = inv(Mii.tocsc()) + inv_Mii.eliminate_zeros() + inv_M_blocks.append(inv_Mii) + inv_M = block_diag(inv_M_blocks) - inv_M = block_diag(inv_M_blocks) self._dual_Hodge_sparse_matrix = inv_M # ============================================================================== diff --git a/psydac/feec/multipatch/plotting_utilities.py b/psydac/feec/multipatch/plotting_utilities.py index af522c6f1..83271ffca 100644 --- a/psydac/feec/multipatch/plotting_utilities.py +++ b/psydac/feec/multipatch/plotting_utilities.py @@ -10,12 +10,12 @@ from collections import OrderedDict from psydac.linalg.utilities import array_to_psydac -from psydac.fem.basic import FemField +from psydac.fem.basic import FemField, FemSpace from psydac.utilities.utils import refine_array_1d from psydac.feec.pull_push import push_2d_h1, push_2d_hcurl, push_2d_hdiv, push_2d_l2 __all__ = ( - 'is_vector_valued', + # 'is_vector_valued', 'get_grid_vals', 'get_grid_quad_weights', 'get_plotting_grid', @@ -28,11 +28,11 @@ # ============================================================================== -def is_vector_valued(u): - # small utility function, only tested for FemFields in multi-patch spaces of the 2D grad-curl sequence - # todo: a proper interface returning the number of components of a general - # FemField would be nice - return u.fields[0].space.is_product +# def is_vector_valued(u: FemField) -> bool: +# # small utility function, only tested for FemFields in multi-patch spaces of the 2D grad-curl sequence +# # todo: a proper interface returning the number of components of a general +# # FemField would be nice +# return u.fields[0].space.is_product # ------------------------------------------------------------------------------ @@ -40,13 +40,25 @@ def is_vector_valued(u): def get_grid_vals(u, etas, mappings_list, space_kind='hcurl'): """ get the physical field values, given the logical field and the logical grid - :param u: FemField + :param u: FemField or callable function :param etas: logical grid :param space_kind: specifies the push-forward for the physical values """ n_patches = len(mappings_list) - vector_valued = is_vector_valued(u) if isinstance( - u, FemField) else isinstance(u[0], (list, tuple)) + print(f' n_patches = {n_patches}') + + if isinstance(u, FemField): + vector_valued = u.is_vector_valued + else: + # then u should be callable + vector_valued = isinstance(u, (list, tuple)) # (MCP 4.02.25) before, it was isinstance(u[0], (list, tuple)) -- this needs to be tested + + # print(f'is_vector_valued(u) = {is_vector_valued(u)}') + # print(f'isinstance(u, FemField) = {isinstance(u, FemField)}') + # print(f'isinstance(u[0], (list, tuple)) = {isinstance(u[0], (list, tuple))}') + + print(f'vector_valued = {vector_valued}') + if vector_valued: # WARNING: here we assume 2D ! u_vals_components = [n_patches * [None], n_patches * [None]] @@ -60,11 +72,11 @@ def get_grid_vals(u, etas, mappings_list, space_kind='hcurl'): uk_field_1 = None if isinstance(u, FemField): if vector_valued: - uk_field_0 = u[k].fields[0] - uk_field_1 = u[k].fields[1] + uk_field_0 = u.patch_field(k).fields[0] + uk_field_1 = u.patch_field(k).fields[1] else: # it would be nice to just write u[k].fields[0] here... - uk_field_0 = u.fields[k] + uk_field_0 = u.patch_field(k) else: # then u should be callable if vector_valued: @@ -242,8 +254,8 @@ def get_patch_knots_gridlines(Vh, N, mappings, plotted_patch=-1): F = [M.get_callable_mapping() for d, M in mappings.items()] if plotted_patch in range(len(mappings)): - grid_x1 = Vh.spaces[plotted_patch].spaces[0].breaks[0] - grid_x2 = Vh.spaces[plotted_patch].spaces[0].breaks[1] + grid_x1 = Vh.patch_space(plotted_patch).spaces[0].breaks[0] + grid_x2 = Vh.patch_space(plotted_patch).spaces[0].breaks[1] x1 = refine_array_1d(grid_x1, N) x2 = refine_array_1d(grid_x2, N) @@ -303,6 +315,7 @@ def plot_field( raise ValueError( 'invalid value for space_kind = {}'.format(space_kind)) + print(f'type(numpy_coeffs) = {type(numpy_coeffs)}, numpy_coeffs.shape = {numpy_coeffs.shape}, Vh.nbasis = {Vh.nbasis}') vh = fem_field if vh is None: if numpy_coeffs is not None: @@ -310,8 +323,9 @@ def plot_field( stencil_coeffs = array_to_psydac(numpy_coeffs, Vh.vector_space) vh = FemField(Vh, coeffs=stencil_coeffs) - mappings = OrderedDict([(P.logical_domain, P.mapping) - for P in domain.interior]) + mappings = domain.mappings + # mappings = OrderedDict([(P.logical_domain, P.mapping) + # for P in domain.interior]) mappings_list = list(mappings.values()) etas, xx, yy = get_plotting_grid(mappings, N=N_vis) @@ -319,13 +333,13 @@ def grid_vals(v): return get_grid_vals( v, etas, mappings_list, space_kind=space_kind) vh_vals = grid_vals(vh) - if plot_type == 'vector_field' and not is_vector_valued(vh): + if plot_type == 'vector_field' and not vh.is_vector_valued: print( "WARNING [plot_field]: vector_field plot is not possible with a scalar field, plotting the amplitude instead") plot_type = 'amplitude' if plot_type == 'vector_field': - if is_vector_valued(vh): + if vh.is_vector_valued: my_small_streamplot( title=title, vals_x=vh_vals[0], @@ -343,7 +357,7 @@ def grid_vals(v): return get_grid_vals( # computing plot_vals_list: may have several elements for several plots if plot_type == 'amplitude': - if is_vector_valued(vh): + if vh.is_vector_valued: # then vh_vals[d] contains the values of the d-component of vh # (as a patch-indexed list) plot_vals = [np.sqrt(abs(v[0])**2 + abs(v[1])**2) @@ -355,7 +369,7 @@ def grid_vals(v): return get_grid_vals( plot_vals_list = [plot_vals] elif plot_type == 'components': - if is_vector_valued(vh): + if vh.is_vector_valued: # then vh_vals[d] contains the values of the d-component of vh # (as a patch-indexed list) plot_vals_list = vh_vals @@ -368,15 +382,23 @@ def grid_vals(v): return get_grid_vals( else: raise ValueError(plot_type) + if isinstance(Vh, FemSpace): + N_gl=20 + gridlines_x1, gridlines_x2 = get_patch_knots_gridlines(Vh, N_gl, mappings, plotted_patch=0) + else: + gridlines_x1, gridlines_x2 = None, None + my_small_plot( title=title, vals=plot_vals_list, titles=subtitles, xx=xx, yy=yy, + gridlines_x1=gridlines_x1, + gridlines_x2=gridlines_x2, surface_plot=surface_plot, cb_min=cb_min, - cb_max=cb_max, + cb_max=cb_max, save_fig=filename, save_vals=False, hide_plot=hide_plot, diff --git a/psydac/fem/basic.py b/psydac/fem/basic.py index d382b0bcf..2de7523ad 100644 --- a/psydac/fem/basic.py +++ b/psydac/fem/basic.py @@ -63,9 +63,32 @@ def is_product( self ): """ Boolean flag that describes whether the space is a product space. If True, an element of this space can be decomposed into separate fields. + + """ + + @property + @abstractmethod + def is_multipatch( self ): + """ + Boolean flag that describes whether the space is a multi-patch space. + + """ + + @property + @abstractmethod + def is_vector_valued( self ): + """ + Boolean flag that describes whether the space is vector-valued. """ + @abstractmethod + def patch_space(self, i): + """ + Return the space of the i-th patch. + (If single-patch space, return self.) + """ + @property @abstractmethod def symbolic_space( self ): @@ -231,7 +254,7 @@ def __init__( self, space, coeffs=None ): coeffs = space.vector_space.zeros() # Case of a vector field, element of a ProductSpace - if space.is_product: + if space.is_product: # (MCP 4.02.25): this should be specified as 'if space.is_vector or space.is_multipatch' fields = tuple(FemField(V, c) for V, c in zip(space.spaces, coeffs)) else: fields = tuple() @@ -265,6 +288,17 @@ def coeffs( self ): def fields(self): return self._fields + @property + def is_vector_valued(self): + return self.space.is_vector_valued + + def patch_field(self, i): + """ return field of patch i (if multipatch), or self if single patch """ + if self.space.is_multipatch: + return self.fields[i] + else: + return self + # ... def __getitem__(self, key): return self._fields[key] diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index cf1a284dd..def1a55fe 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,8 +1,17 @@ import numpy as np -from psydac.linalg.kron import KroneckerDenseMatrix -from psydac.core.bsplines import hrefinement_matrix -from psydac.linalg.stencil import StencilVectorSpace +from sympde.topology import element_of +from sympde.topology.space import ScalarFunction +from sympde.topology.mapping import Mapping +from sympde.calculus import dot +from sympde.expr.expr import LinearForm, integral + +from psydac.api.settings import PSYDAC_BACKENDS + +from psydac.linalg.kron import KroneckerDenseMatrix +from psydac.core.bsplines import hrefinement_matrix +from psydac.linalg.stencil import StencilVectorSpace +from psydac.fem.basic import FemSpace __all__ = ('knots_to_insert', 'knot_insertion_projection_operator') @@ -100,3 +109,63 @@ def knot_insertion_projection_operator(domain, codomain): ops.append(np.eye(d.nbasis)) return KroneckerDenseMatrix(domain.vector_space, codomain.vector_space, *ops) + + +def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='stencil_array'): + """ + return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(Vh)) of a given function f, as a stencil array or numpy array + + Parameters + ---------- + Vh : FemSpace + The discrete space for the dual dofs + + f : + The function used for evaluation + + domain_h : + The discrete domain corresponding to Vh + + backend_language: + The backend used to accelerate the code + + return_format: + The format of the dofs, can be 'stencil_array' or 'numpy_array' + + Returns + ------- + tilde_f: + The dual dofs + """ + + from psydac.api.discretization import discretize + + assert isinstance(Vh, FemSpace) + # 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') + + print(f'Vh.is_vector_valued = {Vh.is_vector_valued} ***** isinstance(v, ScalarFunction) = {isinstance(v, ScalarFunction)}') + if Vh.is_vector_valued: + expr = dot(f,v) + else: + #isinstance(v, ScalarFunction): + expr = f*v + + l = LinearForm(v, integral( V.domain, expr)) + lh = discretize(l, 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 + diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 11ebbc1f2..14e99ca87 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -269,6 +269,18 @@ def vector_space( self ): def is_product(self): return False + @property + def is_multipatch(self): + return False + + @property + def is_vector_valued(self): + return False + + def patch_space(self, i): + if i > 0: raise ValueError('Invalid patch index (> 0) for a single patch space') + return self + @property def symbolic_space( self ): return self._symbolic_space diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index e0367e43d..e616019f3 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -162,6 +162,18 @@ def vector_space(self): def is_product(self): return False + @property + def is_multipatch(self): + return False + + @property + def is_vector_valued(self): + return False + + def patch_space(self, i): + if i > 0: raise ValueError('Invalid patch index (> 0) for a single patch space') + return self + @property def symbolic_space( self ): return self._symbolic_space diff --git a/psydac/fem/vector.py b/psydac/fem/vector.py index 84a29db10..6c2f90169 100644 --- a/psydac/fem/vector.py +++ b/psydac/fem/vector.py @@ -94,6 +94,18 @@ def vector_space(self): def is_product(self): return True + @property + def is_multipatch(self): + return False + + @property + def is_vector_valued(self): + return True + + def patch_space(self, i): + if i > 0: raise ValueError('Invalid patch index (> 0) for a single patch space') + return self + @property def symbolic_space( self ): return self._symbolic_space @@ -415,6 +427,17 @@ def vector_space(self): def is_product(self): return True + @property + def is_multipatch(self): + return True + + @property + def is_vector_valued(self): + return self.patch_space(0).is_vector_valued + + def patch_space(self, i): + return self.spaces[i] + @property def symbolic_space( self ): return self._symbolic_space