Passing a vertex list to lagrangec0d1 is no longer possible #220

Activity

1. krcools commented on Oct 6, 2026

   

   krcools

   on Oct 6, 2026

   Owner

   More actions

   This happened after merging the allocation free code for the generation of lagrange functions of any order. I noticed some of the constructors I used have disappeared as well and needed to do some workarounds to get all test passing. @cmuenger, what was your philosophy here? I noticed quite a bit of code was commented out in lagrange.jl. I assume you just worked this until tests passed for you?

   Anyways, I think the constructor taking a list of vertex-indices should be deprecated in favour of the constructor taking a list of nodes (i.e. the output of skeleton(mesh, 0) or any submesh thereof. For orders > 1 this generalises to a pair of nodes, edges, and faces. Such a triple is a proxy for what I think of as a list of global degrees of freedom. I am plannign to build a concept of interpolator based on such a list. This is needed for e.g. the interpolation of traces in terms of a basis and in the context of higher order single trace methods for composite bodies.

   In short:

   • @cmuenger, can you revive (dual)lagrangec0d1(mesh::Mesh{U,3,T}, nodes::Mesh{U,1,T})
   • @paresula : please prefer this over the version that takes vertexlist indices.
2. paresula commented on Oct 7, 2026

   

   paresula

   on Oct 7, 2026

   ContributorAuthor

   More actions

   Alright yes, I can work with that, thanks.

   I'll leave this open for now in case @cmuenger wants to comment.
3. cmuenger commented on Oct 7, 2026

   

   cmuenger

   on Oct 7, 2026

   More actions

   I tried to generalize the construction of the Lagrange basis to an arbitrary order. I wanted to reduce code duplication caused by separate constructions for the lower orders (lagrangec0d1, lagrangec0d2, lagrangecxd0, etc.). Now, there are only two functions (continuous and discontinuous) for construction for each mesh dimension (1D [segments], 2D [triangle], 3D [tetrahedron]). Plus a few convenience functions that handle Dirichlet boundary conditions by removing global degrees of freedom at the boundary.

   As @krcools mentioned, the generalization to order > 1 requires not only a list of active nodes but also a list of active edges. For tetrahedral meshes, it even needs a list of active faces.

   So the new generalized interface for continuous Lagrange basis functions on triangles is:

   function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{U,3},
                       nodes::CompScienceMeshes.AbstractMesh{U,1},
                       edges::CompScienceMeshes.AbstractMesh{U,2}; order=1) where {U}  
   

   where nodes is a "node mesh" with the active nodes (nodes that have a global degree of freedom associated) , edges is an "edge mesh" with the active edges, and it is implicitly assumed that all faces in the mesh are active for order >= 3.

   For convenience and compatibility, there is also a function that takes a boolean flag to construct a Lagrange basis, which assumes the entire mesh boundary is subjected to a Dirichlet boundary condition (i.e., only the interior DOFs are considered active). This function constructs the "node mesh" of interior nodes and "edge mesh" of interior edges from the mesh before calling the generalized construction function.

   function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any, 3}; dirichlet::Bool=false, order=1)

   There is also a convenience function for mixed Neumann and Dirichlet boundary conditions. However, this function expects the edges associated with Neumann boundary conditions as the second argument.

   lagrangec0d1(mesh, bnd) = lagrangec0(mesh, bnd; order=1)

   @krcools

   • We can add more aliases for lower-order Lagrange basis functions for convenience, but they should all call the generalized construction function to avoid excessive code duplication. Almost everything I have commented out is duplicated code.
     Maybe something like this line could work. ( I did not test it)
     lagrangec0d1(mesh::::CompScienceMeshes.AbstractMesh{U,3}, interiornodes::CompScienceMeshes.AbstractMesh{U,1}) = lagrangec0(mesh, interiornodes, skeleton(mesh,1); order=1)

   @paresula

   • For what have you been using the old function? To define boundary conditions or for something more complex.
4. krcools commented on Oct 7, 2026

   

   krcools

   on Oct 7, 2026

   Owner

   More actions

   All in favour for pruning the forrest of constructors. I am not sure, but it might be necessary to go for something like this in the end:

   function lagrangec0(mesh::Mesh{3}, dofs::Tuple{Mesh{1}, Mesh{2}, Mesh{3}}; order=1)
       ....
   end
   
   dofs = (active_nodes, active_faces, active_cells)
   

   There are applications where the cell indices should refer to the cell buffer of mesh, but not all cells in mesh are active (multi-threading and local constructions for dual bases for example). For the lowest order specialisations are needed that do not require stipulation of edge- and face-based dofs.

   And at some point I would like to have constructors like this:

   X, itpl = lagrangec0(mesh, dofs; order=p, interpolator=true)
   
   Phi = planewave(;direction, wavenumber)
   phi = strace(Phi)
   
   phi_h = itpl(phi)::FEMFunction 
   

   where itpl is the interpolator for the basis and phi_h is a FEMFunction, i.e. a pair of coefficients and a FEM space that represents the interpolated version of the trace phi.

   For the dual spaces the story is even more complicated because in this case specifying boundary conditions is not equivalent to including or exclusing boundary supported degrees of freedom. Another argument is required: a mesh that is a submesh of the boundary of the supporting mesh.
5. paresula commented on Oct 7, 2026

   

   paresula

   on Oct 7, 2026

   ContributorAuthor

   More actions

   Thanks to both! The reasoning makes sense, and the need for active edges and faces at higher order is convincing for going this route instead of the vertex index lists.
   I am not blocked: I can build the node submesh myself, and BEAST.subset on a global space gives me the same groups.

   For context on what I needed it for: different equations on different parts of a closed surface, with the total function still continuous across the whole surface. So I split one continuous P1 space into disjoint groups of basis functions selected by vertex — no boundary involved, which is why the boundary-mesh argument does not cover it.

   I tested what @cmuenger proposed and it gives the same space as the removed method. (Same basis functions, only enumerated in a different order, since the node submesh keeps skeleton order while the old method followed the order of the vertex list. So it would be a welcome convenience in BEAST, though I can also live without it.
6. krcools commented on Oct 7, 2026

   

   krcools

   on Oct 7, 2026

   Owner

   More actions

   Sorry for the extra work but I think it is the sound choice in the end. For example, gmsh sometimes generates vertices that do not appear in any of the cell definitions (the center of a sphere is an important example here). Also when considering higher order panels and functions, a vertex based approach becomes a nightmare pretty quickly...