Based on other comments, I know NGsolve can do curved meshes. However, I am having trouble finding an example, or tutorial online. Perhaps it is buried in the examples in the source code. Is that right?

I just want to be able to take a mesh, and define a higher order space on it, e.g. cubic lagrange, and set the nodal values of those triangles myself.

in Netgen, high order curved elements are represented by a hierarchical basis, which is better conditioned in comparison to a Lagrange basis, for high p. It makes it also easy to vary p over the mesh, i.e. straight edges don’t need high order basis functions at all.

Beside vertices, coefficients are associated with edges in in 2D/3D, and faces in 3D. No element-interior coefficients are needed.

In general, we don’t have the same basis (or even the same space) for the element mapping, and the finite element space, i.e. no isoparametric elements.

For order=2, we get the traditional nodal 6-node trigs / 10-node tets using mesh.SecondOrder().

Maybe you consider to go with mesh.SetDeformation(gridfunction) ?

Ah, ok. So my GridFunction in “SetDeformation” can be any H1 space (any order)? That should be fine.

So a (standard) cubic Lagrange element in 2-D has one interior node. But you don’t have that with hierarchical basis?

So is there a way to loop through all triangles that have an edge on the boundary?

The way I have done this in the past, is I use the Lenoir method of finding the optimal curved triangle that maintains convergence rates. But that uses the standard Lagrange element. I’ll need to think how this would work with hierarchical basis. I guess I could (element-wise) define my own Lagrange element (of whatever order) and interpolate that into the NGsolve space.

So a (standard) cubic Lagrange element in 2-D has one > interior node. But you don’t have that with hierarchical > basis?

with the same order for all nodes we have the cubic bubble, but we have the freedom to switch it off.

Please let us know if you get better results with your geometry interpolation method, in comparison to the projection-based interpolation we use. The projection-based interpolation is optimal in h and p, but we are simplifying a bit here and there.

You switch it off? But you need that interior node. That’s kind of the point of that Lenoir paper.

Yes, correct, we are not using interior basis functions for the element mapping.

Due to the hierarchical basis, we have a smooth blending from the boundary into the element.

This is different from Lagrange elements as in Lenoir’s 86 SINUM paper.

If you have an implementation of Lenoir’s method at hand, I would be curious to see some comparisons.

My bet is that the results will be very similar. In case we are worse, I am happy to learn how to improve NGSolve.