Calculating stiffness matrix for specific cells in mesh

[jmm1828]

I have a simple 2D mesh, say mesh=UnitSquareMesh(2,2), and I want to calculate the stiffness matrix for the bilinear form a(u,v) in my linear variational problem at different cells in the mesh. I can calculate a "global" stiffness matrix as follows:

# create mesh
m = UnitSquareMesh(2,2)
# Create piecewise linear function space
V=FunctionSpace(m,"CG",1)
# Create test and trial functions
v=TestFunction(V)
u=TrialFunction(V)
# create bilinear form a(u,v) = ∫ uv dx
a = u*v*dx
# calculate stiffness matrix
assemble(a)
...

But is there a way to perform this calculation cell by cell? Any help would be much appreciated.

[wence-]

If you make a BrokenElement version of your finite elements then it will have discontinuous (cell-wise) continuity and that will allow you to get cell-wise matrices. Although it won't work if you have facet integrals I think.

But, what do you want to do with these cell-wise matrices? There might well be a better way to achieve what you want.

[jmm1828]

Thank you for this, and would it be possible to assign weights to different edges with this method?
Say for a given cell each edge has a weight and once you calculate an outer product you multiply the matrix by this weight before adding.

[wence-]

This is a little trickier since you need to be able to associate the edge weights with the appropriate cell edge vector. How are you generating and storing the edge weights at the moment?

[jmm1828]

I have not generated the weights yet, they would depend on the length of the edge and the shape regularity of the mesh. I wanted to try some arbitrary weights first

[wence-]

Arbitrary weights are in some sense harder than depending on length and shape regularity since you need to be able to associate them with the correct edges and then pull them out again (which is not trivial because if you store them as a Function they will be evaluated at quadrature points whereas you want to "see" the basis vector).

On the other hand, the length of each relevant edge can be computed with sqrt(dot(cev[0, :], cev[0, :])) (for edge 0, say). I do not think there is a symbolic expression for the shape regularity coefficient of a given cell, however, this is easier. If you have some way of computing it, you can store it in a shape_regularity = Function(FunctionSpace(mesh, "DG", 0)) and then use that in your form.

[wence-]

I guess if you're using 2D elements and you want the condition number of the Jacobian as your regularity metric then you can use that for a 2x2 matrix $J$ that the eigenvalues are:

$$ \text{tr}(J) / 2 \pm \sqrt{\big(\text{tr}(J) / 2\big)^2 - |J|} $$

So

J = Jacobian(mesh)
detJ = JacobianDeterminant(mesh)
m = tr(J) / 2

kappa = (m + sqrt(m**2 - p))**2 / p

Edit: if you go this route, you'll probably want to manually select an appropriate quadrature degree for the stabilisation term (...*dx(degree=degree_of_polynomial_you_want_to_integrate_exactly)) since UFL will otherwise estimate something that is too high.

[jmm1828]

Thank you for this. I just wanted arbitrary weights for testing purposes, but this is probably a better route.