Add 1D Higher Order Lagrange Functions (Cx and C0)#175
Conversation
Polynomials are implemented hard-coded up to order 10 for maximal speed. TODO: The Mie unit test should be adapted once curvilinear meshes are available. So far, we only confirm that higher-order functions do not break the performance. However, since the mesh is not refined the geometrical error starts to dominate quickly.
sbadrian
force-pushed
the
feature/higher_order_4
branch
from
6cc104b to
b1dc8dc
Compare
Greatly improves the speed. I expect further improvements are possible, but this was the most significant type instability that I have found.
| function lagrangec0_dirichlet(mesh; order=1) | ||
|
|
||
| verts = skeleton(mesh, 0) | ||
| detached = trues(numvertices(mesh)) |
There was a problem hiding this comment.
skeleton(mesh, 0) gives the nodes (vertices that are acutally used in the cell definitions). So it is not necessarily true that numvertices(mesh) == length(skeleton(mesh,0)). The latter should be used here
There was a problem hiding this comment.
Are you suggesting to replace detached = trues(numvertices(mesh)) by detached = trues(length(skeleton(mesh,0)))? I think it is correct the way it is, because the goal in the following is to filter out the detached vertices (I am following here the same recipe as at other places in lagrange.jl and I have verified that it works as desired).
| notonbnd[v] = false | ||
| end | ||
|
|
||
| vertexlist = findall(notonbnd .& .!detached) |
There was a problem hiding this comment.
I think this can be expressed as:
setminus(skeleton(mesh,0), skeleton(boundary(mesh),0))
There was a problem hiding this comment.
I agree. In fact, the entire function could be replaced by
function lagrangec0_dirichlet(mesh; order=1)
internalnodes = setminus(skeleton(mesh,0), skeleton(boundary(mesh),0)).cells
internalnodes = sort!([node[1] for node in internalnodes])
lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order)
end
Given that the same mechanism is used in the other lagrange functions, I could replace it everywhere in lagrange.jl.
There was a problem hiding this comment.
Let's leave this all as is for now. I'll make a testitem for the case where I think the two behaviours might differ.
There was a problem hiding this comment.
My currents view is to hide as much as possible the vertices and use the 0-skeleton instead. Sometimes meshes like gmsh add floating (=unused) vertices to the file and so the two concepts are not always interchangeable.
There was a problem hiding this comment.
Yes, definitely, but I think the function was written precisely in this way to filter out the unused vertices when finding the internal nodes.. I think though that the setminus approach you have proposed is indeed much more elegant and clearer. I think we just should do it consistently like this in lagrange.jl as otherwise one would wonder why there are different approaches are used.
Of course, this could be done in a separate PR as this is a general refactoring.
krcools
left a comment
krcools
left a comment
There was a problem hiding this comment.
Dont you need implementations for interpolate in laglocal.jl?
|
It seems that for the momintegral routines the |
|
For some quadrature strategies they are almost certainly required. I even though that the vertex reorderings that are sometimes needed for Sauter-Schwab quadrature to work rely on them. For non-conforming meshes and refined meshes they are also used in general. I suggest we leave as is and implement the first time someone tries one of these more complicated scenarios. |
|
Sounds good to me. I think then the branch is ready to be merged. |
2 participants
Polynomials are implemented hard-coded up to order 10 for maximal speed.
TODO: The Mie unit test should be adapted once
curvilinear meshes are available. So far, we only
confirm that higher-order functions do not break the performance. However, since the mesh is not refined the geometrical error starts to dominate quickly.