Why am i getting the "The point is not inside the domain." error when trying to evaluate solution at the point inside the domain of an imported mesh

General
1

I have constructed a dummy model to illustrate the problem. In the model I am evaluating the field inside a parallel plate capacitor composed of 1m^2 parallel plates located at x=+/-0.5m held at +/-1V. I first used COMSOL to construct a mesh of a 1m^3 cube that I exported to an .stl file. Using the command,

meshio convert Block_Mesh.stl Block_Mesh.xdmf,

I converted the mesh to the .xdmf file format compatible with FeNics, and successfully imported the mesh into my python code, as seen plotted below.
Block_Mesh

My code is as follows,

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

def Boundary(r):
    tolerance = 1E-12
    if np.abs(r[0]-0.5)<tolerance or np.abs(r[0]+0.5)<tolerance:
        return 1
    else:
        return 0

def BC(x):    
    if x<0:
        return -1
    if x>0:
        return 1

mesh = Mesh()
with XDMFFile("Block_Mesh.xdmf") as infile:
    infile.read(mesh)
ax = plot(mesh)
plt.savefig("Block_Mesh.png")

V = FunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)

class MyExpression0(UserExpression):
    def eval(self, value, x):
        value[0] = BC(x[0])

u_d = MyExpression0()

bc = DirichletBC(V, u_d, Boundary)

f = Constant(0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

u = Function(V)

solve(a == L, u, bc)

c=plot(u)
plt.colorbar(c)
plt.savefig("soln.png")
print(u(Point(0,0,0))) #point is clearly in domain

I think this obtains the correct solution, because the plot
soln this is what you’d expect if you plotted the solution for a slice in the x-y plane. Unforturnately, because im running ubuntu via WSL i cannot access the interactive feature of the plot, so I have just been assuming that this is what the plot is doing.
The problem is, even though (0,0,0) is clearly in the domain, the final line where i try to evaluate the solution at this point raises the error
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling “Function::set_allow_extrapolation(true)” on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
Why is this the case? I tried something similar in 2D with an imported mesh and didn’t have this problem. Is there a way to prevent this in 3D?

2

Without having access to the mesh it’s difficult to give advice.

Are you sure you have a volume mesh? I thought stereolithography define surfaces only. If so, the point (0, 0, 0) is clearly outside of your domain.

What if you replace your mesh with the inbuilt BoxMesh of appropriate size as a test?

1 Like
3

My bad, I forgot to upload them. Here are the mesh files exported from comsol in the two available file formats that are compatible with meshio, .stl: Dropbox - Block_Mesh.stl - Simplify your life and .nas Dropbox - Block_Mesh.nas - Simplify your life.

I thought stereolithography define surfaces only. If so, the point (0,0,0)
is clearly outside of your domain.

I’m not hugely knowledgable when it comes to these file formats, so that could explain it. From what I’ve read, nastran files can store volume meshes, however after conversion to .xdmf the python script gives the error
*** Error: Unable to recognise cell type.
*** Reason: Unknown value “tetrahedron_10”.
when trying to read the mesh.

4

This is because the mesh is second order. Legacy Dolfin does not support second order meshes. You would have to use DOLFINx for that functionality.

2 Likes
5

ah alright thanks, so if I download dolfinx I’ll be able to read the converted .nas?

6

You can use meshio to write the second order mesh to xdmf, and read it in dolfinx.

1 Like