This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Constrain 1 DOF

+1 vote

I have a mixed variational problem with a non-trivial null-space (think of the Stokes problem with velocity boundary conditions), which I would like to solve using a direct solver. The degrees of freedom involved in the null-space vector are coming from the DG0 and RT1 elements. In order to remove the null-space I would like to simply constrain one DOF of a DG0 space, which seems to be quite difficult to "catch" using the standard tools. Does anyone know an easy way to do this (which also works with decomposed meshes in MPI environment)?

Here is a simplified example of the issue:

from __future__ import print_function
from dolfin import *

comm    = mpi_comm_world()
mpiRank = MPI.rank(comm)
parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(10,10)
DG0_Element = FiniteElement('DG',mesh.ufl_cell(),0)
V = FunctionSpace(mesh,DG0_Element)
u = TrialFunction(V)
v = TestFunction(V)

# just an example
a = inner(jump(u),jump(v))*dS
L = Constant(1.0)*v*dx

# Question: how to construct a boundary condition in order to remove
# the non-trivial null-space from the problem? 
bcs = []

A,b = assemble_system(a,L,bcs)

null_space = project(Constant(1.0),V)
# test the null-space
tmp = Function(V) 
A.mult(null_space.vector(),tmp.vector())
tmp_norm = norm(tmp.vector(),'l2')
if mpiRank==0:
    print('norm = %e' % tmp_norm)
asked Apr 12, 2017 by ae FEniCS Novice (290 points)

1 Answer

0 votes
 
Best answer

Hi, the following uses the fact that for DG0 term inner(u, v)*dX contributes only to the diagonal of the matrix. You could scale the term by cell volume to get 1 on the diagonal.

cell_f = MeshFunction('size_t', mesh, 0)
if mpiRank == 0:
    cell_f[0] = 1
dX = Measure('dx', domain=mesh, subdomain_data=cell_f, subdomain_id=1)

# just an example
a = inner(jump(u),jump(v))*dS + inner(u, v)*dX
answered Apr 14, 2017 by MiroK FEniCS Expert (80,920 points)
selected Apr 14, 2017 by ae

Great, thanks! Does one need to do anything like 'apply()' after setting the cell-function element on the mpi node 0 because of the ghost mesh?

No, the code as given works in parallel. I'd just like to point out that it is more mathematically sound to seek the solution with zero mean value, rather than pinpointing single dof.

Thanks again! Presently I have the code with Lagrange multipliers implemented, but the system is quite large as is, so I wanted to see, if I could cut corners without having dense rows/columns (or at least rows with many non-zeros) in the matrix.

...