### Bug in printing the solution at a coordinate in case of quad mesh

61

views

0

I have a time dependent problem. I want to print the solution at a certain coordinate in each time step, Here is the code:

I have tried a lot of things like using:

```
from dolfin import *
import numpy as np
from numpy import cos , pi
import mshr
D = 7. * (10**(-11.))
R = 8.31
C_0 = Constant(500.0)
epsilon = 0.025
Temp = 293.0
mu = D / (R*Temp)
F_N = 96485.34
Charge_num=Constant(1.)
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.01)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], (200.0 * (10**(-6))))
left = Left()
top = Top()
right = Right()
bottom = Bottom()
def mesh2(Nx,lx,ly):
m = UnitSquareMesh.create(5*Nx, Nx, CellType.Type_quadrilateral)
x = m.coordinates()
x[:,:] = (x[:,:] - 0.5) * 2.
x[:,:] = 0.5 * (cos(pi * (x[:,:] - 1.) / 2.) + 1.)
x[:,0] = x[:,0]*lx
x[:,1] = x[:,1]*ly
return m
mesh = mesh2(30,0.01,200.*(10**-6))
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 2)
V_c = FunctionSpace(mesh, Element1)
V_phi = FunctionSpace(mesh, Element1)
W_elem = MixedElement([Element1, Element1])
W = FunctionSpace(mesh, W_elem)
z = Function(W)
dz=TrialFunction(W)
c, phi = split(z)
(v_1, v_2) = TestFunctions(W)
dt = 0.04
t = 0
T = 1
W_const = C_0
C_previous = interpolate(W_const, V_c)
Voltage = Constant(0.5)
bc_top = DirichletBC(W.sub(1), Voltage, boundaries, 1)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)
bcs = [bc_top, bc_bottom]
F = dot(grad(phi), grad(v_2)) * dx \
+ c * v_1 * dx \
- (F_N / epsilon) * c * v_2 * dx \
+ dt * D * dot(grad(c), grad(v_1)) * dx \
+ dt * Charge_num * mu * F_N * c * dot(grad(phi), grad(v_1)) * dx \
- C_previous * v_1 * dx + (F_N / epsilon) * C_0 * v_2 * dx
while t <= T:
J = derivative(F, z, dz)
problem = NonlinearVariationalProblem(F, z, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
(c, phi) = z.split(True)
C_previous.assign(c)
t += dt
#Here is where the problem is!!!
print (c(0.005, 0.00))
```

This code works fine and I can visualize the results in Paraview if I remove the last line (**print (c(0.005, 0.00))**). I added this line because I want to print the value of the variable**c**in every time step at coordinate (0.005 , 0.0) but it gives me this error:```
Error: Unable to intersect cell and point.
Reason: Intersection is only implemented for simplex meshes.
```

I have tried a lot of things like using:

**print (z.sub(0)(0.005, 0.0))**but none of them worked. Any idea how to resolve this issue?
Community: FEniCS Project

### 1 Answer

1

based on https://fenicsproject.org/qa/6450/get-field-value-at-mesh-vertex/ your approach should work.. apparently it's not working on the current version..

A workaround could be:

A workaround could be:

```
for index, v in enumerate(mesh.coordinates()):
if v[0] == 0.005 and v[1] == 0.0:
print project(c,V_c).vector().get_local()[index]
```

I think that point evaluation is not well supported for quads and hex meshes. Try with triangular cells.

written
10 days ago by
bleyerj

1

I guess you are right. I tried the simple Poisson demo :

```
from dolfin import *
# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
# Compute solution
u = Function(V)
solve(a == L, u, bc)
print (u(0.5,0.5))
print (u(0.8,0.8))
file = File("poisson.pvd")
file << u
print (u(0.5,0.5))
```

It works fine. Here is the output:```
Solving linear variational problem.
0.252805505996
```

If I change the mesh to quadrilateral by replacing the mesh with:`mesh = UnitSquareMesh.create( 32, 32, CellType.Type_quadrilateral)`

I get this again:```
*** Error: Unable to intersect cell and point.
*** Reason: Intersection is only implemented for simplex meshes.
```

But the question is : How could we do it when we are dealing with quad elements? It was officially announced that the FEniCS supports quad elements but even a simple evaluation of results at a coordinate is not possible! Isn't there really no way to do it?!
written
10 days ago by
Leo

Please login to add an answer/comment or follow this question.

I remember the :

print (c(0.005, 0.00))used to work perfectly in the previous versions of DOLFIN. The version of DOLFIN that I am using is :2017.2.0

I have provided the complete code above. Any idea about how to extract results at an arbitrary point? What could be wrong with that? Is it a BUG or something like that?