[MOAB-dev] Possible bug in 2-D quadratic elements

[Wollaber, Allan B.]

To whom it may concern:

I am a developer for the neutronics part (UNIC) of the SHARP project at
ANL, and as of late I have been using MOAB to read CUBIT 10.2 files in
order to generate UNIC-native mesh objects.  The mesh data in UNIC is
organized by the MATERIAL_SET (block) tags.  Therefore, to bring in the
entities, I have a subroutine (UNIC is written in FORTRAN90) that (1)
calls "getEntitiesRec" to get all the elements in a block, and (2) calls
"getEntArrAdj" to get the vertex connectivities of these entities, among
other things.  For an example, see the code snippet below.

This works as desired for all 3D elements that I've tried and for linear
2-D elements. However, I have encountered a difficulty regarding
quadratic 2-D elements, for instance: TRI6, QUAD8, and QUAD9. The
connectivity array contains only the linear vertices. For example,
although a TRI6 has 6 nodes, the connectivity array for a single TRI6 is
3 integers long and it correctly indexes into the global vertex
coordinates with the correct MBCN ordering.  For a QUAD8, the
connectivity indexing is 4 integers long, and for a QUAD9 (strangely),
it is 5 integers long.  

I've run my code in valgrind and it does not appear to have any memory
leaks or errors with regard to my coding. Also, I've spoken to Tim about
this today, and he thinks it may be a bug.  Any further assistance would
be appreciated.

Thanks,
Allan

Dr. Allan Wollaber
Assistant Nuclear Engineer
Nuclear Engineering Division
Argonne National Laboratory
Argonne, IL 60439-4842


! Variables of interest
iMesh_Instance        ,intent(inout) :: iMesh
IMESH_HANDLE          ,intent(in) :: block
IMESH_HANDLE          ,intent(in) :: blocksHandle
IMESH_INT             ,intent(in) :: Element_Type
! elements
IMESH_HANDLE, allocatable, dimension(:) :: elems
IMESH_INT             :: elemsAllocated, elemsSize
! connectivity
IMESH_HANDLE          :: conn(0:*)
IMESH_HANDLE          :: connPointer
pointer                  (connPointer, conn)
IMESH_INT             :: connAllocated, connSize
! "Offset" information for adjacency 
MeshLibrary_Int       :: pointsInterp
IMESH_INT, allocatable, dimension(:) :: offset 
IMESH_INT             :: offsetAllocated, offsetSize

! --------------------------------------------------------
! Get the elements in this block
! --------------------------------------------------------
allocate(elems(numElements), stat=ios)
elemsAllocated = numElements
elemsSize      = numElements
call iMesh_getEntitiesRec(%val(iMesh), %val(block), %val(Element_Type),
&
                          %val(iMesh_ALL_TOPOLOGIES), %val(1),
&
                          %LOC(elems(1)), elemsAllocated, elemsSize,
err)
call iMeshAssert(err, iMesh, __LINE__, __FILE__, Output_Unit)

! --------------------------------------------------------------------
! Get the number of interpolation points per element and connectivity 
! --------------------------------------------------------------------
!! Need the number of vertices per block for connectivity allocation !!
connAllocated = 0
connPointer = 0
allocate(offset(numElements+1), stat=ios)
offsetAllocated = numElements+1
offsetSize      = numElements+1
call iMesh_getEntArrAdj(%VAL(iMesh), elems, %VAL(elemsSize),
%VAL(iBase_VERTEX), &
     connPointer, connAllocated, connSize, %LOC(offset(1)),
offsetAllocated, offsetSize, err)
call iMeshAssert(err, iMesh, __LINE__, __FILE__, Output_Unit)
pointsInterp = offset(2)-offset(1)


At this point, the array "conn" contains the connectivity only for the
linear vertices and connSize is the number of linear vertices.  This is
despite the global vertex coordinates correctly containing all the
global nodes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mcs.anl.gov/mailman/private/moab-dev/attachments/20090113/8c0e87a4/attachment.html>