dolfin team mailing list archive

[Zeljka Tutek <question107272@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Zeljka Tutek proposed the following answer:
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

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