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

[Wollaber, Allan B.]

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
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?

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.

The SHARP group (with Tim's guidance) is primarily focused on using MOAB
for multiphysics coupling.  Is it possible to expose whatever function
that would rectify this in the MOAB API to C/FORTRAN?  It is acceptable
to me that we may lose compatibility with the iMesh API in the short
term.  

Thanks again,
- Allan  



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

Wollaber, Allan B. wrote:
> 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.
> 

Hi,

You are using the iMesh API to access MOAB data.  The ITAPS workgroup
has not yet reached a consensus on how higher-order nodes will be
handled in this API.  The current behavior, as required by the interface
standard, is to return only the topological connectivity (first-order
nodes) for adjacency queries.  So, unfortunately, this is a feature
rather than a bug (except for the 5-handle result for QUAD9, which is a
bug.)

The normal work-around for this is to use MOAB's native API rather than
the iMesh API.  However, as you are using FORTRAN and MOAB's native API
is C++, I assume that using MOAB's native API is not an option.  I have
just checked a change to MOAB that provides an alternate work-around.
Using the latest version of MOAB, if you pass the string
";HIGHER_ORDER_ADJ" (the leading ';'
is not a typo) to the iMesh_createMesh function then iMesh_getEntArr and
iMesh_getEntAdjArr will return the complete connectivity.

Please note that this will not work with other implementations of the
iMesh implementation, and may no longer be supported by MOAB's
implementation of that API once the ITAPS workgroup has defined a
standard mechanism for obtaining this information.

- jason