How can I get the local mesh size in velocity direction?

Question

For a SUPG stabilization I need the local mesh size in velocity direction. How can I get this mesh size?

Answer 1

If you have an isotropic mesh, you can do

h = CellSize(mesh)

Otherwise you will have to look at what we do in the pragmatic anisotropic mesh adaptation library

def mesh_metric(mesh):
# this function calculates a mesh metric (or perhaps a square inverse of that, see mesh_metric2...)
cell2dof = c_cell_dofs(mesh,TensorFunctionSpace(mesh, "DG", 0))
cells = mesh.cells()
coords = mesh.coordinates()
p1 = coords[cells[:,0],:]
p2 = coords[cells[:,1],:]
p3 = coords[cells[:,2],:]
r1 = p1-p2; r2 = p1-p3; r3 = p2-p3
Nedg = 3
if mesh.geometry().dim() == 3:
  Nedg = 6
  p4 = coords[cells[:,3],:]
  r4 = p1-p4; r5 = p2-p4; r6 = p3-p4
rall = zeros([p1.shape[0],p1.shape[1],Nedg])
rall[:,:,0] = r1; rall[:,:,1] = r2; rall[:,:,2] = r3
if mesh.geometry().dim() == 3:
  rall[:,:,3] = r4; rall[:,:,4] = r5; rall[:,:,5] = r6
All = zeros([p1.shape[0],Nedg**2])
inds = arange(Nedg**2).reshape([Nedg,Nedg])
for i in range(Nedg):
  All[:,inds[i,0]] = rall[:,0,i]**2; All[:,inds[i,1]] = 2.*rall[:,0,i]*rall[:,1,i]; All[:,inds[i,2]] = rall[:,1,i]**2
  if mesh.geometry().dim() == 3:
    All[:,inds[i,3]] = 2.*rall[:,0,i]*rall[:,2,i]; All[:,inds[i,4]] = 2.*rall[:,1,i]*rall[:,2,i]; All[:,inds[i,5]] = rall[:,2,i]**2
Ain = zeros([Nedg*2-1,Nedg*p1.shape[0]])
ndia = zeros(Nedg*2-1)
for i in range(Nedg):
  for j in range(i,Nedg):
      iks1 = arange(j,Ain.shape[1],Nedg)
      if i==0:
          Ain[i,iks1] = All[:,inds[j,j]]
      else:
          iks2 = arange(j-i,Ain.shape[1],Nedg)
          Ain[2*i-1,iks1] = All[:,inds[j-i,j]]
          Ain[2*i,iks2]   = All[:,inds[j,j-i]]
          ndia[2*i-1] = i
          ndia[2*i]   = -i

A = scipy.sparse.spdiags(Ain, ndia, Ain.shape[1], Ain.shape[1]).tocsr()
b = ones(Ain.shape[1])
X = scipy.sparse.linalg.spsolve(A,b)
#set solution
XX = sym2asym(X.reshape([mesh.num_cells(),Nedg]).transpose())
M = Function(TensorFunctionSpace(mesh,"DG", 0))
M.vector().set_local(XX.transpose().flatten()[cell2dof])
return M

You can then do

sqrt(dot(dot(v,mesh_metric(mesh)),v)/(dot(v,v)+DOLFIN_EPS))

which will give you the inverse length in the direction of v as given by a steiner ellipsoid around each tetrahedron.

Comments

Thank you very much for your information. So I should install your mesh adaptation library?
I tried to do that, but on Ubuntu I get some errors of the form: 'METIS_PartMeshNodal was not declared in this scope'. I'm grateful for any hints.

---

No, you do not need to install the software. It should be sufficient to copy the bits you need from the adaptivity.py file.

---

Thank you again very much. Your advices have solved my problem.