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

Thu Jan 15 11:12:02 CST 2009

Jason (Tim),

The connectivity indices are now the correct length for all the elements
that I've tried (HEX8/20/27, QUAD4/8/9, TRI3/6).  This is great! 

I have a new question now, however: At one point in the mesh processing
I must determine the local reference number of a facet (in a sideset
with the NEUMANN_SET tag) with respect to the element it is bounding.
To do this, I use the MBCN_SideNumber routine.  This has worked well for
every 3-D element that I've tried (hexes, tets) and for the linear 2-D
elements (tri3, quad4), but not for the quadratic 2-D elements (tri6,
quad8/9), even before the change that was just made to the repository.
For the quadratic 2-D elements, the local facet number is always
returned as "1".  The following code snippet should help illuminate how
I am doing this:

      ! --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
--  --
      ! Get the reference surface number of the side with respect to the
element
      ! For this we use the "Moab Canonical Numbering" (MBCN) routines
      ! --  --  --  --  --  --  --  --  --  --  --  --  --  --  --  --
--  --
      ! Determine the topology of the adjacent element (region or
surface?)
      call iMesh_getEntTopo(%VAL(iMesh), %VAL(adjElem(1)),
elementTopology, err)
      call iMeshAssert(err, iMesh, __LINE__, __FILE__, Output_Unit)
      ! Get MBCN type equivalent to the iMesh topology 
      call iMesh_MBCNType(%VAL(elementTopology), elementMBCNtype)
      ! Get the vertices in this side
      sideVerticesAllocated = ElementLibrary_MaxVertices  
      sideVerticesSize = 0
      call iMesh_getEntAdj(%VAL(iMesh), %VAL(sides(j)),
%VAL(iBase_VERTEX), &
           %LOC(sideVertices(1)), sideVerticesAllocated,
sideVerticesSize, err)
      call iMeshAssert(err, iMesh, __LINE__, __FILE__, Output_Unit)
      ! Get the connectivity of the vertices for this element
      adjElemConnAllocated   = ElementLibrary_MaxVertices
      adjElemOffsetAllocated = ElementLibrary_MaxVertices
      call iMesh_getEntArrAdj(%VAL(iMesh), adjElem, %VAL(adjElemSize),
&
           %VAL(iBase_VERTEX), %LOC(adjElemConn(1)),
adjElemConnAllocated, adjElemConnSize, &
           %LOC(adjElemOffset(1)), adjElemOffsetAllocated,
adjElemOffsetSize, err)
      call iMeshAssert(err, iMesh, __LINE__, __FILE__, Output_Unit)
      ! Finally, request the MBCN reference side number
      idum = 0
      call MBCN_SideNumberULong(adjElemConn,   %VAL(elementMBCNtype),
sideVertices, &
           %VAL(sideVerticesSize), %VAL(numDimensions-1), tmpInt, idum,
idum)
      ! Convert indexing from 0->1 based / possible type conversion
      sideToReferenceNo(sidesCounter) = tmpInt+1
      if (Debug_Print > 5) then
         write(Output_Unit, '("[NTmesh]...Side ", I9, " has reference
no.",7X,I4,".")') &
                            sidesCounter, tmpInt+1
      end if

Tim, since this is essentially a simple modification to your original
code, do you have any insight regarding this?  Did you ever try the
previous code for 2-D quadratic elements?

I'm not sure how the change that was just made to getEntArrAdj() should
or should not change this code.  Since the MBCN call works with the
higher-order vertices for the 3-D elements, I had hoped that the fix
that was just checked in would correct the local facet numbering as
well, but it didn't.  Any suggestions would again be appreciated.

Thanks,
- Allan

-----Original Message-----
From: Jason Kraftcheck [mailto:kraftche at cae.wisc.edu] 
Sent: Wednesday, January 14, 2009 5:59 PM
To: Wollaber, Allan B.
Cc: moab-dev at mcs.anl.gov; Smith, Micheal A.; Tautges, Timothy J.
Subject: Re: [MOAB-dev] Possible bug in 2-D quadratic elements

The current version of MOAB now has the following relevant changes:

o Remove "HIGHER_ORDER_ADJ" option for iMesh_newMesh o Restore previous
iMesh_getEntArrAdj behavior: return full
   connectivity list for higher-order elements o Fix bug reading in
higher-order face elements from CUB files.
o Fix incorrect connectivity order for tet8, tet14, and hex27
   elements read from CUB files.

The last item above may be an issue if you are working with those
element types (volume elements with mid-face nodes.)  MOAB's node
ordering is not the same as PATRAN or ExodusII for many higher-order
element types.  The CUB file reader was previously storing hex27
connectivity in the PATRAN rather than MOAB order.  Further, it was
storing tet8 and tet14 connectivity in an order that was neither
PATRAN/Exodus nor MOAB (some internal Cubit ordering
apparently.)

- jason

Wollaber, Allan B. wrote:
> Jason,
> 
> I had forgotten that I called both "iMesh_newMesh" and "iMesh_load".
> Thanks for clearing that up.
> 
> I think we can rule out that someone has modified my local copy of 
> MOAB, since I installed it locally (in my home directory) using 
> autoreconf around mid-December.
> 
> Let me know what you find out, or if I can do anything useful while 
> you figure out what's going on from your side.
> 
> Thanks,
> - Allan
> 
> 
> 
> -----Original Message-----
> From: Jason Kraftcheck [mailto:kraftche at cae.wisc.edu]
> Sent: Tuesday, January 13, 2009 4:11 PM
> To: Wollaber, Allan B.
> Cc: moab-dev at mcs.anl.gov; Smith, Micheal A.; Tautges, Timothy J.
> Subject: Re: [MOAB-dev] Possible bug in 2-D quadratic elements
> 
> Wollaber, Allan B. wrote:
>> Jason,
>>
>> Thank you for your prompt reply. 
>>
>> Actually, I never call "iMesh_createMesh", and I can't find this 
>> function in any of the header files.  I'm working out of the trunk of

>> MOAB (grep -ir createmesh *).  Do you possibly mean "iMesh_newMesh"?
>> I
> 
> Yes, I meant the latter.  Apologies for any confusion.
> 
>> would prefer to just use "iMesh_load", since I am just loading the 
>> mesh directly from the CUBIT file.  Would it be necessary to create a

>> new mesh instance and copy in the mesh loaded from CUBIT to achieve 
>> this? If so, I'm actually not sure how best to do that.  Is there a 
>> "copy constructor" somewhere?
> 
> You need to call iMesh_newMesh to create the "iMesh_instance" object 
> which is later passed to iMesh_load.  The name "newMesh" is a bit 
> misleading, as it doesn't create a mesh (new or otherwise), but rather

> an instance of the database.
> 
>> Also, will using this setting affect the success I've had thus far 
>> with the quadratic 3-D elements? It is only the 2-D elements that 
>> I've
> 
>> had trouble with.
>>
> 
> It shouldn't.  But it appears that the problem is more complicated 
> that I had at first thought.  Apparently the version of MOAB you are 
> using is always returning the higher-order nodes from getEntArrAdj 
> (regardless of any
> setting.)  I don't know if the MOAB behavior has changed or if someone

> has modified your copy of MOAB.  So the problem you are seeing is a 
> bug rather than an API issue.  I'm looking into this right now.  So 
> far the problem appears to be that the 'cub' file reader in MOAB 
> looses the higher-order nodes on the edges of face elements (thus 5 
> nodes for a
> quad9.)
> 
> 
> - jason
> 

--
"A foolish consistency is the hobgoblin of little minds" - Ralph Waldo
Emerson

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mcs.anl.gov/mailman/private/moab-dev/attachments/20090115/031e4fc1/attachment.htm>