firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
513 stars 160 forks source link

Extracting DOF #1413

Closed Alexey-Voronin closed 4 years ago

Alexey-Voronin commented 5 years ago

Hello developers,

Is there a way to obtain the degree of freedom (dof) to node number (or coordinate) mappings in firedrake? I have spent some time looking at FunctionSpace object but it's not clear to me how it's done.

Any suggestions?

dham commented 5 years ago

The short answer is FunctionSpace.cell_node_map .

The slightly longer answer is that you almost certainly shouldn’t need to access that, because Firedrake provides mechanisms for doing most things that you might want to do with Function values without manually accessing the maps. Could you please explain what it is you actually need to do?

Alexey-Voronin commented 5 years ago

Thank you for your reply, David.

I agree that firedrake does an excellent job of handling those concerns for you. The reason I would like to have access to that info is that I am working on multigrid coarsening and relaxation techniques. This info comes in handy in understanding what's going on.

Couple more questions - my V.cell_node_map() looks like this - Map(Set((18, 18, 18), 'Cells'), Set((49, 49, 49), 'set_4'), 6, None, 'None_cell_node') I under stand that I can access the first and second set via _iterset and _toset, respecticvely. Both of the sets are of type pyop2.base.Set. What is not clear to me is how to actually access that information. Any advice? pyop2 documentation is rather terse on this topic.

wence- commented 5 years ago

So you can get at the process-local values by accessing the .values property of a map. V.cell_node_map().values. Can you explain a little more about what coarsening and relaxation you want to do, it is possible that there are better ways of doing this than by manually inspecting these data structures.

Alexey-Voronin commented 5 years ago

Is there a way to access coordinates too, not just the values? It seems like this information stored in the pyop2.base.Set, on which I cannot operate without a custom kernel routine.

My research in this topic is still in early phases but we are currently exploring algebraic relaxation and coarsening techniques used for solving coupled systems of PDEs. An example of such would be Braess-Sarazin smoothing and AMG (Prokopenko & Tuminaro, 2017)

wence- commented 5 years ago

Sets don't contain any data, only some size information (and ownership of points). Section 4 of this paper https://arxiv.org/pdf/1501.01809.pdf explains a little bit about the data model.

The coordinate field describing the mesh is available as mesh.coordinates. This contains the dofs (basis coefficients). The relationship between cells of the mesh and the coordinate field is again available through mesh.coordinates.cell_node_map() with the array of indirections available with the .values property.

If what you need is topological relationships, this is quite a fragile way to get at it though.

Alexey-Voronin commented 5 years ago

Thank you for the link!

Could you please clarify how one should go about using mesh.coordinates? It looks like it returns a '.Function containing the coordinates of this mesh', but it's unclear what arguments it expects.

It's Fragile in what sense?

wence- commented 5 years ago

Apologies for the radio silence. I mean "fragile" in the sense that, these coordinates are not always associated with vertices of the mesh (e.g. for high-order meshes, or periodic meshes).

To access the coordinates associated with a cell, you use

c = mesh.coordinates
# relationship between cells and coordinate dofs
mapping = c.cell_node_map().values
# dofs
dofs = c.dat.data_ro
# dofs for a given cell
dofs_cell = dofs[mapping[0, ...], ...]

Skimming the referenced paper, it looks like you need the location of the nodes. For point evaluation functionals (which the Q2-Q1 discretisation has), you can compute this by doing:

mesh = ...
X = SpatialCoordinate(mesh)
Q2loc = interpolate(X, VectorFunctionSpace(mesh, "Q", 2))
# this is a function, the values are
Q2loc.dat.data_ro

Q1loc = interpolate(X, FunctionSpace(mesh, "Q", 1))

I guess since you're doing AMG, you'll have some mechanism for transferring these locations to the selected coarse points.

Does this (finally?) answer the initial question?

Alexey-Voronin commented 5 years ago

Yes, this should do the trick. Thank you for the response.