Physical location of basis functions

[bensworth]

Suppose I have a 2d DG function space with nodal basis. How can I get the global geometric location of the matrix DOFs?

In this particular instance, I want to construct a sparse matrix that maps a DG vector to a continuous vector by averaging the overlapping DOFs of adjacent elements. Of course this can be done with FEM projecitons, but for the preconditioning concepts I'm working on I really need it in operator form and as sparse as possible.

FWIW this has come up before as well for me; geometric information can be useful in developing preconditioners (which is a lot of what I work on), but it seems very hard to access this information?

Thanks!

[ScottMacLachlan]

For any space with a nodal basis, can't you just interpolate the spatial coordinates onto that basis? Something like

(x,y) = SpatialCoordinate(mesh)
V = FunctionSpace(mesh,"DG",2)
xfunc = Function(V).interpolate(x)
yfunc = Function(V).interpolate(y)

Then you can look at xfunc.dat.data and yfunc.dat.data to get the coordinates. I don't know how robust that is, but I've done this before with CG spaces and gotten what I needed from it.

[bensworth]

That's funny, thanks Scott. I had more or less tried that and had a mountain of errors, I was thinking as due to the DG space. The problemw as actually due to a VectorFunctionSpace, which the above does not work on. For DG it is fine.

[wence-]

Hmm, interpolate(SpatialCoordinate(mesh), VectorFunctionSpace(mesh, "DG", 1)) works fine for me.

[wence-]

Note that this approach of interpolateing the SpatialCoordinate works for any space where the nodes are point evaluations. If nodes are moments (for example BDM) then the geometric location isn't really very well defined, I think.

[bensworth]

Hm, I suppose this is my confusion though; neither of the following work:

X = interpolate(SpatialCoordinate(mesh), FunctionSpace(mesh, "DG", 1))
*** RuntimeError: Expression of length 1 required, got length 2

(x,y) = SpatialCoordinate(mesh)
xfunc = Function(VectorFunctionSpace(mesh,"DG",degree)).interpolate(x)
yfunc = Function(VectorFunctionSpace(mesh,"DG",degree)).interpolate(y)
*** RuntimeError: Expression of length 2 required, got length 1

although both of the following work:

X = interpolate(SpatialCoordinate(mesh), VectorFunctionSpace(mesh, "DG", 1))

(x,y) = SpatialCoordinate(mesh)
xfunc = Function(FunctionSpace(mesh,"DG",degree)).interpolate(x)
yfunc = Function(FunctionSpace(mesh,"DG",degree)).interpolate(y)

[wence-]

So SpatialCoordinate(mesh) is a representation of the physical coordinate at the quadrature points, so it has shape (geometric_dimension,). If you want to interpolate it somewhere, then you need to be interpolating into a functionspace with matching shape. FunctionSpace(mesh, "DG", 1) gives you a scalar function space (with shape ()). VectorFunctionSpace(mesh, "DG", 1) gives you a vector-valued function space with geometric_dimension copies of a scalar space at each node (this is what you need for interpolating the spatial coordinates).

Conversely, when you destructure x, y = SpatialCoordinate(mesh), x is the X component of the physical coordinates (with shape ()) and similarly for y. So now you can't interpolate one of those into a VectorFunctionSpace (because that needs an expression with shape (geometric_dimension,) and you're providing (). Make sense?

[bensworth]

In case anybody is interested, here is code that constructs sparse scipy CSR matrices to average across DG interfaces in vector form, as well as a projection onto the space of continuous functions in a discrete sense. Note, I have not accounted carefully for BCs here.

# Get ordered nodes of mesh
(x,y) = SpatialCoordinate(mesh)
Vdg = FunctionSpace(mesh,"DG",degree)
xfunc = Function(Vdg).interpolate(x)
yfunc = Function(Vdg).interpolate(y)
nodes = np.vstack((xfunc.dat.data.T, yfunc.dat.data)).T
num_dg_dofs = nodes.shape[0]

# Get sorting of mesh nodes by x and then y
inds1 = nodes[:,0].argsort()
inds2 = nodes[inds1,1].argsort(kind='mergesort')
sort_inds = inds1[inds2]

# Find repeated nodes; repeated nodes are zero when we take the
# difference of nodes_sorted[1:] - nodes_sorted[0:]. The groups
# of indices are the intervals in node_diff, e.g., sucessive
# entries [4,10] in node_diff -> sort_inds[4:10] correspond to
# matching geometric nodes.
node_diff = 1 + np.where( np.sum(np.abs(nodes[sort_inds[1:]] - nodes[sort_inds[0:-1]]),axis=1) > 0)[0]
if node_diff[-1]==num_dg_dofs:
    node_diff = np.concatenate(([0],node_diff))
else:
    node_diff = np.concatenate(([0],node_diff,[num_dg_dofs]))
num_unique = len(node_diff)-1

# Construct CSR matrices to restrict and interpolate between 
# discontinuous space by averaging over shared nodes, or
# interpoalting by value to all shared nodes.
data = np.ones((num_dg_dofs,))
R = sp.csr_matrix((data,sort_inds,node_diff),shape=[num_unique,num_dg_dofs])
P = R.T

# Normalize rows of R; might be better way to do this, but
# this should be pretty efficient
row_sums = 1.0*(node_diff[1:] - node_diff[0:-1])
row_indices, col_indices = R.nonzero()
R.data /= row_sums[row_indices]

# Form projection onto R(P)
diag = R*P
diag.data = 1.0/diag.data
proj = P*diag*R

[colinjcotter]

The matrix action is available via the Projector class.

…

________________________________ From: Ben Southworth ***@***.***> Sent: 17 February 2022 16:44 To: firedrakeproject/firedrake ***@***.***> Cc: Subscribed ***@***.***> Subject: [firedrakeproject/firedrake] Getting coordiantes of matrix DOFs (Discussion #2354) Suppose I have a 2d DG function space with nodal basis. How can I get the global geometric location of the matrix DOFs? In this particular instance, I want to construct a sparse matrix that maps a DG vector to a continuous vector by averaging the overlapping DOFs of adjacent elements. Of course this can be done with FEM projecitons, but for the preconditioning concepts I'm working on I really need it in operator form and as sparse as possible. FWIW this has come up before as well for me; geometric information can be useful in developing preconditioners (which is a lot of what I work on), but it seems very hard to access this information? Thanks! — Reply to this email directly, view it on GitHub<#2354>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4SN77MICFXXR4VW4H3U3UQW7ANCNFSM5OVHM2MA>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[wence-]

@colinjcotter that computes the L2 projector which is not what @bensworth wants, I think. I thought you would be able to construct P like so:

V = FunctionSpace(mesh, "P", degree)
multiplicity = Function(V)
for c in multiplicity.cell_node_map().values:
    multiplicity.dat.data_with_halos[c] += 1

intp = Interpolator(TestFunction(V)/multiplicity, Vdg)
Pprime = intp.callable().handle # PETSc matrix represented the averaging prolongation

The rows all come out with right values, AFAICT, but the numbering is not the same as @bensworth's code.

[wence-]

e.g.:

In [65]: R.todense().round(2)
Out[65]: 
array([[0.  , 1.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.33, 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.33, 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.17, 0.  , 0.  , 0.17, 0.17, 0.  , 0.  , 0.  , 0.17, 0.  , 0.17, 0.  , 0.  , 0.  , 0.17, 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.33, 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.33, 0.  , 0.33, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 1.  ]])

In [66]: intp = Interpolator(TestFunction(Vdg)/f, V)

In [67]: intp.callable().handle[:,:].round(2)
Out[67]: 
array([[0.  , 1.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.33, 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.33, 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.17, 0.  , 0.  , 0.17, 0.17, 0.  , 0.  , 0.  , 0.17, 0.  , 0.17, 0.  , 0.  , 0.  , 0.17, 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.33, 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.33, 0.  , 0.  , 0.33, 0.  , 0.33, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 1.  ]])

[wence-]

That might be because @bensworth sorts via the coordinates, so the inverse permutation applied to the rows and columns would probably get things out.

[bensworth]

Interesting, thanks @wence- . It seems like Firedrake can usually do what I want, but I still have to bug you guys to figure out how. The ordering is a bit puzzling, I sort the coordinates to find overlap, but when constructing P/R I use the original ordering (i.e., the ordering in my code I think should be right, and I have confirmed that it is definitely the kernel operator I am looking for). Will play with this some more.

[wence-]

I guess it's not fully obvious that an Interpolator can deliver back the operator (rather than just the action). Hopefully many of these things will improve with the addition of symbolic-level proper dual space support (and an honest-to-goodness Interp operator in UFL)