I am having trouble extracting my results at certain points on my mesh. I understand that there are 2 methods to do this (i) at() function and (ii) VertexOnlyMesh. Both dont work for me. While the former raises
Point Not in Domain Error, the latter raises
TypeError: an integer is required. Here is a minimal example of trying both methods out:
from firedrake import *
from firedrake.petsc import PETSc
from firedrake.mesh import VertexOnlyMesh
from collections import OrderedDict
import numpy as np
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)
a = inner(grad(u), grad(v)) * dx
L = Constant(0) * v * dx
bc1 = DirichletBC(V, 0.0, 1)
bc2 = DirichletBC(V, 1.0, 2)
sol = Function(V)
problem1 = LinearVariationalProblem(a, L, sol, bcs=[bc1, bc2], constant_jacobian=True)
solver1 = LinearVariationalSolver(problem1, solver_parameters={"ksp_type": "cg", "pc_type": "jacobi"})
solver1.solve()
print(sol.dat.data_ro)
crd=[0.5256,0.5051]
num_points = 25
äs = np.linspace(0.0, 1.0, num_points + 1)
X, Y = np.meshgrid(äs, äs)
xs = np.vstack((X.flatten(), Y.flatten())).T
vm=1
pe=1
#Test VertexOnlyMesh Generation
if vm:
point_cloud = VertexOnlyMesh(mesh, xs)
#Test Point Evaluation of Function
if pe:
print(sol.at(crd))