DOLFIN

locally defined source returns zero?

Asked by Achim Schroll

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:
Last query:
Last reply:
Revision history for this message
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

Revision history for this message
Achim Schroll (achim-simula-deactivatedaccount) said : 		#2

Thanks Zeljka Tutek, that solved my question.