locally defined source returns zero?

Asked by Achim Schroll on 2010-04-13

Hi!
Why does my locally defined source "g" return just zero, while the same source integrated over the whole domain (replace "dx1" by "dx" in line 62) gives a nontrivial result? (The actual equations are simplified and made up from some other "real problem".)
Any suggestions how to include a local source? Achim.

from dolfin import *
import sys

logging(False)

# linear, coupled system; implicit Euler

k = 0.1 # time step
n = 10 # number of timesteps

# create mesh and function space
mesh = UnitSquare(20, 20)
V = FunctionSpace(mesh, "CG", 1)
W = V + V

# define the source domain
class WellDomain(SubDomain):
    def inside(self, x, on_boundary):
        return bool((0.3 <= x[0] and x[0] <= 0.7 and 0.3 <= x[1] and x[1] <= 0.7))

# mark the subdomain for the well
well = WellDomain()
subdomains = MeshFunction("uint",mesh, mesh.topology().dim())
well.mark(subdomains, 1)

# define integrals over the well domain
dx1 = dx(1)

# area of the well
f1 = Constant(mesh, 1.0)
wellarea = assemble(f1*dx1, cell_domains=subdomains, mesh=mesh)

print ""
print "test solver: dt = %g, n = %g, wellarea = %g" %(k, n, wellarea)

# define trial & testfunctions
(phi, psi) = TestFunctions(W)
(v, w) = TrialFunctions(W)

# define primal functions
s1 = Function(V); h1 = Function(V)#; f = Function(V)

# set initial condition
s0 = project(Constant(mesh, 0.0), V)
h0 = project(Constant(mesh, 0.0), V)

# define source term
temp = Expression("sin(pi*x[0])", V=V)
g = project(temp, V)

t = 0.0

plot1 = plot(s0)
plot2 = plot(h0)

# go timestepping
for i in xrange(1,n+1,1): # i=1,2,...,n

# implicit Euler:
    F = v*phi*dx - s0*phi*dx + k*( (v-w)*phi )*dx \
      + w*psi*dx - h0*psi*dx + k*( (w-v)*psi )*dx \
      + k*g*psi*dx1

# sort out left- & right hand side
    L = lhs(F); R = rhs(F)

# construct the variational problem
    problem = VariationalProblem(L, R)

# solve the discrete problem
    (s1, h1) = problem.solve().split()

    t = t+k

    print "t = %g, |s| = %g, |h| = %g" \
          %(t, sqrt(assemble(s1*s1*dx, mesh=mesh)), sqrt(assemble(h1*h1*dx, mesh=mesh)) )

    s0.assign(s1)
    h0.assign(h1)

    plot1.update(s0)
    plot2.update(h0)

interactive()

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Zeljka Tutek
Solved:
2010-04-14
Last query:
2010-04-14
Last reply:
2010-04-14
Best Zeljka Tutek (zeljkat) said : #1

Achim!
I think that you should define the coefficient m (which enters k*g*psi*m*dx) like this
...
import numpy
...
V0 = FunctionSpace(mesh, 'DG', 0)
m = Function(V0)

class WellDomain0(SubDomain):
    def inside(self, x, on_boundary):
 return True if (0.3 <= x[0] and x[0] <= 0.7 and 0.3 <= x[1] and x[1]<= 0.7) else False
class WellDomain1(SubDomain):
    def inside(self, x, on_boundary):
 return False if (0.3 <= x[0] and x[0] <= 0.7 and 0.3 <= x[1] and x[1]<= 0.7) else True
# mark the subdomain for the well
well0 = WellDomain0()
well1 = WellDomain1()
subdomains = MeshFunction("uint",mesh, mesh.topology().dim())
well0.mark(subdomains, 0)
well1.mark(subdomains, 1)

m_values = [1.0, 0.0] # values of coefficient m in the two subdomain
help = numpy.asarray(subdomains.values(), dtype=numpy.int32)
m.vector()[:] = numpy.choose(help, m_values)
# define integrals over the well domain
dx1 = dx(0)
dx2 = dx(1)

# area of the well
f1 = Constant(1.0)
area0 = assemble(f1*dx1, cell_domains=subdomains, mesh=mesh)
area1 = assemble(f1*dx2, cell_domains=subdomains, mesh=mesh)
....
Zeljka

Thanks Zeljka Tutek, that solved my question.