Hi
I want to evaluate a variable at a particular point of a domain. My domain is a rectangle. In a part of my code where I have defined my mesh (inside a class) I have moved the mesh coordinates upward (in the y-direction):
def mesh(self):
.
.
m = RectangleMesh(Point(0,0), Point(1, 1),self.Nx, self.Ny)
x = m.coordinates()
.
.
for x in m.coordinates(): #moving mesh in the y-direction
x[1] += 0.000002
return m
I want to evaluate the variable "c" (in my problem) at the upper edge (left corner).
print c(0.0,202.0 * (10**(-6)))
I get this error:
*** Reason: The point is not inside the domain...
While the y=202.0 * (10**(-6)) is obviously inside the domain!! (It is the upper edge of the domain)!
In addition when I do:
print c(0.0,201.9999 * (10**(-6)))
I do not get the error anymore and it gives me the value of "c" at this point but I need to get this value at the coordinate:
(0.0,202.0 * (10**(-6)))
Here is the complete code:
from dolfin import *
from numpy import cos, pi
import numpy as np
D = 1. * (10**(-9))
R = 8.31
C_0 = 1.
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.03)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.000002)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], (202. * (10**(-6))))
left = Left()
top = Top()
right = Right()
bottom = Bottom()
class MyMesh:
def __init__(self, lx,ly, Nx,Ny):
self.lx = lx
self.ly = ly
self.Nx = Nx
self.Ny = Ny
def mesh(self):
m = RectangleMesh(Point(0,0), Point(1, 1),self.Nx, self.Ny)
x = m.coordinates()
x[:,:] = (x[:,:] - 0.5) * 2.
x[:,:] = 0.5 * (cos(pi * (x[:,:] - 1.) / 2.) + 1.)
x[:,0] = x[:,0]*self.lx
x[:,1] = x[:,1]*self.ly
for x in m.coordinates():
x[1] += 0.000002
return m
mesh0=MyMesh(0.03,200.*(10 ** (-6)), 500,30)
mesh1 = mesh0.mesh()
boundaries = FacetFunction("size_t", mesh1)
boundaries.set_all(0)
left.mark(boundaries, 2)
top.mark(boundaries, 1)
right.mark(boundaries, 4)
bottom.mark(boundaries, 3)
CG1_elem = FiniteElement("CG", mesh1.ufl_cell(), 1)
V_c = FunctionSpace(mesh1, CG1_elem)
V_phi = FunctionSpace(mesh1, CG1_elem)
W_elem = MixedElement([CG1_elem, CG1_elem])
W = FunctionSpace(mesh1, W_elem)
z = Function(W)
dz=TrialFunction(W)
c, phi = split(z)
(v_1, v_2) = TestFunctions(W)
dt = 0.01
t = 0
T = 0.1
# Previous solution
W_const = Constant(1.0)
C_previous = interpolate(W_const, V_c)
bc_top = DirichletBC(W.sub(1), Constant(2.0), boundaries, 1)
bc_bottom = DirichletBC(W.sub(1), 0.0, boundaries, 3)
bcs = [bc_top, bc_bottom]
# Variational form For nonlinear solver
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
plot(c,key="c", interactive=False)
print c(0.0,202.0 * (10**(-6)))
Any idea how to fix this issue?
Thanks!