kinnala / scikit-fem

Simple finite element assemblers
https://scikit-fem.readthedocs.io
BSD 3-Clause "New" or "Revised" License
500 stars 79 forks source link

Interpolator with ElementVector #973

Open HelgeGehring opened 1 year ago

HelgeGehring commented 1 year ago
import numpy as np
from skfem import *

mesh = MeshTri()
basis = Basis(mesh, ElementTriP1())
basis.interpolator(basis.zeros())(np.array([[[0, 0]]]).swapaxes(0, -1))

works, but

import numpy as np
from skfem import *

mesh = MeshTri()
basis = Basis(mesh, ElementVector(ElementTriP1()))
basis.interpolator(basis.zeros())(np.array([[[0, 0]]]).swapaxes(0, -1))

doesn't.

ValueError: row, column, and data array must all be the same length
kinnala commented 1 year ago

This should be fixed first:

basis.probes(np.array([[0], [0]]))

which leads to the same error.

kinnala commented 1 year ago

Obviously fixing this requires detection of the tensorial order of the field. This is kind of already done in Basis.refinterp by passing test values to Element.gbasis and checking the shape of the output. This should be moved to Basis._detect_tensor_order and used both in Basis.refinterp and Basis.probes.

Nevertheless, fixing this requires some thorough thinking how Basis.probes should be generalized. Since Basis.probes returns a matrix and the different elements in the range of the linear operator described by the matrix are already describing the different points passed to Basis.probes, perhaps Basis.probes needs to return multiple matrices instead?

kinnala commented 1 year ago

Another option is to use Basis._detect_tensor_order already in Basis.interpolator and skip the reuse of Basis.probes if the output of Basis._detect_tensor_order is higher than scalar field, and rely on a different implementation.