# UMFPACK reports that the matrix being solved is singular.

Asked by Nina S. on 2013-09-21

I created a divided domain with the stokes-equation in the first subdomain and the mixed-poisson-equation (darcy) in the second subdomain. I work with the UnitSquare and the subdomain 1 should be the interval from 0 to 0,5 and the subdomain 2 from 0,5 to 1.

But now i get the following error:

Solving linear variational problem.
UMFPACK problem related to call to numeric
*** Warning: UMFPACK reports that the matrix being solved is singular.
UMFPACK problem related to call to solve
*** Warning: UMFPACK reports that the matrix being solved is singular.

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)

# boundary condition
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)
problem = LinearVariationalProblem(a, L, w, bcp)
solver = LinearVariationalSolver(problem)
solver.solve()
(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)
problem = LinearVariationalProblem(a, L, w, bcp)
solver = LinearVariationalSolver(problem)
solver.solve()
(sigma, u, p) = w.split()

# plot

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

## Question information

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