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

Wollaber, Allan B. awollaber at anl.gov
Thu Jan 15 12:53:14 CST 2009


Actually, HEX27 elements also have the same issue, although HEX20
elements and TET10 elements do not.  Before the change, the HEX27
elements did not have this issue.  
Also, to clarify the code snippet below "tmpInt" has a value of 1 at the
return of the call to MBCN_SideNumberULong().  In our numbering system,
we interpret this as "2".  For QUAD8/9, TRI6, and HEX27 I find that all
of the local reference numbers are "2", which implies that the return
value is always "1".  

It may be that a new counterpart to getEntArrAdj() is required -- one
that returns only the "linear" vertices (the previous behavior that is
ITAPS API compliant) and one that returns all of them (the new behavior
in place now).  Is something like this feasible?

- Allan

> _____________________________________________ 
> From: 	Wollaber, Allan B.  
> Sent:	Thursday, January 15, 2009 11:12 AM
> To:	'Jason Kraftcheck'
> Cc:	moab-dev at mcs.anl.gov; Smith, Micheal A.; Tautges, Timothy J.
> Subject:	RE: [MOAB-dev] Possible bug in 2-D quadratic elements
> 
> 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/8923c06f/attachment.htm>


More information about the moab-dev mailing list