diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 41896bbb9..35dafac3f 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -36,9 +36,10 @@ jobs: run: | sudo apt update sudo apt install graphviz - - name: Install Python dependencies + - name: Install Python dependencies (including PSYDAC) run: | - python -m pip install -r docs/requirements.txt + pip install . + pip install -r docs/requirements.txt - name: Make the sphinx doc run: | rm -rf docs/source/modules/STUBDIR diff --git a/docs/source/modules/feec.multipatch.rst b/docs/source/modules/feec.multipatch.rst index 64bf49a12..0ab0f32ad 100644 --- a/docs/source/modules/feec.multipatch.rst +++ b/docs/source/modules/feec.multipatch.rst @@ -12,6 +12,5 @@ feec.multipatch multipatch.multipatch_domain_utilities multipatch.non_matching_operators multipatch.operators - multipatch.plotting_utilities multipatch.utilities multipatch.utils_conga_2d diff --git a/psydac/cad/tests/test_geometry.py b/psydac/cad/tests/test_geometry.py index 0f0f1a93d..7a75c72ca 100644 --- a/psydac/cad/tests/test_geometry.py +++ b/psydac/cad/tests/test_geometry.py @@ -222,6 +222,11 @@ def test_export_nurbs_to_hdf5(ncells, degree): assert np.allclose(pcoords1[..., :domain.dim], pcoords2, 1e-15, 1e-15) + jacobian1 = new_pipe.gradient(u=eta1, v=eta2)[..., :domain.dim, :] + jacobian2 = np.array([[mapping.jacobian(e1, e2) for e2 in eta2] for e1 in eta1]) + + assert np.allclose(jacobian1, jacobian2, 1e-13, 1e-13) + #============================================================================== @pytest.mark.parametrize( 'ncells', [[8,8], [12,12], [14,14]] ) @pytest.mark.parametrize( 'degree', [[2,2], [3,2], [2,3], [3,3], [4,4]] ) diff --git a/psydac/core/field_evaluation_kernels.py b/psydac/core/field_evaluation_kernels.py index e75854d98..e2e69210d 100644 --- a/psydac/core/field_evaluation_kernels.py +++ b/psydac/core/field_evaluation_kernels.py @@ -1,10 +1,239 @@ import numpy as np from typing import TypeVar +__all__ = ( + 'eval_field_1d_once', + 'eval_field_2d_once', + 'eval_field_3d_once', + 'eval_fields_1d_irregular_no_weights', + 'eval_fields_1d_irregular_weighted', + 'eval_fields_1d_no_weights', + 'eval_fields_1d_weighted', + 'eval_fields_2d_irregular_no_weights', + 'eval_fields_2d_irregular_weighted', + 'eval_fields_2d_no_weights', + 'eval_fields_2d_weighted', + 'eval_fields_3d_irregular_no_weights', + 'eval_fields_3d_irregular_weighted', + 'eval_fields_3d_no_weights', + 'eval_fields_3d_weighted', + 'eval_jac_det_2d', + 'eval_jac_det_2d_weights', + 'eval_jac_det_3d', + 'eval_jac_det_3d_weights', + 'eval_jac_det_irregular_2d', + 'eval_jac_det_irregular_2d_weights', + 'eval_jac_det_irregular_3d', + 'eval_jac_det_irregular_3d_weights', + 'eval_jacobians_2d', + 'eval_jacobians_2d_weights', + 'eval_jacobians_3d', + 'eval_jacobians_3d_weights', + 'eval_jacobians_inv_2d', + 'eval_jacobians_inv_2d_weights', + 'eval_jacobians_inv_3d', + 'eval_jacobians_inv_3d_weights', + 'eval_jacobians_inv_irregular_2d', + 'eval_jacobians_inv_irregular_2d_weights', + 'eval_jacobians_inv_irregular_3d', + 'eval_jacobians_inv_irregular_3d_weights', + 'eval_jacobians_irregular_2d', + 'eval_jacobians_irregular_2d_weights', + 'eval_jacobians_irregular_3d', + 'eval_jacobians_irregular_3d_weights', + 'pushforward_2d_hcurl', + 'pushforward_2d_hdiv', + 'pushforward_2d_l2', + 'pushforward_3d_hcurl', + 'pushforward_3d_hdiv', + 'pushforward_3d_l2', +) + T = TypeVar('T', float, complex) + # ============================================================================= # Field evaluation functions # ============================================================================= + +# ----------------------------------------------------------------------------- +# 0: Evaluation of single 3d field at single point +# ----------------------------------------------------------------------------- +def eval_field_3d_once(local_coeffs : 'T[:,:,:]', + local_bases_1: 'float[:]', + local_bases_2: 'float[:]', + local_bases_3: 'float[:]', + work_1d: 'T[:]', + work_2d: 'T[:,:]') -> T: + """ + Evaluate a scalar field in 3D, given the basis values along each direction. + + Compute the point value of a field, given the values of the non-zero 1D + basis functions at the point location. We assume that the 3D basis has a + tensor-product structure, hence each 3D basis function can be written as + the product of three 1D basis functions: + + phi(x1, x2, x3) = phi_1(x1) * phi_2(x2) * phi_3(x3). + + The algorithm uses sum factorization to compute the reduction. This needs + work arrays in 1D and 2D, which are provided to this function rather than + allocated at each function call. + + Parameters + ---------- + local_coeffs: ndarray[float|complex] + Active (local) coefficients of the field in all directions. + + local_bases_1: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x1. + + local_bases_2: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x2. + + local_bases_3: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x3. + + work_1d: ndarray[float|complex] + x1-dependent work array for the partial reduction in (x2, x3). + + work_2d: ndarray[float|complex] + (x1, x2)-dependent work array for the partial reduction in x3. + + Returns + ------- + res : float | complex + The scalar value of the field at the location of interest. + + """ + n1, n2, n3 = local_coeffs.shape + + # Verify that local bases have the correct size + assert local_bases_1.size == n1 + assert local_bases_2.size == n2 + assert local_bases_3.size == n3 + + # Verify that provided work arrays are large enough + assert work_1d.shape[0] >= n1 + assert work_2d.shape[0] >= n1 + assert work_2d.shape[1] >= n2 + + # Obtain zero of correct type (float or complex) from input array + c0 = local_coeffs[0, 0, 0] + zero = c0 - c0 + + # Compute reduction using sum factorization in 3D + res = zero + for i1 in range(n1): + work_1d[i1] = zero + for i2 in range(n2): + work_2d[i1, i2] = zero + for i3 in range(n3): + work_2d[i1, i2] += local_coeffs[i1, i2, i3] * local_bases_3[i3] + work_1d[i1] += work_2d[i1, i2] * local_bases_2[i2] + res += work_1d[i1] * local_bases_1[i1] + + return res + + +def eval_field_2d_once(local_coeffs : 'T[:,:]', + local_bases_1: 'float[:]', + local_bases_2: 'float[:]', + work_1d: 'T[:]') -> T: + """ + Evaluate a scalar field in 2D, given the basis values along each direction. + + Compute the point value of a field, given the values of the non-zero 1D + basis functions at the point location. We assume that the 2D basis has a + tensor-product structure, hence each 2D basis function can be written as + the product of three 1D basis functions: + + phi(x1, x2) = phi_1(x1) * phi_2(x2). + + The algorithm uses sum factorization to compute the reduction. This needs + a 1D work array, which is provided to this function rather than allocated + at each function call. + + Parameters + ---------- + local_coeffs: ndarray[float|complex] + Active (local) coefficients of the field in all directions. + + local_bases_1: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x1. + + local_bases_2: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x2. + + work_1d: ndarray[float|complex] + x1-dependent work array for the partial reduction in x2. + + Returns + ------- + res : float | complex + The scalar value of the field at the location of interest. + + """ + n1, n2 = local_coeffs.shape + + # Verify that local bases have the correct size + assert local_bases_1.size == n1 + assert local_bases_2.size == n2 + + # Verify that provided work arrays are large enough + assert work_1d.shape[0] >= n1 + + # Obtain zero of correct type (float or complex) from input array + c0 = local_coeffs[0, 0] + zero = c0 - c0 + + # Compute reduction using sum factorization in 2D + res = zero + for i1 in range(n1): + work_1d[i1] = zero + for i2 in range(n2): + work_1d[i1] += local_coeffs[i1, i2] * local_bases_2[i2] + res += work_1d[i1] * local_bases_1[i1] + + return res + + +def eval_field_1d_once(local_coeffs : 'T[:]', + local_bases_1: 'float[:]') -> T: + """ + Evaluate a scalar field in 1D, given the basis values. + + Compute the point value of a field, given the values of the non-zero 1D + basis functions at the point location. + + Parameters + ---------- + local_coeffs: ndarray[float|complex] + Active (local) coefficients of the field in all directions. + + local_bases_1: ndarray[float]. + Active (local) 1D-basis functions values at the point of evaluation x1. + + Returns + ------- + res : float | complex + The scalar value of the field at the location of interest. + + """ + n1, = local_coeffs.shape + + # Verify that local bases have the correct size + assert local_bases_1.size == n1 + + # Obtain zero of correct type (float or complex) from input array + c0 = local_coeffs[0] + zero = c0 - c0 + + # Compute reduction as a simple scalar product + res = zero + for i1 in range(n1): + res += local_coeffs[i1] * local_bases_1[i1] + + return res + # ----------------------------------------------------------------------------- # 1: Regular tensor grid without weight # ----------------------------------------------------------------------------- diff --git a/psydac/core/tests/test_field_evaluation_kernel.py b/psydac/core/tests/test_field_evaluation_kernel.py index 09e1d666e..2702a367c 100644 --- a/psydac/core/tests/test_field_evaluation_kernel.py +++ b/psydac/core/tests/test_field_evaluation_kernel.py @@ -411,17 +411,15 @@ def test_regular_evaluations(knots, ldim, degree, npts_per_cell): f_direct = np.array([space_h.eval_fields([e1], field) for e1 in regular_grid[0]]) # Weighted - f_direct_w = np.array([np.array(space_h.eval_fields([e1], field, weights=weight)) - / np.array(space_h.eval_fields([e1], weight)) - for e1 in regular_grid[0]]) + f_direct_w = np.array([space_h.eval_fields([e1], field, weights=weight) + for e1 in regular_grid[0]]) if ldim == 2: # No weights f_direct = np.array([[space_h.eval_fields([e1, e2], field) for e2 in regular_grid[1]] for e1 in regular_grid[0]]) # Weighted - f_direct_w = np.array([[np.array(space_h.eval_fields([e1, e2], field, weights=weight)) - / np.array(space_h.eval_fields([e1, e2], weight)) + f_direct_w = np.array([[space_h.eval_fields([e1, e2], field, weights=weight) for e2 in regular_grid[1]] for e1 in regular_grid[0]]) @@ -433,8 +431,7 @@ def test_regular_evaluations(knots, ldim, degree, npts_per_cell): for e1 in regular_grid[0]]) # Weighted - f_direct_w = np.array([[[np.array(space_h.eval_fields([e1, e2, e3], field, weights=weight)) - / np.array(space_h.eval_fields([e1, e2, e3], weight)) + f_direct_w = np.array([[[space_h.eval_fields([e1, e2, e3], field, weights=weight) for e3 in regular_grid[2]] for e2 in regular_grid[1]] for e1 in regular_grid[0]]) @@ -571,17 +568,15 @@ def test_irregular_evaluations(knots, ldim, degree, npts): f_direct = np.array([space_h.eval_fields([e1], field) for e1 in irregular_grid[0]]) # Weighted - f_direct_w = np.array([np.array(space_h.eval_fields([e1], field, weights=weight)) - / np.array(space_h.eval_fields([e1], weight)) - for e1 in irregular_grid[0]]) + f_direct_w = np.array([space_h.eval_fields([e1], field, weights=weight) + for e1 in irregular_grid[0]]) if ldim == 2: # No weights f_direct = np.array([[space_h.eval_fields([e1, e2], field) for e2 in irregular_grid[1]] for e1 in irregular_grid[0]]) # Weighted - f_direct_w = np.array([[np.array(space_h.eval_fields([e1, e2], field, weights=weight)) - / np.array(space_h.eval_fields([e1, e2], weight)) + f_direct_w = np.array([[space_h.eval_fields([e1, e2], field, weights=weight) for e2 in irregular_grid[1]] for e1 in irregular_grid[0]]) @@ -593,8 +588,7 @@ def test_irregular_evaluations(knots, ldim, degree, npts): for e1 in irregular_grid[0]]) # Weighted - f_direct_w = np.array([[[np.array(space_h.eval_fields([e1, e2, e3], field, weights=weight)) - / np.array(space_h.eval_fields([e1, e2, e3], weight)) + f_direct_w = np.array([[[space_h.eval_fields([e1, e2, e3], field, weights=weight) for e3 in irregular_grid[2]] for e2 in irregular_grid[1]] for e1 in irregular_grid[0]]) @@ -692,3 +686,7 @@ def test_pushforwards_hcurl(ldim): pushforward_3d_hcurl(field_to_push, inv_jacobians, out) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) + +if __name__ == '__main__': + test_regular_jacobians('bent_pipe.h5', 2) + test_regular_evaluations([np.sort(np.concatenate((np.zeros(3), np.random.random(9), np.ones(3)))) for i in range(2)], 2, [2] * 2,2) diff --git a/psydac/fem/basic.py b/psydac/fem/basic.py index 09a1466fb..11a5224ad 100644 --- a/psydac/fem/basic.py +++ b/psydac/fem/basic.py @@ -19,7 +19,6 @@ class FemSpace( metaclass=ABCMeta ): Generic Finite Element space V. A unique basis is associated to a FemSpace, i.e. FemSpace = Span( basis ) - """ #----------------------------------------- @@ -99,12 +98,11 @@ def axis_spaces(self): Return the axis spaces (self if univariate) as a tuple. """ - #--------------------------------------- # Abstract interface: evaluation methods #--------------------------------------- @abstractmethod - def eval_field( self, field, *eta, weights=None): + def eval_field(self, field, *eta, weights=None): """ Evaluate field at location(s) eta. @@ -113,21 +111,20 @@ def eval_field( self, field, *eta, weights=None): field : FemField Field object (element of FemSpace) to be evaluated. - eta : list of float or numpy.ndarray + eta : float | numpy.ndarray[float] Evaluation point(s) in logical domain. - weights : StencilVector, optional + weights : Vector, optional Weights of the basis functions, such that weights.space == field.coeffs.space. Returns ------- - value : float or numpy.ndarray + value : float | complex | ndarray[float] | ndarray[complex] Field value(s) at location(s) eta. - """ @abstractmethod - def eval_field_gradient( self, field, *eta , weights=None): + def eval_field_gradient(self, field, *eta, weights=None): """ Evaluate field gradient at location(s) eta. @@ -139,17 +136,15 @@ def eval_field_gradient( self, field, *eta , weights=None): eta : list of float or numpy.ndarray Evaluation point(s) in logical domain. - weights : StencilVector, optional + weights : Vector, optional Weights of the basis functions, such that weights.space == field.coeffs.space. Returns ------- - value : float or numpy.ndarray + value : ndarray[float] | ndarray[complex] Value(s) of field gradient at location(s) eta. - """ - #---------------------- # Concrete methods #---------------------- @@ -252,7 +247,6 @@ class FemField: coeffs : psydac.linalg.basic.Vector (optional) Vector of coefficients in finite element basis (by default assume zero vector). - """ def __init__( self, space, coeffs=None ): @@ -290,7 +284,6 @@ def coeffs( self ): Coefficients are stored into one element of the vector space in 'self.space.coeff_space', which is topologically associated to the finite element space. - """ return self._coeffs @@ -320,14 +313,14 @@ def __getitem__(self, key): return self._fields[key] # ... - def __call__( self, *eta , weights=None): + def __call__(self, *eta , weights=None): """Evaluate weighted field at location identified by logical coordinates eta.""" - return self._space.eval_field( self, *eta , weights=weights) + return self._space.eval_field(self, *eta , weights=weights) # ... - def gradient( self, *eta , weights=None): + def gradient(self, *eta , weights=None): """Evaluate gradient of weighted field at location identified by logical coordinates eta.""" - return self._space.eval_field_gradient( self, *eta , weights=weights) + return self._space.eval_field_gradient(self, *eta , weights=weights) # ... def divergence(self, *eta, weights=None): diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 644f08cc4..112645ecc 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -9,7 +9,6 @@ import numpy as np import itertools import h5py -import os from types import MappingProxyType @@ -23,14 +22,15 @@ from psydac.fem.partitioning import create_cart, partition_coefficients from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.core.bsplines import (find_span, - basis_funs, - basis_funs_1st_der, - basis_ders_on_quad_grid, +from psydac.core.bsplines import (basis_ders_on_quad_grid, elements_spans, cell_index, basis_ders_on_irregular_grid) +from psydac.core.bsplines_kernels import (find_span_p, + basis_funs_p, + basis_funs_1st_der_p) + from psydac.core.field_evaluation_kernels import (eval_fields_1d_no_weights, eval_fields_1d_irregular_no_weights, eval_fields_1d_weighted, @@ -42,10 +42,70 @@ eval_fields_3d_no_weights, eval_fields_3d_irregular_no_weights, eval_fields_3d_weighted, - eval_fields_3d_irregular_weighted) + eval_fields_3d_irregular_weighted, + eval_field_3d_once, + eval_field_2d_once, + eval_field_1d_once) __all__ = ('TensorFemSpace',) +#=============================================================================== +def eval_field_nd_once(coeffs, *args): + """ + Evaluate a scalar field of a tensor-product space in any number of dimensions. + + Uses sum factorization for greater efficiency. Does not allocate temporary + arrays, instead uses work arrays given as function arguments. + + Parameters + ---------- + coeffs : ndarray[float|complex] + Active (local) degrees of freedom of the field, i.e. the coefficients + of the non-zero N-dimensional basis functions, given as an N-dimensional + array. + + *args : tuple[Any] + Essentially, args = (*bases, *works). See Notes. + + Notes + ----- + bases : tuple[ndarray[float]] + Active (local) 1D basis functions along each of the N dimensions. + + works : tuple[ndarray[float|complex]] + N-1 work arrays of increasing rank (from 1 to N-1) for partial reductions. + """ + # Extract useful information from the coeffs array + ldim = coeffs.ndim + dtype = coeffs.dtype + shape = coeffs.shape + + # Extract values of 1D basis functions along each direction and check arrays + bases = args[:ldim] + + for basis, n in zip(bases, coeffs.shape): + assert basis.ndim == 1 + assert basis.size == n + assert basis.dtype == float + + # Extract work arrays and check them + works = (np.array(0, dtype=dtype), *args[ldim:]) + + assert len(works) == ldim + for i, work in enumerate(works): + ndim = work.ndim + assert ndim == i + assert all(s >= n for s, n in zip(work.shape, shape[:ndim])) + assert work.dtype == dtype + + # Compute reduction using sum factorization in N dimensions + res = coeffs + for basis, work in zip(bases[::-1], works[::-1]): + np.dot(res, basis, out=work) + res = work + + return res.item() + #=============================================================================== class TensorFemSpace(FemSpace): """ @@ -96,6 +156,7 @@ def __init__(self, domain_decomposition, *spaces, coeff_space=None, cart=None, d coeff_space = StencilVectorSpace(cart, dtype=dtype) # Store some info + self._ldim = sum(V.ldim for V in spaces) self._domain_decomposition = domain_decomposition self._spaces = spaces self._dtype = dtype @@ -132,6 +193,21 @@ def __init__(self, domain_decomposition, *spaces, coeff_space=None, cart=None, d # Store information about nested grids self.set_refined_space(self.ncells, self) + # Select evaluation function based on dimensionality + if self.ldim == 1: + self.eval_field_once = eval_field_1d_once + elif self.ldim == 2: + self.eval_field_once = eval_field_2d_once + elif self.ldim == 3: + self.eval_field_once = eval_field_3d_once + else: + self.eval_field_once = eval_field_nd_once + + # Create work arrays based on dimensionality + nbasis = [p+1 for p in self.degree] + work = [np.zeros(nbasis[:d], dtype=dtype) for d in range(1, self.ldim)] + self._work = tuple(work) + #-------------------------------------------------------------------------- # Abstract interface: read-only attributes #-------------------------------------------------------------------------- @@ -139,7 +215,7 @@ def __init__(self, domain_decomposition, *spaces, coeff_space=None, cart=None, d def ldim(self): """ Parametric dimension. """ - return sum([V.ldim for V in self.spaces]) + return self._ldim @property def periodic(self): @@ -153,10 +229,6 @@ def periodic(self): # spaces in self.spaces seem to be returning a single scalar value. return tuple(V.periodic for V in self.spaces) - @property - def domain_decomposition(self): - return self._domain_decomposition - @property def mapping(self): # [YG, 28.03.2025]: not clear why there should be no mapping here... @@ -172,15 +244,19 @@ def coeff_space(self): return self._coeff_space @property - def symbolic_space( self ): - return self._symbolic_space + def is_multipatch(self): + return False @property - def interfaces( self ): - return self._interfaces_readonly + def is_vector_valued(self): + return False + + @property + def symbolic_space(self): + return self._symbolic_space @symbolic_space.setter - def symbolic_space( self, symbolic_space ): + def symbolic_space(self, symbolic_space): assert isinstance(symbolic_space, BasicFunctionSpace) self._symbolic_space = symbolic_space @@ -196,37 +272,111 @@ def component_spaces(self): def axis_spaces(self): return self._spaces + #-------------------------------------------------------------------------- + # Other properties + #-------------------------------------------------------------------------- @property - def is_multipatch(self): - return False + def interfaces(self): + return self._interfaces_readonly @property - def is_vector_valued(self): - return False + def domain_decomposition(self): + return self._domain_decomposition #-------------------------------------------------------------------------- # Abstract interface: evaluation methods #-------------------------------------------------------------------------- - def eval_field( self, field, *eta, weights=None): + def eval_field(self, field, *eta, weights=None): + """ + Evaluate a single field at a single point location (eta1, eta2, ...). - assert isinstance( field, FemField ) - assert field.space is self - assert len( eta ) == self.ldim - if weights: - assert weights.space == field.coeffs.space + Parameters + ---------- + field : FemField + The field to be evaluated, element of this space (`self`). - bases = [] - index = [] + *eta : Iterable[float] + The logical coordinates of the point of evaluation. - # Necessary if vector coeffs is distributed across processes - if not field.coeffs.ghost_regions_in_sync: - field.coeffs.update_ghost_regions() + weights : Vector + The weights of `field`, necessary if it is a NURBS (optional). + They must belong to the coefficient space of `self`. + + Returns + ------- + float | complex + The value of the field at the point of interest. + """ + values_list = self.eval_fields_one_point([field], *eta, weights=weights) + + return values_list[0] + + # ... + def eval_field_gradient(self, field, *eta, weights=None): + """ + Evaluate the gradient of a single field at the point (eta1, eta2, ...). + + Parameters + ---------- + field : FemField + The field to be evaluated, element of this space (`self`). + + *eta : Iterable[float] + The logical coordinates of the point of evaluation. - for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): + weights : FemField + The weights of `field`, necessary if it is a NURBS (optional). + Returns + ------- + tuple[float | complex] + The gradient of the field at the point of interest. + """ + gradients_list = self.eval_fields_gradient_one_point([field], *eta, weights=weights) + + return gradients_list[0] + + #-------------------------------------------------------------------------- + # Private utility methods + #-------------------------------------------------------------------------- + def _compute_index_and_bases(self, *eta, deriv=False, rational=False): + + assert len(eta) == self.ldim + assert isinstance(deriv, bool) + assert isinstance(rational, bool) + + # ------------------------------- + # Local storage for this function + # ------------------------------- + # [a] List of N 1D arrays w/ values of non-zero 1D basis functions + # along each direction + if not hasattr(self, '_tmp_bf'): + self._tmp_bf = [np.zeros(s.degree + 1, dtype=float) for s in self.spaces] + # + # [b] List of N 1D arrays w/ 1st derivative of non-zero 1D basis + # functions along each direction + if not hasattr(self, '_tmp_bf1d') and deriv: + self._tmp_bf1d = [np.zeros(s.degree + 1, dtype=float) for s in self.spaces] + # ------------------------------- + + # Shortcuts + bases_0 = self._tmp_bf + bases_1 = self._tmp_bf1d if deriv else [] + + index = [] + w_index = [] + + for i, (x, xlim, space) in enumerate(zip(eta, self.eta_lims, self.spaces)): + + # Shortcuts knots = space.knots degree = space.degree - span = find_span( knots, degree, x ) + + # All kernels expect the coordinate 'x' to be a Python float + fx = float(x) + + # Compute span, i.e. index of last non-zero B-spline + span = find_span_p(knots, degree, fx) #-------------------------------------------------# # Fix span for boundaries between subdomains # @@ -234,48 +384,193 @@ def eval_field( self, field, *eta, weights=None): # TODO: Use local knot sequence instead of global # # one to get correct span in all situations # #-------------------------------------------------# - if x == xlim[1] and x != knots[-1-degree]: + if fx == xlim[1] and fx != knots[-1 - degree]: span -= 1 #-------------------------------------------------# - basis = basis_funs( knots, degree, x, span) + + # Compute values of non-zero B-splines, and 1st derivative if required + basis_funs_p(knots, degree, fx, span, out = bases_0[i]) + if deriv: + basis_funs_1st_der_p(knots, degree, fx, span, out = bases_1[i]) # If needed, rescale B-splines to get M-splines if space.basis == 'M': - basis *= space.scaling_array[span-degree : span+1] + scaling = space.scaling_array[span - degree : span + 1] + bases_0[i] *= scaling + if deriv: + bases_1[i] *= scaling # Determine local span - wrap_x = space.periodic and x > xlim[1] + wrap_x = space.periodic and fx > xlim[1] loc_span = span - space.nbasis if wrap_x else span - bases.append( basis ) - index.append( slice( loc_span-degree, loc_span+1 ) ) - # Get contiguous copy of the spline coefficients required for evaluation - index = tuple( index ) - coeffs = field.coeffs[index].copy() + # Collect indices of coefficients corresponding to non-zero basis functions + start = self.coeff_space.starts[i] + pad = self.coeff_space.pads [i] + shift = self.coeff_space.shifts[i] + index.append(slice(loc_span - degree - start + shift * pad, + loc_span + 1 - start + shift * pad)) + if rational: + w_index.append(slice(loc_span - degree, loc_span + 1)) + + return tuple(index), tuple(w_index), bases_0, bases_1 + + # ... + def eval_fields_one_point(self, fields, *eta, weights=None): + + for field in fields: + assert isinstance(field, FemField) + assert field.space is self + if weights: + assert weights.space == field.coeffs.space + if not field.coeffs.ghost_regions_in_sync: + # Necessary if vector coeffs is distributed across processes + field.coeffs.update_ghost_regions() + + assert len(eta) == self.ldim + + # ------------------------------- + # Local storage for this function + # ------------------------------- + # [a] One N-dim. array w/ coefficients of non-zero N-dim. tensor-product basis functions + if not hasattr(self, '_tmp_coeffs'): + self._tmp_coeffs = np.zeros([s.degree + 1 for s in self.spaces], dtype=self.dtype) + # + # [b] Same as [a], but for the weights' coefficients + if not hasattr(self, '_tmp_w_coeffs') and weights: + self._tmp_w_coeffs = np.zeros([s.degree + 1 for s in self.spaces], dtype=float) + # ------------------------------- + + # The following tuples contain, for each direction: + # + # bases: values of the non-zero 1D basis functions at coordinate eta_i + # index: indices (as slices) of the coefficients corresponding to + # the non-zero 1D basis functions. + # w_index: same as `index`, but for obtaining NURBS weights. + index, w_index, bases, _ = self._compute_index_and_bases(*eta, + rational=(weights is not None)) + + # If NURBS, evaluate weights field if weights: - coeffs *= weights[index] + w_coeffs = self._tmp_w_coeffs + w_coeffs[:] = weights[w_index] + w_value = self.eval_field_once(w_coeffs, *bases, *self._work) + + # Evaluate all fields at the same point + coeffs = self._tmp_coeffs + values_list = [] + for field in fields: + + # Get contiguous copy of the spline coefficients required for evaluation + coeffs[:] = field.coeffs._data[index] + + # If NURBS, rescale field coefficients by the weights' coefficients + if weights: + coeffs *= w_coeffs + + # Evaluate field + value = self.eval_field_once(coeffs, *bases, *self._work) + + # If NURBS, normalize result by the value of the weights field + if weights: + value /= w_value + + # Append field value to list + values_list.append(value) + + return values_list + + # ... + def eval_fields_gradient_one_point(self, fields, *eta, weights=None): + """ + Evaluate the gradient of multiple fields at the point (eta1, eta2, ...) + + Parameters + ---------- + field : FemField + The field to be evaluated, element of this space (`self`). - # Evaluation of multi-dimensional spline - # TODO: optimize + *eta : Iterable[float] + The logical coordinates of the point of evaluation. - # Option 1: contract indices one by one and store intermediate results - # - Pros: small number of Python iterations = ldim - # - Cons: we create ldim-1 temporary objects of decreasing size + weights : FemField + The weights of `field`, necessary if it is a NURBS (optional). + + Returns + ------- + list[tuple[float | complex]] + The gradient of all fields at the point of interest. + """ + for field in fields: + assert isinstance(field, FemField) + assert field.space is self + assert len(eta) == self.ldim + + # ------------------------------- + # Local storage for this function + # ------------------------------- + # [a] One N-dim. array w/ coefficients of non-zero N-dim. tensor-product basis functions + if not hasattr(self, '_tmp_coeffs'): + self._tmp_coeffs = np.zeros([s.degree + 1 for s in self.spaces], dtype=self.dtype) + # + # [b] Same as [a], but for the weights' coefficients + if not hasattr(self, '_tmp_w_coeffs') and weights: + self._tmp_w_coeffs = np.zeros([s.degree + 1 for s in self.spaces], dtype=float) + # ------------------------------- + + # The following tuples contain, for each direction: # - res = coeffs - for basis in bases[::-1]: - res = np.dot( res, basis ) + # bases_0: values of the non-zero 1D basis functions at coordinate eta_i + # bases_1: 1st derivative of the non-zero 1D basis functions at eta_i + # index: indices (as slices) of the coefficients corresponding to + # the non-zero 1D basis functions. + # w_index: same as `index`, but for obtaining NURBS weights. + index, w_index, bases_0, bases_1 = self._compute_index_and_bases(*eta, + deriv=True, rational=(weights is not None)) -# # Option 2: cycle over each element of 'coeffs' (touched only once) -# # - Pros: no temporary objects are created -# # - Cons: large number of Python iterations = number of elements in 'coeffs' -# # -# res = 0.0 -# for idx,c in np.ndenumerate( coeffs ): -# ndbasis = np.prod( [b[i] for i,b in zip( idx, bases )] ) -# res += c * ndbasis - return res + # If NURBS, evaluate weights field and its gradient + if weights: + w_coeffs = self._tmp_w_coeffs + w_coeffs[:] = weights[w_index] + w_inv_value = 1.0 / self.eval_field_once(w_coeffs, *bases_0, *self._work) + + # Evaluate each component of the gradient + w_grad = [] + for d in range(self.ldim): + bases = [(bases_1[d] if i==d else bases_0[i]) for i in range(self.ldim)] + w_res = self.eval_field_once(w_coeffs, *bases, *self._work) + w_grad.append(w_res) + + # Evaluate all fields at the same point + coeffs = self._tmp_coeffs + gradients_list = [] + for field in fields: + + # Get contiguous copy of the spline coefficients required for evaluation + coeffs[:] = field.coeffs._data[index] + + # If NURBS, rescale field coefficients by the weights' coefficients + if weights: + coeffs *= weights[w_index] + + # Evaluate the field + value = self.eval_field_once(coeffs, *bases_0, *self._work) + + # Evaluate each component of the gradient + grad = [] + for d in range(self.ldim): + bases = [(bases_1[d] if i==d else bases_0[i]) for i in range(self.ldim)] + grad_d = self.eval_field_once(coeffs, *bases, *self._work) + # If NURBS, apply chain rule of differentiation + if weights: + grad_d = (grad_d - w_grad[d] * value * w_inv_value) * w_inv_value + grad.append(grad_d) + + # Append field gradient (as a tuple) to list + gradients_list.append(tuple(grad)) + + return gradients_list # ... def preprocess_regular_tensor_grid(self, grid, der=0, overlap=0): @@ -612,65 +907,6 @@ def eval_fields_irregular_tensor_grid(self, grid, *fields, weights=None, overlap return out_fields - # ... - def eval_field_gradient(self, field, *eta, weights=None): - - assert isinstance(field, FemField) - assert field.space is self - assert len(eta) == self.ldim - - bases_0 = [] - bases_1 = [] - index = [] - - for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): - - knots = space.knots - degree = space.degree - span = find_span( knots, degree, x ) - #-------------------------------------------------# - # Fix span for boundaries between subdomains # - #-------------------------------------------------# - # TODO: Use local knot sequence instead of global # - # one to get correct span in all situations # - #-------------------------------------------------# - if x == xlim[1] and x != knots[-1-degree]: - span -= 1 - #-------------------------------------------------# - basis_0 = basis_funs(knots, degree, x, span) - basis_1 = basis_funs_1st_der(knots, degree, x, span) - - # If needed, rescale B-splines to get M-splines - if space.basis == 'M': - scaling = space.scaling_array[span-degree : span+1] - basis_0 *= scaling - basis_1 *= scaling - - # Determine local span - wrap_x = space.periodic and x > xlim[1] - loc_span = span - space.nbasis if wrap_x else span - - bases_0.append( basis_0 ) - bases_1.append( basis_1 ) - index.append( slice( loc_span-degree, loc_span+1 ) ) - - # Get contiguous copy of the spline coefficients required for evaluation - index = tuple( index ) - coeffs = field.coeffs[index].copy() - if weights: - coeffs *= weights[index] - - # Evaluate each component of the gradient using algorithm described in "Option 1" above - grad = [] - for d in range( self.ldim ): - bases = [(bases_1[d] if i==d else bases_0[i]) for i in range( self.ldim )] - res = coeffs - for basis in bases[::-1]: - res = np.dot( res, basis ) - grad.append( res ) - - return grad - # ... def integral(self, f, *, nquads=None): diff --git a/psydac/mapping/discrete.py b/psydac/mapping/discrete.py index b80532d94..22bfbd430 100644 --- a/psydac/mapping/discrete.py +++ b/psydac/mapping/discrete.py @@ -133,11 +133,11 @@ def from_control_points(cls, tensor_space, control_points): # Abstract interface #-------------------------------------------------------------------------- def __call__(self, *eta): - return [map_Xd(*eta) for map_Xd in self._fields] + return self.space.eval_fields_one_point(self._fields, *eta) # ... def jacobian(self, *eta): - return np.array([map_Xd.gradient(*eta) for map_Xd in self._fields]) + return np.array(self.space.eval_fields_gradient_one_point(self._fields, *eta)) # ... def jacobian_inv(self, *eta): @@ -873,20 +873,30 @@ def from_control_points_weights(cls, tensor_space, control_points, weights): #-------------------------------------------------------------------------- # Abstract interface #-------------------------------------------------------------------------- + # def __call__(self, *eta): + # map_W = self._weights_field + # w = map_W(*eta) + # Xd = [map_Xd(*eta , weights=map_W.coeffs) for map_Xd in self._fields] + # return np.asarray(Xd) / w + + # # ... + # def jacobian(self, *eta): + # map_W = self._weights_field + # w = map_W(*eta) + # grad_w = np.array(map_W.gradient(*eta)) + # v = np.array([map_Xd(*eta, weights=map_W.coeffs) for map_Xd in self._fields]) + # grad_v = np.array([map_Xd.gradient(*eta, weights=map_W.coeffs) for map_Xd in self._fields]) + # return grad_v / w - v[:, None] @ grad_w[None, :] / w**2 + def __call__(self, *eta): - map_W = self._weights_field - w = map_W(*eta) - Xd = [map_Xd(*eta , weights=map_W.coeffs) for map_Xd in self._fields] - return np.asarray(Xd) / w + return self.space.eval_fields_one_point(self._fields, *eta, + weights=self._weights_field.coeffs) # ... def jacobian(self, *eta): - map_W = self._weights_field - w = map_W(*eta) - grad_w = np.array(map_W.gradient(*eta)) - v = np.array([map_Xd(*eta, weights=map_W.coeffs) for map_Xd in self._fields]) - grad_v = np.array([map_Xd.gradient(*eta, weights=map_W.coeffs) for map_Xd in self._fields]) - return grad_v / w - v[:, None] @ grad_w[None, :] / w**2 + jac = self.space.eval_fields_gradient_one_point(self._fields, *eta, + weights=self._weights_field.coeffs) + return np.array(jac) #-------------------------------------------------------------------------- # Fast evaluation on a grid