Shape function error for quadratic hexahedral elements?

[schwarz-e]

I produced a cylinder out of quadratic hexahedral elements and attempted to run a simple inner wall pressure condition simulation. However, on the post-processing step it exits with a "ERROR: Error in computing shape functions" which appears to stem from the GETNNX subroutine.

Is there a way to work around this? Is this an issue with the simulation (appears to converge) or just the post-processing step? I've attached the input files below for reference.

test.zip
?

[vvedula22]

@schwarz-e The mesh looks weird to me. It seems to have only 3 cells but the connectivity appears to have been messed up. How did you create the hex mesh? If you have a trilinear brick mesh, you can use the mesh-converter tool to convert to either a 20-noded quadratic brick or a 27-noded triquadratic brick element.

[schwarz-e]

I created it using a pyvista script, so it doesn't start as a trilinear brick mesh. It's a 20-noded quadratic brick. This example only has 3 cells (was testing a simple case). What's weird about the connectivity? I followed the numbering system outlined in the vtk documentation: https://raw.githubusercontent.com/Kitware/vtk-examples/gh-pages/src/Testing/Baseline/Cxx/GeometricObjects/TestIsoparametricCellsDemo.png

[schwarz-e]

Looking closer, it appears the calculated local coordinates calculated in GETXI via Newton iterations converges to a value that is not within the element bounds which causes this issue.

[vvedula22]

@schwarz-e Yes, that happens because you have a bent geometry around the edge nodes of the quadratic brick element. During post-processing, we project pressure field (and perhaps some other quantities) to trilinear basis and interpolate the solution at the edge nodes using trilinear basis functions. This is done to avoid a checkered/spotty field where you will see high values on the edge nodes compared to the corners.

In the current case, while projecting to trilinear brick element, the bent geometry is lost and you will project onto a prismatic geometry (triangular cross-section). During the interpolation step, however, the edge nodes are outside the `triangular' element, and therefore, the code cannot find a suitable basis to interpolate. Is there a reason for choosing this bent configuration?

[schwarz-e]

Yes, the reason is to best approximate curved reference geometry. In addition, it is part of my formulation where I apply calculated displacement solutions onto the mesh and then run that as the reference geometry in the next iteration, so even if mesh starts out as brick elements with nodes added to make it quadratic, it will eventually travel to become bent.

[vvedula22]

You can either refine the mesh and add more cells, or comment the section of the post-processing code (POST.f) that does the projection and interpolation steps on certain quantities like pressure, stresses, etc. Let me know if you run into any issue.

[schwarz-e]

I am also having issues with the GETNNX call in SETBC.f in the SETBCDIRWL. I found that using a weakly couple dirichlet boundary condition was the best way to impose axially clamped (effective direction) motion at inlet faces. Is there a good workaround for this?

[schwarz-e]

With Dirichlet BCs on the HEX20 faces, depending on if effective direction is applied the flow direction flips, and in strongly applied effective direction on the dirichlet face the flow profile does not remain parabolic and instead warps with the motion of the inlet face.

[vvedula22]

Are you doing a fluid or FSI simulation? Weakly coupled Dirichlet BCs are implemented for fluid domains only. May not work for solid domains or if you are running a pure solid mechanics simulation. The current implementation doesn't not handle virtual stress for hyperelastic materials. Might only work for small strain linear elastic materials.

Yes, it is not surprising. The derivatives are computed using the Newton method and therefore, fail in your case that has edge nodes lying outside of the main element.

[vvedula22]

With Dirichlet BCs on the HEX20 faces, depending on if effective direction is applied the flow direction flips, and in strongly applied effective direction on the dirichlet face the flow profile does not remain parabolic and instead warps with the motion of the inlet face.

The effective direction doesn't take the face normal into account. In a regular simulation (without effective direction), the sign of the flow is computed based on the sign of (u.n). The user should then provide a value that makes physical sense (negative value for inflow, positive value for outflow). On the other hand, if you use effective direction, the code doesn't use normal anymore. It simply uses the user-provided value in that direction.

If you don't want the flow to warp with the motion of the inlet face (which typically happens in FSI simulations), you should fix the inlet and outlet faces for the mesh equation.