Creating a FunctionSpace changes Mesh

[havogt]

What's the design rationale for having the halo on the function space constructor and not on the mesh constructor, since it changes the mesh? It's quite unexpected.

My pattern is:

atlas::StructuredGrid structuredGrid = atlas::Grid("O32");
atlas::MeshGenerator::Parameters generatorParams;
generatorParams.set("triangulate", true);
generatorParams.set("angle", -1.0);
auto mesh = atlas::StructuredMeshGenerator generator(generatorParams);
atlas::functionspace::NodeColumns fs_nodes_(mesh, atlas::option::levels(nb_levels)
     /*| atlas::option::halo(1)*/); // <--

std::cout << mesh.nodes().size() << std::endl;

[wdeconinck]

There are a number of mesh::actions that need to be applied that may be required only when creating the functionspace::NodeColumns
Other mesh actions may be required when creating functionspace::EdgeColumns
I understand it is unexpected. I have struggled to find a good way.

At the point of mesh-creation you may not yet have the logic to what you want to do with it, as that may depend on a configurable numerical method chosen.

One should use NodeColumns::size(), which would be what you expect.
Creating subsequently another NodeColumns with smaller halo should also return what you expect.

[havogt]

Maybe a function like

atlas::functionspace::NodeColumns enableNodeColumns(mesh, [options]);

would make it more obvious that the mesh changes.

At the point of mesh-creation you may not yet have the logic to what you want to do with it, as that may depend on a configurable numerical method chosen.

Are there useful operations you could already do with the mesh, before you have a functionsspace on it? Why do you want to create the mesh before you have all required information on how you want to construct it?
Of course, the code I am dealing with are very simple standalone examples, so I don't see the complication you might encounter when implementing a full model.

[wdeconinck]

I have considered your proposed function enableNodeColumns versus a NodeColumns constructor, and in the end the API would be the same, bar a different name.
The paradigm of having it in a constructor is that the functionspace object is always created correctly and does not depend on a disjoint creator function.

The mesh should be seen as a more heavy mutable datastructure in support for various function spaces. The Finite Volume scheme e.g. requires 2 functionspaces: NodeColumns, functionspace::EdgeColumns.
Atlas provides a 'fvm::Method', which taken a mesh and options, constructs both function spaces within its Method object. Again the mesh will be adapted.
I may have to document all this better and more clearly

All this said, the functionspace::NodeColumns::size() function, when created with the right halo option, should return you just the number of points in the mesh including the right halo.
Halo's in the mesh are currently always constructed incrementally so node numbering for each halo is following the previous halo nodes. This allows for multiple functionspace::NodeColumns instances with different halos to use the same mesh

[havogt]

I have considered your proposed function enableNodeColumns versus a NodeColumns constructor, and in the end the API would be the same, bar a different name.

There is a big difference: enableNodeColumns documents its side-effect in the name. A side-effect on an argument of a constructor is maximally hidden.

All this said, the functionspace::NodeColumns::size() function, when created with the right halo option, should return you just the number of points in the mesh including the right halo.
Halo's in the mesh are currently always constructed incrementally so node numbering for each halo is following the previous halo nodes. This allows for multiple functionspace::NodeColumns instances with different halos to use the same mesh

I think I understand the concept of functionspace much better now. I was thinking of it as an object which just creates fields. But it is more, it also defines the iteration space on a mesh. In my use-case I don't want to create fields, just iterate over all nodes and access neighbors. This should be done via the function space, not the Mesh.nodes(), right?

If that's right, maybe Mesh.nodes() (and all other methods of Mesh which expose the side effect of the functionspace constructor) could be made private? Then I would be less concerned with this pattern.

[havogt]

And what about the functions to get metric fields, like Mesh.nodes().field("dual_volumes")? Can I get this field via a functionspace? Because with my new understanding of functionspace, I should be able to get this field with different number of halos, depending on the function space.