diff --git a/psydac/feec/evaluate_dofs_k.py b/psydac/feec/evaluate_dofs_k.py new file mode 100644 index 000000000..4f5ed42b4 --- /dev/null +++ b/psydac/feec/evaluate_dofs_k.py @@ -0,0 +1,313 @@ +# +# Not yet pyccelizable since pyccel doesn't support Callable. +# See https://github.com/pyccel/pyccel/issues/1900 +# +from pyccel.decorators import template +from collections.abc import Callable + +#============================================================================== +# 1D DEGREES OF FREEDOM +#============================================================================== + +@template(name='T', types=[float, complex]) +def evaluate_dofs_1d_0form( + intp_x1:'T[:]', # interpolation points + F:'T[:]', # array of degrees of freedom (intent out) + f:Callable[['T[:]'], 'T'] # input scalar function (callable) + ): + + (n1,) = F.shape + + for i1 in range(n1): + F[i1] = f(intp_x1[i1]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_1d_1form( + quad_x1:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', # quadrature weights + F:'T[:]', # array of degrees of freedom (intent out) + f:Callable[['T[:,:]'], 'T'] # input scalar function (callable) + ): + + k1 = quad_x1.shape[1] + + n1, = F.shape + for i1 in range(n1): + F[i1] = 0.0 + for g1 in range(k1): + F[i1] += quad_w1[i1, g1] * f(quad_x1[i1, g1]) + +#============================================================================== +# 2D DEGREES OF FREEDOM +#============================================================================== + +@template(name='T', types=[float, complex]) +def evaluate_dofs_2d_0form( + intp_x1:'T[:]', intp_x2:'T[:]', # interpolation points + F:'T[:,:]', # array of degrees of freedom (intent out) + f:Callable[['T[:]', 'T[:]'], 'T'] # input scalar function (callable) + ): + + n1, n2 = F.shape + + for i1 in range(n1): + for i2 in range(n2): + F[i1, i2] = f(intp_x1[i1], intp_x2[i2]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_2d_1form_hcurl( + intp_x1:'T[:]' , intp_x2:'T[:]' , # interpolation points + quad_x1:'T[:,:]', quad_x2:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', # quadrature weights + F1:'T[:,:]', F2:'T[:,:]', # arrays of degrees of freedom (intent out) + f1:Callable[['T[:,:]', 'T[:]' ], 'T'], # input scalar functions (callable) + f2:Callable[['T[:]' , 'T[:,:]'], 'T'], # + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + + n1, n2 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + F1[i1, i2] = 0.0 + for g1 in range(k1): + F1[i1, i2] += quad_w1[i1, g1] * f1(quad_x1[i1, g1], intp_x2[i2]) + + n1, n2 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + F2[i1, i2] = 0.0 + for g2 in range(k2): + F2[i1, i2] += quad_w2[i2, g2] * f2(intp_x1[i1], quad_x2[i2, g2]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_2d_1form_hdiv( + intp_x1:'T[:]', intp_x2:'T[:]' , # interpolation points + quad_x1:'T[:,:]', quad_x2:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', # quadrature weights + F1:'T[:,:]', F2:'T[:,:]', # arrays of degrees of freedom (intent out) + f1:Callable[['T[:]' , 'T[:,:]'], 'T'], # input scalar functions (callable) + f2:Callable[['T[:,:]', 'T[:]' ], 'T'], # + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + + n1, n2 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + F1[i1, i2] = 0.0 + for g2 in range(k2): + F1[i1, i2] += quad_w2[i2, g2] * f1(intp_x1[i1], quad_x2[i2, g2]) + + n1, n2 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + F2[i1, i2] = 0.0 + for g1 in range(k1): + F2[i1, i2] += quad_w1[i1, g1] * f2(quad_x1[i1, g1], intp_x2[i2]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_2d_2form( + quad_x1:'T[:,:]', quad_x2:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', # quadrature weights + F:'T[:,:]', # array of degrees of freedom (intent out) + f:Callable[['T[:,:]', 'T[:,:]'], 'T'], # input scalar function (callable) + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + + n1, n2 = F.shape + for i1 in range(n1): + for i2 in range(n2): + F[i1, i2] = 0.0 + for g1 in range(k1): + for g2 in range(k2): + F[i1, i2] += quad_w1[i1, g1] * quad_w2[i2, g2] * \ + f(quad_x1[i1, g1], quad_x2[i2, g2]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_2d_vec( + intp_x1:'T[:]', intp_x2:'T[:]' , # interpolation points + F1:'T[:,:]', F2:'T[:,:]', # array of degrees of freedom (intent out) + f1:Callable[['T[:]', 'T[:]'], 'T'], # input scalar function (callable) + f2:Callable[['T[:]', 'T[:]'], 'T'], # + ): + + n1, n2 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + F1[i1, i2] = f1(intp_x1[i1], intp_x2[i2]) + + n1, n2 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + F2[i1, i2] = f2(intp_x1[i1], intp_x2[i2]) + + +#============================================================================== +# 3D DEGREES OF FREEDOM +#============================================================================== +@template(name='T', types=[float, complex]) +def evaluate_dofs_3d_0form( + intp_x1:'T[:]', intp_x2:'T[:]', intp_x3:'T[:]', # interpolation points + F:'T[:,:,:]', # array of degrees of freedom (intent out) + f:Callable[['T[:]','T[:]','T[:]'], 'T'] # input scalar function (callable) + ): + + n1, n2, n3 = F.shape + + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F[i1, i2, i3] = f(intp_x1[i1], intp_x2[i2], intp_x3[i3]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_3d_1form( + intp_x1:'T[:]', intp_x2:'T[:]', intp_x3:'T[:]', # interpolation points + quad_x1:'T[:,:]', quad_x2:'T[:,:]', quad_x3:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', quad_w3:'T[:,:]', # quadrature weights + F1:'T[:,:,:]', F2:'T[:,:,:]', F3:'T[:,:,:]', # arrays of degrees of freedom (intent out) + f1:Callable[['T[:,:]','T[:]' ,'T[:]' ], 'T'], # input scalar functions (callable) + f2:Callable[['T[:]' ,'T[:,:]','T[:]' ], 'T'], + f3:Callable[['T[:]' ,'T[:]' ,'T[:,:]'], 'T'] + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + k3 = quad_x3.shape[1] + + n1, n2, n3 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F1[i1, i2, i3] = 0.0 + for g1 in range(k1): + F1[i1, i2, i3] += quad_w1[i1, g1] * \ + f1(quad_x1[i1, g1], intp_x2[i2], intp_x3[i3]) + + n1, n2, n3 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F2[i1, i2, i3] = 0.0 + for g2 in range(k2): + F2[i1, i2, i3] += quad_w2[i2, g2] * \ + f2(intp_x1[i1], quad_x2[i2, g2], intp_x3[i3]) + + n1, n2, n3 = F3.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F3[i1, i2, i3] = 0.0 + for g3 in range(k3): + F3[i1, i2, i3] += quad_w3[i3, g3] * \ + f3(intp_x1[i1], intp_x2[i2], quad_x3[i3, g3]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_3d_2form( + intp_x1:'T[:]', intp_x2:'T[:]', intp_x3:'T[:]' , # interpolation points + quad_x1:'T[:,:]', quad_x2:'T[:,:]', quad_x3:'T[:,:]' , # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', quad_w3:'T[:,:]' , # quadrature weights + F1:'T[:,:,:]', F2:'T[:,:,:]', F3:'T[:,:,:]', # arrays of degrees of freedom (intent out) + f1:Callable[['T[:]' , 'T[:,:]', 'T[:,:]'], 'T'] , # input scalar functions (callable) + f2:Callable[['T[:,:]', 'T[:]', 'T[:,:]'], 'T'] , # + f3:Callable[['T[:,:]', 'T[:,:]', 'T[:]' ], 'T'] # + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + k3 = quad_x3.shape[1] + + n1, n2, n3 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F1[i1, i2, i3] = 0.0 + for g2 in range(k2): + for g3 in range(k3): + F1[i1, i2, i3] += quad_w2[i2, g2] * quad_w3[i3, g3] * \ + f1(intp_x1[i1], quad_x2[i2, g2], quad_x3[i3, g3]) + + n1, n2, n3 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F2[i1, i2, i3] = 0.0 + for g1 in range(k1): + for g3 in range(k3): + F2[i1, i2, i3] += quad_w1[i1, g1] * quad_w3[i3, g3] * \ + f2(quad_x1[i1, g1], intp_x2[i2], quad_x3[i3, g3]) + + n1, n2, n3 = F3.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F3[i1, i2, i3] = 0.0 + for g1 in range(k1): + for g2 in range(k2): + F3[i1, i2, i3] += quad_w1[i1, g1] * quad_w2[i2, g2] * \ + f3(quad_x1[i1, g1], quad_x2[i2, g2], intp_x3[i3]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_3d_3form( + quad_x1:'T[:,:]', quad_x2:'T[:,:]', quad_x3:'T[:,:]', # quadrature points + quad_w1:'T[:,:]', quad_w2:'T[:,:]', quad_w3:'T[:,:]', # quadrature weights + F:'T[:,:,:]', # array of degrees of freedom (intent out) + f:Callable[['T[:,:]','T[:,:]','T[:,:]'], 'T'], # input scalar function (callable) + ): + + k1 = quad_x1.shape[1] + k2 = quad_x2.shape[1] + k3 = quad_x3.shape[1] + + n1, n2, n3 = F.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F[i1, i2, i3] = 0.0 + for g1 in range(k1): + for g2 in range(k2): + for g3 in range(k3): + F[i1, i2, i3] += \ + quad_w1[i1, g1] * quad_w2[i2, g2] * quad_w3[i3, g3] * \ + f(quad_x1[i1, g1], quad_x2[i2, g2], quad_x3[i3, g3]) + +#------------------------------------------------------------------------------ +@template(name='T', types=[float, complex]) +def evaluate_dofs_3d_vec( + intp_x1:'T[:]', intp_x2:'T[:]', intp_x3:'T[:]', # interpolation points + F1:'T[:,:,:]', F2:'T[:,:,:]', F3:'T[:,:,:]', # arrays of degrees of freedom (intent out) + f1:Callable[['T[:]', 'T[:]', 'T[:]'], 'T'] , # input scalar functions (callable) + f2:Callable[['T[:]', 'T[:]', 'T[:]'], 'T'] , # + f3:Callable[['T[:]', 'T[:]', 'T[:]'], 'T'] # + ): + + # evaluate input functions at interpolation points (make sure that points are in [0, 1]) + n1, n2, n3 = F1.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F1[i1, i2, i3] = f1(intp_x1[i1], intp_x2[i2], intp_x3[i3]) + + n1, n2, n3 = F2.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F2[i1, i2, i3] = f2(intp_x1[i1], intp_x2[i2], intp_x3[i3]) + + n1, n2, n3 = F3.shape + for i1 in range(n1): + for i2 in range(n2): + for i3 in range(n3): + F3[i1, i2, i3] = f3(intp_x1[i1], intp_x2[i2], intp_x3[i3]) \ No newline at end of file diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index bdccc0e87..3792fe9f6 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -14,12 +14,16 @@ from psydac.utilities.utils import roll_edges +from psydac.feec.evaluate_dofs_k import (evaluate_dofs_1d_0form, evaluate_dofs_1d_1form, + evaluate_dofs_2d_0form, evaluate_dofs_2d_1form_hcurl, evaluate_dofs_2d_1form_hdiv, evaluate_dofs_2d_2form, + evaluate_dofs_3d_0form, evaluate_dofs_3d_1form, evaluate_dofs_3d_2form, evaluate_dofs_3d_3form, + evaluate_dofs_2d_vec, evaluate_dofs_3d_vec) + from abc import ABCMeta, abstractmethod -__all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2', - 'evaluate_dofs_1d_0form', 'evaluate_dofs_1d_1form', - 'evaluate_dofs_2d_0form', 'evaluate_dofs_2d_1form_hcurl', 'evaluate_dofs_2d_1form_hdiv', 'evaluate_dofs_2d_2form', - 'evaluate_dofs_3d_0form', 'evaluate_dofs_3d_1form', 'evaluate_dofs_3d_2form', 'evaluate_dofs_3d_3form') + + +__all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2') #============================================================================== class GlobalProjector(metaclass=ABCMeta): @@ -636,277 +640,4 @@ def __call__(self, fun): finite element space). This is also a real-valued vector function in the logical domain. """ - return super().__call__(fun) - -#============================================================================== -# 1D DEGREES OF FREEDOM -#============================================================================== - -def evaluate_dofs_1d_0form(intp_x1, F, f): - (n1,) = F.shape - - for i1 in range(n1): - F[i1] = f(intp_x1[i1]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_1d_1form( - quad_x1, # quadrature points - quad_w1, # quadrature weights - F, # array of degrees of freedom (intent out) - f # input scalar function (callable) - ): - - k1 = quad_x1.shape[1] - - n1, = F.shape - for i1 in range(n1): - F[i1] = 0.0 - for g1 in range(k1): - F[i1] += quad_w1[i1, g1] * f(quad_x1[i1, g1]) - -#============================================================================== -# 2D DEGREES OF FREEDOM -#============================================================================== - -def evaluate_dofs_2d_0form(intp_x1, intp_x2, F, f): - n1, n2 = F.shape - - for i1 in range(n1): - for i2 in range(n2): - F[i1, i2] = f(intp_x1[i1], intp_x2[i2]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_2d_1form_hcurl( - intp_x1, intp_x2, # interpolation points - quad_x1, quad_x2, # quadrature points - quad_w1, quad_w2, # quadrature weights - F1, F2, # arrays of degrees of freedom (intent out) - f1, f2 # input scalar functions (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - - n1, n2 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - F1[i1, i2] = 0.0 - for g1 in range(k1): - F1[i1, i2] += quad_w1[i1, g1] * f1(quad_x1[i1, g1], intp_x2[i2]) - - n1, n2 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - F2[i1, i2] = 0.0 - for g2 in range(k2): - F2[i1, i2] += quad_w2[i2, g2] * f2(intp_x1[i1], quad_x2[i2, g2]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_2d_1form_hdiv( - intp_x1, intp_x2, # interpolation points - quad_x1, quad_x2, # quadrature points - quad_w1, quad_w2, # quadrature weights - F1, F2, # arrays of degrees of freedom (intent out) - f1, f2 # input scalar functions (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - - n1, n2 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - F1[i1, i2] = 0.0 - for g2 in range(k2): - F1[i1, i2] += quad_w2[i2, g2] * f1(intp_x1[i1], quad_x2[i2, g2]) - - n1, n2 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - F2[i1, i2] = 0.0 - for g1 in range(k1): - F2[i1, i2] += quad_w1[i1, g1] * f2(quad_x1[i1, g1], intp_x2[i2]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_2d_2form( - quad_x1, quad_x2, # quadrature points - quad_w1, quad_w2, # quadrature weights - F, # array of degrees of freedom (intent out) - f, # input scalar function (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - - n1, n2 = F.shape - for i1 in range(n1): - for i2 in range(n2): - F[i1, i2] = 0.0 - for g1 in range(k1): - for g2 in range(k2): - F[i1, i2] += quad_w1[i1, g1] * quad_w2[i2, g2] * \ - f(quad_x1[i1, g1], quad_x2[i2, g2]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_2d_vec( - intp_x1, intp_x2, # interpolation points - F1, F2, # array of degrees of freedom (intent out) - f1, f2, # input scalar function (callable) - ): - - n1, n2 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - F1[i1, i2] = f1(intp_x1[i1], intp_x2[i2]) - - n1, n2 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - F2[i1, i2] = f2(intp_x1[i1], intp_x2[i2]) - - -#============================================================================== -# 3D DEGREES OF FREEDOM -#============================================================================== - -def evaluate_dofs_3d_0form(intp_x1, intp_x2, intp_x3, F, f): - n1, n2, n3 = F.shape - - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F[i1, i2, i3] = f(intp_x1[i1], intp_x2[i2], intp_x3[i3]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_3d_1form( - intp_x1, intp_x2, intp_x3, # interpolation points - quad_x1, quad_x2, quad_x3, # quadrature points - quad_w1, quad_w2, quad_w3, # quadrature weights - F1, F2, F3, # arrays of degrees of freedom (intent out) - f1, f2, f3 # input scalar functions (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - k3 = quad_x3.shape[1] - - n1, n2, n3 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F1[i1, i2, i3] = 0.0 - for g1 in range(k1): - F1[i1, i2, i3] += quad_w1[i1, g1] * \ - f1(quad_x1[i1, g1], intp_x2[i2], intp_x3[i3]) - - n1, n2, n3 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F2[i1, i2, i3] = 0.0 - for g2 in range(k2): - F2[i1, i2, i3] += quad_w2[i2, g2] * \ - f2(intp_x1[i1], quad_x2[i2, g2], intp_x3[i3]) - - n1, n2, n3 = F3.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F3[i1, i2, i3] = 0.0 - for g3 in range(k3): - F3[i1, i2, i3] += quad_w3[i3, g3] * \ - f3(intp_x1[i1], intp_x2[i2], quad_x3[i3, g3]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_3d_2form( - intp_x1, intp_x2, intp_x3, # interpolation points - quad_x1, quad_x2, quad_x3, # quadrature points - quad_w1, quad_w2, quad_w3, # quadrature weights - F1, F2, F3, # arrays of degrees of freedom (intent out) - f1, f2, f3 # input scalar functions (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - k3 = quad_x3.shape[1] - - n1, n2, n3 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F1[i1, i2, i3] = 0.0 - for g2 in range(k2): - for g3 in range(k3): - F1[i1, i2, i3] += quad_w2[i2, g2] * quad_w3[i3, g3] * \ - f1(intp_x1[i1], quad_x2[i2, g2], quad_x3[i3, g3]) - - n1, n2, n3 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F2[i1, i2, i3] = 0.0 - for g1 in range(k1): - for g3 in range(k3): - F2[i1, i2, i3] += quad_w1[i1, g1] * quad_w3[i3, g3] * \ - f2(quad_x1[i1, g1], intp_x2[i2], quad_x3[i3, g3]) - - n1, n2, n3 = F3.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F3[i1, i2, i3] = 0.0 - for g1 in range(k1): - for g2 in range(k2): - F3[i1, i2, i3] += quad_w1[i1, g1] * quad_w2[i2, g2] * \ - f3(quad_x1[i1, g1], quad_x2[i2, g2], intp_x3[i3]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_3d_3form( - quad_x1, quad_x2, quad_x3, # quadrature points - quad_w1, quad_w2, quad_w3, # quadrature weights - F, # array of degrees of freedom (intent out) - f, # input scalar function (callable) - ): - - k1 = quad_x1.shape[1] - k2 = quad_x2.shape[1] - k3 = quad_x3.shape[1] - - n1, n2, n3 = F.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F[i1, i2, i3] = 0.0 - for g1 in range(k1): - for g2 in range(k2): - for g3 in range(k3): - F[i1, i2, i3] += \ - quad_w1[i1, g1] * quad_w2[i2, g2] * quad_w3[i3, g3] * \ - f(quad_x1[i1, g1], quad_x2[i2, g2], quad_x3[i3, g3]) - -#------------------------------------------------------------------------------ -def evaluate_dofs_3d_vec( - intp_x1, intp_x2, intp_x3, # interpolation points - F1, F2, F3, # array of degrees of freedom (intent out) - f1, f2, f3, # input scalar function (callable) - ): - - # evaluate input functions at interpolation points (make sure that points are in [0, 1]) - n1, n2, n3 = F1.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F1[i1, i2, i3] = f1(intp_x1[i1], intp_x2[i2], intp_x3[i3]) - - n1, n2, n3 = F2.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F2[i1, i2, i3] = f2(intp_x1[i1], intp_x2[i2], intp_x3[i3]) - - n1, n2, n3 = F3.shape - for i1 in range(n1): - for i2 in range(n2): - for i3 in range(n3): - F3[i1, i2, i3] = f3(intp_x1[i1], intp_x2[i2], intp_x3[i3]) - + return super().__call__(fun) \ No newline at end of file