NGSolve / ngsolve

Netgen/NGSolve is a high performance multiphysics finite element software. It is widely used to analyze models from solid mechanics, fluid dynamics and electromagnetics. Due to its flexible Python interface new physical equations and solution algorithms can be implemented easily.
https://ngsolve.org/
GNU Lesser General Public License v2.1
436 stars 79 forks source link

Vector Field from 2D meshes Not Properly Recognized In ParaView #43

Closed Alex-Vasile closed 2 years ago

Alex-Vasile commented 2 years ago

If I run a 2D solve and export my vector field in .vtu format, ParaView will not allow me to use it as the input for vector filters.

I can see the x and y components of my vector field, as a well as the magnitude, but I cannot use it as the input for e.g. the Glyphs filter.

What I have to do every time is create a Calculator filter which takes my vector u as input are returns:

v = u_x * iHat + u_y * jHat.

I think it is because the vector field being outputted from NGSolve is missing a z component (or whichever component is not on the 2D plane that's being solved.

schruste commented 2 years ago

Yes. Your observation is right.

The question is of course if you see this as a flaw of the VTK-output or a flaw of ParaView. You can simply extend your CoefficientFunction with a zero when exporting so that ParaView obtains a 3D vector or you do the Calculator-filter on the paraview side. The VTKOutput is general w.r.t. the dimension of fields, so that I don't want to generically fill up a 2D field with a zero in the third component.

Best, Christoph

Alex-Vasile commented 2 years ago

Hi Christoph,

You're right about whether this is a vtk-ouput or a praview issue. I agree that it's a paraview issue.

There is even an issue open for it on their issue tracker https://gitlab.kitware.com/paraview/paraview/-/issues/12677

However, that issue is 10 years old. So I do not have high hopes of it being fixed anytime soon.

So, while I agree this problem is not cause by NGSolve, I am still asking if NGSolve could add the functionality to fill in those zeros for vtk output as a compatability feature.

Alex

schruste commented 2 years ago

Hi Alex,

I don't like exporting so many additional zeros where the data is obviously not necessary. It increases the file sizes unneccessarily. A simple fix for you could be that you simply put a wrapper around your CoefficientFunctions:

def Make3DVec(u):
  if u.dim == 2:
    return CF((u[0],u[1],0))
  else:
    return u

This makes sure that you never export 2D fields to VTK, but it still allows me to only export 2D fields when I don't want to have the additional fill-up with zeros.

Best, Christoph

Alex-Vasile commented 2 years ago

The problem I have with that solution is that it assumes that u is on the x-y plane. Does u, as a gridfunction/coefficientfunction, store the relevant information to figure out which plane (x-y, y-z, or x-z) it lives one?

schruste commented 2 years ago

Well, it's a 2D vector field, it only has the two components u[0] and u[1]. The vector field u originally does not know that it should be embedded into a 3D space.

The workaround that I posted before adds a zero z-component, resulting in a 3D vector that lives the x-y plane, but you could of course also put the zero somewhere else ending up in a vector field that is in the y-z or x-z plane.

Alex-Vasile commented 2 years ago

Right.

Perhaps I didn't explain myself properly.

I wish to write this Make3DVec function to be general and able to support all 3 planes.

To do that, I need to know which plane u is on as it will change which how I construct the 3D coefficient function.

Does the grid function/coefficient function provided to Make3DVec contain enough information to for my to identify which of the three planes it lies on?

schruste commented 2 years ago

No. The CoefficientFunction is only 2D. It does not know anything about its potential embedding into 3D. The user needs to provide that information in this case. I am actually not sure which use case you are talking about where the CF is 2D but you need the third component to interprete the vector field.

Alex-Vasile commented 2 years ago

Is this closed as "won't address", or closed as "addressed with a commit"?

ChrLackner commented 2 years ago

This is closed as the behavior is explained and seems to default to standard vtk usage. We won't write a lot of 0 entries (although you can if you want, as Christoph explained) into the vtk output as it can be simply fixed with a paraview filter. Also if you have a 2D CF you have to somehow define yourself what the plane embedding into 3D space is. I do not see an open question here.