pyccel / psydac

Python 3 library for isogeometric analysis
https://pyccel.github.io/psydac/
MIT License
52 stars 18 forks source link

evaluate_fields method not working for FEM spaces in 1D derham sequence #382

Closed wbarham closed 7 months ago

wbarham commented 8 months ago

The evaluate_fields method may not be implemented for the FEM spaces in the 1D derham sequence.

As a minimal example:

import numpy as np

from sympde.topology           import Line
from sympde.topology           import Derham
from mpi4py                    import MPI
from psydac.api.discretization import discretize

ncells   = 10
degree   = 2
periodic = False

domain   = Line('Omega', bounds=(0, 1))
derham   = Derham(domain)

domain_h = discretize(domain, ncells=[ncells], periodic=[periodic])
derham_h = discretize(derham, domain_h, degree=[degree])

P0, P1 = derham_h.projectors(nquads=[degree+2])

f0 = P0(lambda x : np.cos(2 * np.pi * x) )
f1 = P1(lambda x : np.cos(2 * np.pi * x) )

eval_grid = np.linspace(0, 1, num=51, endpoint=True)

derham_h.V0.eval_fields((eval_grid,), f0)
derham_h.V1.eval_fields((eval_grid,), f1)

I get the following error message:

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[28], line 25
     21 f1 = P1(lambda x : np.cos(2 * np.pi * x) )
     23 eval_grid = np.linspace(0, 1, num=51, endpoint=True)
---> 25 derham_h.V0.eval_fields((eval_grid,), f0)
     26 derham_h.V1.eval_fields((eval_grid,), f1)

File ~/virtualenvs/pyccel_env/lib/python3.10/site-packages/psydac/fem/tensor.py:446, in TensorFemSpace.eval_fields(self, grid, weights, npts_per_cell, overlap, *fields)
    443 # Case 2. 1D array of coordinates and no npts_per_cell is given
    444 # -> grid is tensor-product, but npts_per_cell is not the same in each cell
    445 elif grid[0].ndim == 1 and npts_per_cell is None:
--> 446     out_fields = self.eval_fields_irregular_tensor_grid(grid, *fields, weights=weights, overlap=overlap)
    447     return [np.ascontiguousarray(out_fields[..., i]) for i in range(len(fields))]
    449 # Case 3. 1D arrays of coordinates and npts_per_cell is a tuple or an integer
    450 # -> grid is tensor-product, and each cell has the same number of evaluation points

File ~/virtualenvs/pyccel_env/lib/python3.10/site-packages/psydac/fem/tensor.py:606, in TensorFemSpace.eval_fields_irregular_tensor_grid(self, grid, weights, overlap, *fields)
    603         eval_fields_3d_irregular_weighted(*npoints, *degree, *cell_indexes, *global_basis,
    604                                           *global_spans, glob_arr_coeffs, global_weight_coeff, out_fields)
    605 else:
--> 606     raise NotImplementedError("1D not Implemented")
    608 return out_fields

NotImplementedError: 1D not Implemented
yguclu commented 7 months ago

Fixed by #394