dolfin team mailing list archive

[Achim Schroll <question107272@xxxxxxxxxxxxxxxxxxxxx>]

New question #107272 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/107272

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()

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.