Probing field returns value error

[fanoway]

Tryin to probe a solved problem return the following error.

Traceback (most recent call last):
  File "C:\Users\user name\Downloads\sfep_model\app.py", line 202, in <module>
    model()
  File "C:\Users\User Name\Downloads\sfep_model\app.py", line 188, in model
    fig, results = probe_results(ii, t, probe, labels[ii])
  File "C:\Users\User Name\Downloads\sfep_model\app.py", line 98, in probe_results
    pars, vals = probe(t)
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\probes.py", line 271, in __call__
    return self.probe(variable, **kwargs)
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\probes.py", line 313, in probe
    vals, ref_coors, cells, status = ev(
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\variables.py", line 2022, in evaluate_at
    out = self.field.evaluate_at(coors, source_vals,
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\common\fields.py", line 413, in evaluate_at
    ref_coors, cells, status = get_ref_coors(self, coors,
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\common\global_interp.py", line 350, in get_ref_coors
    return get_ref_coors_general(field, coors, close_limit=close_limit,
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\common\global_interp.py", line 266, in get_ref_coors_general
    potential_cells, offsets = get(coors, cmesh, centroids=centroids,
  File "C:\Users\User Name\miniconda3\envs\heat_model\lib\site-packages\sfepy\discrete\common\global_interp.py", line 163, in get_potential_cells
    ips = kdtree.query_ball_point(centroid, radii[ic], p=nm.inf)
  File "_ckdtree.pyx", line 940, in scipy.spatial._ckdtree.cKDTree.query_ball_point
  File "_ckdtree.pyx", line 1619, in scipy.spatial._ckdtree.num_points
ValueError: x must consist of vectors of length 3 but has shape (2,)

Here is my source code

import meshio
import numpy as np
import pygmsh
from matplotlib import pyplot as plt
from sfepy.base.base import IndexedStruct
from sfepy.discrete import (
    Equation,
    Equations,
    FieldVariable,
    Integral,
    Material,
    Problem,
)
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.discrete.fem import FEDomain, Field, Mesh
from sfepy.discrete.probes import LineProbe
from sfepy.postprocess.viewer import Viewer
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.terms import Term

def create_mesh():
    with pygmsh.geo.Geometry() as geom:
        geom.add_polygon(
            [
                [0.0, 0.0],
                [0, 1],
                [1, 1],
                [1.2, 0.7],
                [1.2, 0],
            ],
            mesh_size=0.1,
        )
        mesh = geom.generate_mesh()

    return mesh

def gen_probes(problem, x_range, y_range):
    """
    Define a line probe probe.
    """

    # Use enough points for higher order approximations.
    n_point = 1000

    p0, p1 = np.array([x_range[0], 0.5 * sum(y_range), 0.0]), np.array(
        [x_range[1], 0.5 * sum(y_range), 0.0]
    )
    line = LineProbe(p0, p1, n_point, share_geometry=True)
    # Workaround current probe code shortcoming.
    line.set_options(close_limit=0.5)

    probes = [line]

    pars, points = line.get_points()

    print(f"{points = }")
    labels = ["%s -> %s" % (p0, p1)]

    return probes, labels

def probe_results(ax_num, t, probe, label):
    """
    Probe the results using the given probe and plot the probed values.
    """

    results = {}

    pars, vals = probe(t)
    results["t"] = (pars, vals)

    fig = plt.figure(1)

    ax = plt.subplot(2, 2, 2 * ax_num + 1)
    ax.cla()
    pars, vals = results["t"]
    ax.plot(pars, vals, label=r"$T$", lw=1, ls="-", marker="+", ms=3)
    dx = 0.05 * (pars[-1] - pars[0])
    ax.set_xlim(pars[0] - dx, pars[-1] + dx)
    ax.set_ylabel("temperature")
    ax.set_xlabel("probe %s" % label, fontsize=8)
    ax.legend(loc="best", fontsize=10)

    return fig, results

def model():
    mesh = Mesh.from_file("working_mesh.vtk")

    copper_diffusion = 111  # mm**2/s

    domain = FEDomain("domain", mesh)
    x_range = domain.get_mesh_bounding_box()[:, 0]
    y_range = domain.get_mesh_bounding_box()[:, 1]

    eps = 1e-8 * (x_range[1] - x_range[0])

    omega = domain.create_region("Omega", "all")

    Left = domain.create_region(
        "Gamma1", "vertices in x < %.10f" % (x_range[0] + eps), "facet"
    )
    Right = domain.create_region(
        "Gamma2", "vertices in x > %.10f" % (x_range[1] - eps), "facet"
    )

    field = Field.from_args("temperature", np.float64, 1, omega, approx_order=2)

    t = FieldVariable("t", "unknown", field)
    s = FieldVariable("s", "test", field, primary_var_name="t")

    cu = Material("cu", diffusivity=copper_diffusion * np.eye(2))

    integral = Integral("i", order=2)

    diffusion = Term.new(
        "dw_diffusion(cu.diffusivity, s, t)", integral, omega, cu=cu, s=s, t=t
    )

    eq = Equation("thermal", diffusion)
    eqs = Equations([eq])

    temp1 = EssentialBC("temp1", Left, {"t.0": 2.0})
    temp2 = EssentialBC("temp2", Right, {"t.0": -2.0})

    ls = ScipyDirect({})
    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem("heat", equations=eqs)

    pb.save_regions_as_groups("regions")
    view = Viewer("regions.vtk")
    view(is_wireframe=True)

    pb.set_bcs(ebcs=Conditions([temp1, temp2]))
    pb.set_solver(nls)

    status = IndexedStruct()
    variables = pb.solve(
        status=status,
    )

    print("\n######")
    print("Nonlinear solver status:\n", nls_status)
    print("\n######")
    print("Stationary solver status:\n", status)
    print("######")

    out = variables.create_output()

    pb.save_state("thermal.vtk", out=out)
    view = Viewer("thermal.vtk")
    view(is_wireframe=True, is_scalar_bar=True)

    probes, labels = gen_probes(pb, x_range, y_range)

    for ii, probe in enumerate(probes):
        fig, results = probe_results(ii, t, probe, labels[ii])

if __name__ == "__main__":

    mesh = create_mesh()

    meshio.write(
        "working_mesh.vtk",
        meshio.Mesh(points=mesh.points, cells={"triangle": mesh.cells[1].data}),
    )

    model()

[fanoway]

I've done some digging and it seems the line 184 in probes.py is returning an array of length 2 vectors, shape = (n,2) while the KDTree method is expecting 3 dimensional vectors.

[fanoway]

I've discovered that this is because my mesh was planar (2D) and my line probe was specified in 3 dimensions.