dolfin team mailing list archive

[Ole Elvetun <question152507@xxxxxxxxxxxxxxxxxxxxx>]

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.