Knowing knot positions for a THBSplineBasis

[Milloupe]

Hi all,

I'm using G+Smo to manage THBSplines, in order to avoid using a tensor basis. My main problem is that I could not find a way to know the exact knot positions of a refined THBSpline basis.

What I have done:

• Create rather sparse knot vectors
• Build a tensor basis using these knot vectors
• Refine this basis in specific places

Now I want tot do some interpolation, so I need to know the spline values in many points. I know that for this, the anchors function gives a list of interpolation points, but for my application it is important that the interpolation points are at the same position as the basis knots.

However, when I want to know the basis' knot positions, I can't find a correct way to do it, because the function knot(lvl, dim, i) gives you knot positions along each direction, but since we don't have a tensor basis it is not clear which 2D points are constructed with these 1D knot vectors.

Am I missing something here? Otherwise, I suppose I will just compute by hand what the refined knot positions are.

Thank you very much, Best, Denis

[filiatra]

Hi, what is the the result you wish to have? Can you provide a small example? Say you have a square and one knot at 0.5 at level 0 for u and v, and then level one is the refinement of th element [0, 0.5] in both u and v. Then which are the anchors that you want to get?

[Milloupe]

Hello,

Thank you for the quick answer. You're right, it needs some details. Ideally, I want the anchors to be both the 0-level knots (whose positions I know) and the 1-level knots added by the refinement (whose positions I can't seem to access). So, in your example : (0.5, 0.5) and whatever the knots created by the refinement are.

In the first place, I thought I understood how the refinement took place, and that I could compute the 1-level knot positions by hand, but when I used the anchors function, I saw that there were many more points in the refined zone than I expected.

Why this is complicated in my case is that the zone I refine is the sides of a rectangle, without refining the center, so the refined zone is not convex. This is why the 1-level basis is not a tensor basis, and knowing the knot positions in each direction does not tell me directly which are the added knots.

So my goal is either to :

• find a way for the basis to "tell" me which knots it contains on each level or
• understand the refinement process better so that I compute the refined positions by hand.

However, if this is too complicated, I think I can refine the whole rectangle, at least to see if it does not add too many points in non-significant areas.

Thank you very much, Best, Denis

[filiatra]

A relatively easy way I can see to get the knots/points you want is: Iterate over all the elements of the mesh (there is a domain iterator that does that) and for every element add the upper corner to a list, this would look something like:

    typename gsBasis<>::domainIter domIter = thbbasis.makeDomainIterator();
   for (; domIter.good(); domIter.next() )    // loop over all elements
   {
     // The following provide the lower and upper corner of the current element
     // that you can grab here:
     domIter.lowerCorner();
     domIter.upperCorner();
  }

Note that the points on lower and upper will result in having each point twice, so you would probably use only upperCorner() to get distinct points.

Otherwise: how do you plan to get as many points as the number of basis functions so that you obtain a well-defined interpolation?

[filiatra]

I think this is a nice function to have, if you try it and it works we can add it as a new member function in the spline class.

[Milloupe]

Hello,

Thank you again for the super fast answer. The following function works just fine:

     gsHDomainIterator<real_t, 2> domIter(base);
     std::vector< gsVector<real_t> > lower_corners;
     int nb_points = 0;
    for (; domIter.good(); domIter.next() )    // loop over all elements
    {
     lower_corners.push_back(domIter.lowerCorner());
     nb_points++;
    }

I'm not an expert in c++, so this code might not be up to standards, but it does pretty much what we expect. I think I will need the upper corners too, because of the difference it makes on the upper- and right-most boundaries of the different mesh zones.

Now, I need to compute the splines on the knot positions specifically because I want to have more values than splines, in order to add periodicity conditions. For this, I define the spline domain as larger than the unit cell I want to compute, and then "wrap" the overflowing splines.

More precisely, with an actual example:

Here in blue are all interpolation points (which are also knot positions). In orange are the "overflowing" interpolation points for which I have found a corresponding point within the unit cell (e.g. a point whose position is (x-period_x, y))

As you can see, using only the lower corners leaves blanks around the refined zone, and therefore leads to missing correspondences with the overflowing points, so I'll try using the upper corners as well, and see if it works better.

[Milloupe]

For the sake of completeness, here's the same image but using only upper corners:

[filiatra]

thnx for the explainations, can yopu also show the mesh (it is possible by creating a gsMesh using a constructor: https://gismo.github.io/classgismo_1_1gsMesh.html)

[Milloupe]

I would gladly do it but how can it be represented? When I create it the avaiblable functions are just two different addline's (and displaying it simply shows the number of edges, vertices and faces)

[filiatra]

I would gladly do it but how can it be represented? When I create it the avaiblable functions are just two different addline's (and displaying it simply shows the number of edges, vertices and faces)

Then you can output it to paraview: gsWrite(mesh,"mesh");

[Milloupe]

Sorry for not answering earlier, I do not easily have access to Paraview, it should change soon, I'll keep you updated as soon as possible. In the meantime, I've had another issue with what I was implementing: I have written the part of the code that makes sure to take into account the periodicity, as mentioned above. With that, even if taking all the upperCorner point positions of a non-refined mesh gives more points than there are functions, taking into account the periodicity removes the extra points (by "wrapping" them back to the beginning of my unit cell) and I have a square invertible matrix[point, function]. However, when I try to do the same thing with a refined mesh, even when the refined part is well within my unit cell and has no influence on the "wrapped" points, I find that the refinement adds more points than it adds functions, so I end up with a rectangular matrix[point, function] (more lines, i.e. points, than columns, i.e. functions). (I'll provide images and numbers in following messages) Do you know why this is? And because I know that my explanation is a bit fuzzy, do you have any reference about the refinement process, aside digging in the source code itself, so I can figure out what exactly is happening? Thank you very much in advance.

[Milloupe]

This is the refined mesh I am using, with the colored points being there to show explicitly what is being "wrapped" around due to periodicity

And this is the same type of refined mesh, with all points being in blue, and orange crosses are added to points that I keep to have as many points as functions. To find which points I keep, I simply see, for each point, if the matrix[point, function] is still full-rank when I remove the point. If it is, the point gives degenerate information, and I do not keep it.

[dominik-mokris-mtu]

Hi @Milloupe, is this still relevant or can we close the issue?

[Milloupe]

Hi, sorry for not answering earlier. This was part of a PhD work, and I since moved on a new project, I'm not sure if I'll ever have the time to work on it again; I close the issue. Thak you!

[dominik-mokris-mtu]

This was part of a PhD work

Cool! If you end up using G+Smo for your thesis or papers, let us know, please. It is always nice to read about nice results to which the library has contributed.