New question #152507 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/152507
I'm working on a second-order elliptic equation, but I have some strange
effects regarding my solution, so I went back to the demo on solution of
poisson's equation with neumann bc, and the integral over the domain equal to
zero for uniqueness property.
This demo works fine of course, but when I incorporate some SubDomains, and
split the unitcircle in three different domains, there is trouble. Just to
control that I have done it right, I do use all the same coefficients in all
subdomains, and when I do solve this modified version, I get the same solution.
But, if I call the same code several times, without actually leaving python for
each call, I get different results each time.
My guess is that some information, which affects the solver, gets stored from
one call to the next. Does anyone have any thoughts on what might be wrong?
Here is the code:
from dolfin import *
# Create mesh and define function space
mesh = UnitCircle(50)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
#Separate the mesh into different regions using the SubDomain class
class dom1(SubDomain):
def inside(self, x, on_boundary):
return True if x[0]*x[0]+x[1]*x[1] <= 0.2*0.2 \
else False
class dom2(SubDomain):
def inside(self, x, on_boundary):
return True if (x[0]*x[0]+x[1]*x[1] >= 0.2*0.2 and \
x[0]*x[0]+x[1]*x[1] <= 0.4*0.4) else False
class dom3(SubDomain):
def inside(self, x, on_boundary):
return True if x[0]*x[0]+x[1]*x[1] >= 0.4*0.4 \
else False
#Use the MeshFunction to assign a number to each cell,
#telling them which part of the domain they're in
subdomains = MeshFunction('uint', mesh, 2)
# Mark subdomains with numbers 0, 1 and 2
subdomain0 = dom1()
subdomain0.mark(subdomains, 0)
subdomain1 = dom2()
subdomain1.mark(subdomains, 1)
subdomain2 = dom3()
subdomain2.mark(subdomains, 2)
A1 = as_matrix ([[1.0, 0.0], [0.0, 1.0]])
# Define variational problem
(v, d) = TestFunctions(W)
(u, c) = TrialFunction(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = (inner(grad(v), A1*grad(u)) + v*c + d*u)*dx(0) + \
(inner(grad(v), A1*grad(u)) + v*c + d*u)*dx(1) + \
(inner(grad(v), A1*grad(u)) + v*c + d*u)*dx(2)
L = v*f*dx + v*g*ds
# Compute solution
A = assemble(a, cell_domains = subdomains)
b = assemble(L)
U = Function(W)
solve(A, U.vector(), b)
(u, c) = U.split(True)
print u(0,1)
So, when I'm in, say iPython types : %run tutorial.py, several times, I get a
different value at point (0,1) each time.
--
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
