2nd order Lagrange Prism in writing with pumi_mesh_write

[eisungy]

Dear SCOREC colleagues,

Hi =)

While I was trying to create a PRISM type curved element, I encoutered a seg fault in writing it into VTK using pumi_mesh_write.

My source code looks as below. I couldn't see a seg fault with PUMI_TET, and was able to see the curved tet in Paraview. But PUMI_PRISM causes the runtime error, the seg fault.

#include <mpi.h>
#include <pumi.h>

#include <cassert>
#include <iostream>
int main(int argc, char** argv)
{
        MPI_Init(&argc, &argv);
        pumi_start();

        pGeom g = pumi_geom_load("null", "null");
        pMesh mesh = pumi_mesh_create(g, 3, false);

        double points[6][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1},{1,0,1},{0,1,1}};
        pMeshEnt vertices[6];

        for (int i = 0; i < 6; ++i)
                vertices[i] = pumi_mesh_createVtx(mesh, NULL, points[i]);
        pumi_mesh_createElem(mesh, NULL, PUMI_PRISM, vertices);

        pumi_mesh_freeze(mesh);
        pumi_mesh_verify(mesh);


        // ES: after mesh is created, I'd like to change its order
        pumi_mesh_setShape(mesh, pumi_shape_getLagrange(2));

        double coord[3];
        pMeshIter mit = mesh->begin(1);
        pMeshEnt e;
        while(e = mesh->iterate(mit))
        {
            pumi_node_getCoord(e, 0, coord);   // 0th node on an edge
            coord[1] += 0.1;
            pumi_node_setCoord(e, 0, coord);
        }
        mesh->end(mit);

        pumi_mesh_write(mesh, "oneprism", "vtk");
        pumi_mesh_delete(mesh);

        pumi_finalize();
        MPI_Finalize();
}

gdb shows below.

Program received signal SIGSEGV, Segmentation fault.
apf::countElementNodes (n=0x8f15c0, e=0x5) at /home/esyoon/core/apf/apfVtk.cc:411
411	  return n->getShape()->getEntityShape(n->getMesh()->getType(e))->countNodes();
(gdb) 

Could you tell me if I'm using APIs in a wrong way and how I can correct this issue, please?

Thank you in advance!

[seegyoung]

Hi Eisung, I will take a look and get back to you.

[eisungy]

Hello, Seegyoung. Thank you for taking care of this issue. 🙏

[cwsmith]

@mortezah

[mortezah]

@eisungy does pumi_mesh_write write the mesh without error before changing the mesh shape?

[mortezah]

@eisungy nevermind my previous question.

We do not have the implementation of quadratic prisms shape functions in PUMI (see here https://github.com/SCOREC/core/blob/c2996745a3e98304dcbe30a34903d329210dcc23/apf/apfShape.cc#L538), and therefore getEntityShape returns null. It would be good if you can check in the debugger that the type of entity in

Program received signal SIGSEGV, Segmentation fault.
apf::countElementNodes (n=0x8f15c0, e=0x5) at /home/esyoon/core/apf/apfVtk.cc:411
411	  return n->getShape()->getEntityShape(n->getMesh()->getType(e))->countNodes();
(gdb) 

is actually a prisim.

[mortezah]

That being said, based on the error code it seems that the only reason getEntityShape is called is to then get the number of nodes associated with the entity which is known.

To make things work for vtk write, I guess the following should work. Implement the class for prism similar to, for example, the tet here (https://github.com/SCOREC/core/blob/c2996745a3e98304dcbe30a34903d329210dcc23/apf/apfShape.cc#L481)

Implement countNodes() which should return the total number of nodes and the rest of the functions should have only the following line

fail("not implemented for PRISMS!");

[mortezah]

@eisungy @cwsmith @seegyoung I have started working on this in PR #362, and here is the latest update:

So far I have been able to fix the problem with writeVtk. That is, now the mesh is written to vtk without a segmentation file. However, the node locations in the vtk file seem incorrect. I think this has to do with how we number nodes on the quadratic prism and how vtk expect them to be and most likely those are different. Fixing this is relatively straightforward however at the moment my schedule is pretty full. I can get back to working on this when my schedule opens up a bit.

[eisungy]

Hi @mortezah , thank you for your effort and time. Also I appreciate your pointing out where I should check and read in apf. I'll take a look at your change in #362 . 👍

[eisungy]

I'm looking at the remaining issue. Just for a later reference,

cell represents a parabolic, 18-node isoparametric wedge

vtkBiQuadraticQuadraticWedge is a concrete implementation of vtkNonLinearCell to represent a three-dimensional, 18-node isoparametric biquadratic wedge. The interpolation is the standard finite element, biquadratic-quadratic isoparametric shape function plus the linear functions. The cell includes a mid-edge node. The ordering of the 18 points defining the cell is point ids (0-5,6-15, 16-18) where point ids 0-5 are the six corner vertices of the wedge; followed by nine midedge nodes (6-15) and 3 center-face nodes. Note that these midedge nodes correspond lie on the edges defined by (0,1), (1,2), (2,0), (3,4), (4,5), (5,3), (0,3), (1,4), (2,5), and the center-face nodes are laying in quads 16-(0,1,4,3), 17-(1,2,5,4) and (2,0,3,5).

from https://vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html

One more ...