dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Apr 12, 2011 at 08:54:38AM -0000, Ole Elvetun wrote:
> 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)

The problem is in this part. It's very likely that (as a result of
round-off errors) you are not marking all cells of the mesh. Some will
be marked 0, 1 or 2 and the rest will have "random" values.

Replace the above by

subdomains.set_all(0)
subdomain1.mark(subdomains, 1)
subdomain2.mark(subdomains, 2)

and it should work.

--
Anders


> 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.
>