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


61
views
0
11 days ago by
Leo  
I have a time dependent problem. I want to print the solution at a certain coordinate in each time step, Here is the code:

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
10 days ago by
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:

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]​
Thanks but it gives me this Syntax error:

SyntaxError: invalid character in identifier
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?
written 10 days ago by Leo  
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.

Similar posts:
Search »