# problems with plotting subdomains

Asked by Nina S. on 2013-09-13

I created a divided domain with the stokes-equation in the first subdomain and the mixed-poisson-equation (darcy) in the second subdomain.
But now i get the following error when i want to plot the solution:

UMFPACK problem related to call to solve
Warning: UMFPACK reports that the matrix being solved is singular
...
at the end:
assert vmax>=vmin, "empty range, please specify vmin and/or vmax"
Assertion error: empty range, please specify vmin and/or vmax

Can anyone help?
Thanks!

Here is the code:

#-*- coding: utf-8 -*-
from dolfin import *
import numpy as np

# Define mesh
mesh = UnitSquare(32,32)

#Subdomain 1
# Gitter übergeben
subdomains = CellFunction("uint", mesh)

# Klasse des Teilgebiets
class Domain_1(SubDomain):
def inside(self, x, on_boundary):
return between(x, (0, 0.5)) # Koordinatenangabe des Teilgebiets

# Objekt der Klasse erstellen
sub_domain1 = Domain_1()
sub_domain1.mark(subdomains,0)

# Definition Funktionenräume
U = FunctionSpace(mesh, "CG", 2)
V = FunctionSpace(mesh, "CG", 1)
W = U*V

# Definition Trial- und Testfunktion
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Randbedingungen
p_in = 1
p_out = 0
noslip = DirichletBC(W.sub(0), (0),
"on_boundary && \
(x <= DOLFIN_EPS | x >= 0.5-DOLFIN_EPS)")

inflow = DirichletBC(W.sub(1), p_in, "x <= 0.0 + DOLFIN_EPS*1000")
outflow = DirichletBC(W.sub(1), p_out, "x >= 0.5 - DOLFIN_EPS*1000")

bcp = [noslip,inflow, outflow]

# Definition f
f = Expression("0")

# Variationsformulierung
L = inner(f,v)*dx(0)

# Lösung berechnen
w = Function(W)
solve(a == L, w, bcp)
(u, p) = w.split()

# Subdomain 2
# Gitter übergeben
subdomains = CellFunction("uint", mesh)

# Klasse des Teilgebiets
class Domain_2(SubDomain):
def inside(self,x,on_boundary):
return between(x, (0.5,1.0)) # Koordinatenangabe des Teilgebiets

# Objekt der Klasse erstellen
sub_domain2 = Domain_2()
sub_domain2.mark(subdomains,1)

# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
CG = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([BDM, DG, CG])

# Define trial and test functions
(sigma, u, p) = TrialFunctions(W)
(tau, v, q) = TestFunctions(W)

#Define pressure boundary condition
p_in = 1
p_out = 0
noslip = DirichletBC(W.sub(1), (0),
"on_boundary && \
(x <= 0.5 + DOLFIN_EPS | x >= 1.0-DOLFIN_EPS)")

inflow = DirichletBC(W.sub(2), p_in, "x <= 0.5 + DOLFIN_EPS*1000")
outflow = DirichletBC(W.sub(2), p_out, "x >= 1.0 - DOLFIN_EPS*1000")
bcp = [noslip,inflow, outflow]

# Define f

#f = Expression("0")

f = Expression("10*exp(-(pow(x - 0.5, 2) + pow(x - 0.5, 2)) / 0.02)")

# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(1) + inner(p,q)*dx(1) + u*q*dx(1)
L = f*v*dx(1)

# Compute solution
w = Function(W)
solve(a == L, w, bcp)
(sigma, u, p) = w.split()

# plot
plot(u, axes = True, interactive=True, title = "u")
plot(p, axes = True, interactive=True, title = "p")
plot(sigma, axes = True, interactive=True, title = "sigma")

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-09-13