Use of MeshFunction for BC yields singular matrix

On Wednesday August 31 2011 15:16:06 Marie Rognes wrote:
> Question #169724 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/169724
>
> Status: Open => Answered
>
> Marie Rognes proposed the following answer:
>
> On 08/31/11 21:50, Allan Leal wrote:
> > New question #169724 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/169724
> >
> > Hello,
> >
> > Below is a code that I slightly changed from the one Marie wrote in:
> > https://answers.launchpad.net/dolfin/+question/169325
> >
> > The code solves a pure Neumann problem with mixed finite element. When I
> > use MeshFunction to mark all the boundary facets and then apply the
> > boundary condition with this MeshFunction instance, the resulting
> > coefficient matrix of the system is singular. What is wrong here?
>
> Short story: add
>
> boundaries.set_all(1) # Initialize all to some value, here 1
>
> after your definition of boundaries.
>
> A MeshFunction("...", mesh, mesh.topology().dim() - 1) represents a map
> from facets (edges in the 2d case) to values. So there is a value for
> every edge in the mesh (not just the edges on the boundary). The values
> are not initialized by default, so it is usually a good idea to
> initialize them all to some value first.
>
> When debugging these things, a usually fruitful approach is to reduce
> the size of your mesh and inspect the quantities involved, try for
> instance:
>
> info(boundaries, True)

Another way to debug meshfunction values is to convert it to a NumPy array.
Various things can be looked at. For instance:

  for marker in set(boundaries.array()):
      print "Num:", marker, sum(boundaries.array()==marker)

You loose the connection to facets in the above example, but it will give you
some information about what the markers are and how many of each.

Johan

> --
> Marie
>
> > Allan
> >
> > -------------------------------
> > from dolfin import *
> >
> > # Create mesh
> > mesh = UnitSquare(32, 32)
> >
> > # Define function spaces and mixed (product) space
> > BDM = FunctionSpace(mesh, "BDM", 1)
> > DG = FunctionSpace(mesh, "DG", 0)
> > R = FunctionSpace(mesh, "R", 0)
> > W = MixedFunctionSpace([BDM, DG, R])
> >
> > # Define trial and test functions
> > (sigma, u, c) = TrialFunctions(W)
> > (tau, v, d) = TestFunctions(W)
> >
> > # Define source function
> > f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) /
> > 0.02)")
> >
> > # Define variational form
> > a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v + c*v + d*u)*dx
> > L = - f*v*dx
> >
> > # Define the boundary of the mesh
> >
> > class Boundary(SubDomain):
> > def inside(self, x, on_boundary):
> > return on_boundary
> >
> > # Define the boundary regions of the mesh
> > boundaries = MeshFunction('uint', mesh, mesh.topology().dim()-1)
> >
> > Gamma = Boundary()
> > Gamma.mark(boundaries, 0)
> >
> > # bc = DirichletBC(W.sub(0), (0.0, 0.0), "on_boundary") # From Marie code
> > # bc = DirichletBC(W.sub(0), (0.0, 0.0), Gamma) #<<<< It works with
> > SubDomain instance! bc = DirichletBC(W.sub(0), (0.0, 0.0), boundaries,
> > 0) #<<<< Don't work!
> >
> > # Compute solution
> > problem = VariationalProblem(a, L, bc)
> >
> > w = Function(W)
> > w = problem.solve()
> >
> > (sigma, u, c) = w.split()

Hi Marie,

Thanks. Just a remark for future readers. Adding the line

boundaries.set_all(1)

right before marking it with any SubDomain instance, e.g.

boundaries.set_all(1) # Initialize all to some value, here 1
Gamma = Boundary()
Gamma.mark(boundaries, 0)

 worked perfectly. It does not work if not in this order.

Thanks Marie Rognes, that solved my question.

On Thursday September 1 2011 00:05:52 Allan Leal wrote:
> Question #169724 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/169724
>
> Allan Leal posted a new comment:
> Hi Marie,
>
> Thanks. Just a remark for future readers. Adding the line
>
> boundaries.set_all(1)

Another remark:

  boundaries = MeshFunction('uint', mesh, mesh.topology().dim()-1)
  boundaries.set_all(1)

can be shortened to:

  boundaries = FacetFunction('uint', mesh, 1)

Cheers!

Johan

> right before marking it with any SubDomain instance, e.g.
>
> boundaries.set_all(1) # Initialize all to some value, here 1
> Gamma = Boundary()
> Gamma.mark(boundaries, 0)
>
> worked perfectly. It does not work if not in this order.

Very nice Johan.

Allan