[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jun 6 11:24:27 CDT 2013
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/adea7c874b79/
Changeset: adea7c874b79
Branch: None
User: jhu
Date: 2013-06-06 18:21:13
Summary: Added user's guide document to doxygen.
Affected #: 1 file
diff --git a/src/moab/moabUG.h b/src/moab/moabUG.h
new file mode 100644
index 0000000..cb244a8
--- /dev/null
+++ b/src/moab/moabUG.h
@@ -0,0 +1,1062 @@
+/*! \page the User's Guide (MOAB 4.6)
+=20
+ \subpage team=20
+=20
+ \subpage contents
+=20
+ \subpage figures
+=20
+ \subpage tables
+=20
+ \subpage musthave
+
+ \subpage building
+
+ \page team MOAB team members
+ <h2>The MOAB Team, including: </h2>
+=20
+ - Timothy J. Tautges (Argonne National Lab, Univ Wisconsin-Madison)=20
+ - Iulian Grindeanu (Argonne National Lab)=20
+ - Rajeev Jain (Argonne National Lab)
+ - Xiabing Xu (Argonne National Lab)
+
+
+ <h2>Emeritus members:</h2>
+=20
+ - Jason A. Kraftcheck
+ - Brandon M. Smith
+ - Hong-Jun Kim
+ - Jim Porter
+=20
+ \page contents Table of Contents
+=20
+ \ref introduction =20
+
+ \ref interface =20
+
+ \ref twoone =20
+
+ \ref twotwo =20
+
+ \ref twothree =20
+
+ \ref twofour =20
+
+ \ref api =20
+
+ \ref services =20
+
+ \ref fourone =20
+
+ \ref fourtwo =20
+
+ \ref fourthree =20
+
+ \ref fourfour =20
+
+ \ref fourfive =20
+
+ \ref foursix =20
+
+ \ref parallel =20
+
+ \ref fiveone =20
+
+ \ref fivetwo =20
+
+ \ref fivethree =20
+
+ \ref fivefour =20
+
+ \ref applications =20
+
+ \ref implementation =20
+
+ \ref representation =20
+
+ \ref element =20
+
+ \ref nineone =20
+
+ \ref ninetwo =20
+
+ \ref ninethree =20
+
+ \ref performance =20
+
+ \ref conclusions =20
+
+ \ref references=20
+
+ \section introduction 1. Introduction
+
+In scientific computing, systems of partial differential equations (PDEs) =
are solved on computers. One of the most widely used methods to solve PDEs=
numerically is to solve over discrete neighborhoods or =E2=80=9Celements=
=E2=80=9D of the domain. Popular discretization methods include Finite Dif=
ference (FD), Finite Element (FE), and Finite Volume (FV). These methods r=
equire the decomposition of the domain into a discretized representation, w=
hich is referred to as a =E2=80=9Cmesh=E2=80=9D. The mesh is one of the fu=
ndamental types of data linking the various tools in the analysis process (=
mesh generation, analysis, visualization, etc.). Thus, the representation =
of mesh data and operations on those data play a very important role in PDE=
-based simulations.
+=20
+MOAB is a component for representing and evaluating mesh data. MOAB can s=
tore structured and unstructured mesh, consisting of elements in the finite=
element =E2=80=9Czoo=E2=80=9D, along with polygons and polyhedra. The fun=
ctional interface to MOAB is simple, consisting of only four fundamental da=
ta types. This data is quite powerful, allowing the representation of most=
types of metadata commonly found on the mesh. Internally MOAB uses array-=
based storage for fine-grained data, which in many cases provides more effi=
cient access, especially for large portions of mesh and associated data. M=
OAB is optimized for efficiency in space and time, based on access to mesh =
in chunks rather than through individual entities, while also versatile eno=
ugh to support individual entity access.
+
+The MOAB data model consists of the following four fundamental types: mesh=
interface instance, mesh entities (vertex, edge, tri, etc.), sets, and tag=
s. Entities are addressed through handles rather than pointers, to allow t=
he underlying representation of an entity to change without changing the ha=
ndle to that entity. Sets are arbitrary groupings of mesh entities and oth=
er sets. Sets also support parent/child relationships as a relation distin=
ct from sets containing other sets. The directed graph provided by set par=
ent/child relationships is useful for embedding graphs whose nodes include =
collections of mesh entities; this approach has been used to represent a wi=
de variety of application-specific data, including geometric model topology=
, processor partitions, and various types of search trees. Tags are named =
data which can be assigned to the mesh as a whole, individual entities, or =
sets. Tags are a mechanism for attaching data to individual entities, and =
sets are a mechanism for describing relations between entities; the combina=
tion of these two mechanisms is a powerful yet simple interface for represe=
nting metadata or application-specific data.
+
+Various mesh-related tools are provided with MOAB or can be used directly =
with MOAB. These tools can be used for mesh format translation (mbconvert)=
, mesh skinning (Skinner class), solution transfer between meshes (MBCouple=
r tool), ray tracing and other geometric searches (OrientedBoxTreeTool, Ada=
ptiveKDTree), visualization (vtkMOABReader tool), and relation between mesh=
and geometric models (the separately-packed Lasso tool). These tools are =
described later in this document.
+
+MOAB is written in the C++ programming language, with applications interac=
ting with MOAB mostly through its moab::Interface class. All of the MOAB f=
unctions and classes are isolated in and accessed through the moab namespac=
e<sup>1</sup>. The remainder of this report gives class and function names =
without the =E2=80=9Cmoab::=E2=80=9D namespace qualification; unless otherw=
ise noted, the namespace qualifier should be added to all class and functio=
n names referenced here. MOAB also implements the iMesh interface, which i=
s specified in C but can be called directly from other languages. Almost a=
ll of the functionality in MOAB can be accessed through the iMesh interface=
. MOAB is developed and supported on the Linux and MacOS operating systems=
, as well as various HPC operating systems. MOAB can be used on parallel c=
omputing systems as well, including both clusters and high-end parallel sys=
tems like IBM BG/P and Cray systems. MOAB is released under a standard LGP=
L open source software license.
+
+MOAB is used in several ways in various applications. MOAB serves as the =
underlying mesh data representation in several scientific computing applica=
tions [1]. MOAB can also be used as a mesh format translator, using reader=
s and writers included in MOAB. MOAB has also been used as a bridge to cou=
ple results in multi-physics analysis and to link these applications with o=
ther mesh services [2].
+
+The remainder of this report is organized as follows. Section 2, =E2=80=
=9CGetting Started=E2=80=9D, provides a few simple examples of using MOAB t=
o perform simple tasks on a mesh. Section 3 discusses the MOAB data model =
in more detail, including some aspects of the implementation. Section 4 su=
mmarizes the MOAB function API. Section 5 describes some of the tools incl=
uded with MOAB, and the implementation of mesh readers/writers for MOAB. S=
ection 6 describes how to build MOAB-based applications. Section 7 contain=
s a brief description of MOAB=E2=80=99s relation to the iMesh mesh interfac=
e. Sections 8 and 9 discuss MOAB's representations of structured and spect=
ral element meshes, respectively. Section 10 gives helpful hints for acces=
sing MOAB in an efficient manner from applications. Section 11 gives a con=
clusion and future plans for MOAB development. Section 12 gives references=
cited in this report.
+
+Several other sources of information about MOAB may also be of interest to=
readers. Meta-data conventions define how sets and /or tags are used toge=
ther to represent various commonly-used simulation constructs; conventions =
used by MOAB are described in Ref [4], which is also included in the MOAB s=
ource distribution. This document is maintained separately from this docum=
ent, since it is expected to change over time. The MOAB project maintains =
a wiki [5], which links to most MOAB-related information. MOAB also uses s=
everal mailing lists [6],[7] for MOAB-related discussions. Potential users=
are encouraged to interact with the MOAB team using these mailing lists.
+
+<sup>1</sup> Non-namespaced names are also provided for backward compatibi=
lity, with the =E2=80=9CMB=E2=80=9D prefix added to the class or variable n=
ame.
+
+ \ref contents
+
+ \section interface 2. MOAB Data Model
+The MOAB data model describes the basic types used in MOAB and the languag=
e used to communicate that data to applications. This chapter describes th=
at data model, along with some of the reasons for some of the design choice=
s in MOAB.
+
+ \ref contents
+
+ \subsection twoone 2.1. MOAB Interface
+MOAB is written in C++. The primary interface with applications is throug=
h member functions of the abstract base class Interface. The MOAB library =
is created by instantiating Core, which implements the Interface API. Mult=
iple instances of MOAB can exist concurrently in the same application; mesh=
entities are not shared between these instancesi<sup>1</sup>. MOAB is mos=
t easily viewed as a database of mesh objects accessed through the instance=
. No other assumptions explicitly made about the nature of the mesh stored=
there; for example, there is no fundamental requirement that elements fill=
space or do not overlap each other geometrically.
+=20
+<sup>1</sup> One exception to this statement is when the parallel interfac=
e to MOAB is used; in this case, entity sharing between instances is handle=
d explicitly using message passing. This is described in more detail in Se=
ction 5 of this document.
+
+ \ref contents
+
+ \subsection twotwo 2.2. Mesh Entities
+MOAB represents the following topological mesh entities: vertex, edge, tri=
angle, quadrilateral, polygon, tetrahedron, pyramid, prism, knife, hexahedr=
on, polyhedron. MOAB uses the EntityType enumeration to refer to these ent=
ity types (see Table 1). This enumeration has several special characterist=
ics, chosen intentionally: the types begin with vertex, entity types are gr=
ouped by topological dimension, with lower-dimensional entities appearing b=
efore higher dimensions; the enumeration includes an entity type for sets (=
described in the next section); and MBMAXTYPE is included at the end of thi=
s enumeration, and can be used to terminate loops over type. In addition t=
o these defined values, the an increment operator (++) is defined such that=
variables of type EntityType can be used as iterators in loops.
+MOAB refers to entities using =E2=80=9Chandles=E2=80=9D. Handles are impl=
emented as long integer data types, with the four highest-order bits used t=
o store the entity type (mesh vertex, edge, tri, etc.) and the remaining bi=
ts storing the entity id. This scheme is convenient for applications becau=
se:
+- Handles sort lexicographically by type and dimension; this can be useful=
for grouping and iterating over entities by type.
+- The type of an entity is indicated by the handle itself, without needing=
to call a function.
+- Entities allocated in sequence will typically have contiguous handles; t=
his characteristic can be used to efficiently store and operate on large li=
sts of handles.
+.
+
+This handle implementation is exposed to applications intentionally, becau=
se of optimizations that it enables, and is unlikely to change in future ve=
rsions.
+
+ \subsection tableone Table 1: Values defined for the EntityType enumerat=
ed type.
+<table border=3D"1">
+<tr>
+<td>MBVERTEX =3D 0</td>
+<td>MBPRISM</td>
+</tr>
+<tr>
+<td>MBEDGE</td>
+<td>MBKNIFE</td>
+</tr>
+<tr>
+<td>MBTRI</td>
+<td>MBHEX</td>
+</tr>
+<tr>
+<td>MBQUAD</td>
+<td>MBPOLYHEDRON</td>
+</tr>
+<tr>
+<td>MBPOLYGON</td>
+<td>MBENTITYSET</td>
+</tr>
+<tr>
+<td>MBTET</td>
+<td>MBMAXTYPE</td>
+</tr>
+<tr>
+<td>MBPYRAMID</td>
+<td></td>
+</tr>
+</table>
+
+MOAB defines a special class for storing lists of entity handles, named Ra=
nge. This class stores handles as a series of (start_handle, end_handle) s=
ubrange tuples. If a list of handles has large contiguous ranges, it can b=
e represented in almost constant size using Range. Since entities are typi=
cally created in groups, e.g. during mesh generation or file import, a high=
degree of contiguity in handle space is typical. Range provides an interf=
ace similar to C++ STL containers like std::vector, containing iterator dat=
a types and functions for initializing and iterating over entity handles st=
ored in the range. Range also provides functions for efficient Boolean ope=
rations like subtraction and intersection. Most API functions in MOAB come=
in both range-based and vector-based variants. By definition, a list of e=
ntities stored in an Range is always sorted, and can contain a given entity=
handle only once. Range cannot store the handle 0 (zero).
+
+Typical usage of an Range object would look like:
+
+\code
+using namespace moab;
+ int my_function(Range &from_range) {
+ int num_in_range =3D from_range.size();
+ Range to_range;
+ Range::iterator rit;
+ for (rit =3D from_range.begin(); rit !=3D from_range.end(); rit++) {
+ EntityHandle this_ent =3D *rit;
+ to_range.insert(this_ent);
+ }
+ }
+\endcode
+
+Here, the range is iterated similar to how std::vector is iterated.
+
+ \subsection adjacencies 2.2.1. Adjacencies & AEntities=20
+
+The term adjacencies is used to refer to those entities topologically conn=
ected to a given entity, e.g. the faces bounded by a given edge or the vert=
ices bounding a given region. The same term is used for both higher-dimens=
ional (or bounded) and lower-dimensional (or bounding) adjacent entities. =
MOAB provides functions for querying adjacent entities by target dimension,=
using the same functions for higher- and lower-dimension adjacencies. By =
default, MOAB stores the minimum data necessary to recover adjacencies betw=
een entities. When a mesh is initially loaded into MOAB, only entity-verte=
x (i.e. =E2=80=9Cdownward=E2=80=9D) adjacencies are stored, in the form of =
entity connectivity. When =E2=80=9Cupward=E2=80=9D adjacencies are request=
ed for the first time, e.g. from vertices to regions, MOAB stores all verte=
x-entity adjacencies explicitly, for all entities in the mesh. Non-vertex =
entity to entity adjacencies are never stored, unless explicitly requested =
by the application.
+
+In its most fundamental form, a mesh need only be represented by its verti=
ces and the entities of maximal topological dimension. For example, a hexa=
hedral mesh can be represented as the connectivity of the hex elements and =
the vertices forming the hexes. Edges and faces in a 3D mesh need not be e=
xplicitly represented. We refer to such entities as =E2=80=9CAEntities=E2=
=80=9D, where =E2=80=98A=E2=80=99 refers to =E2=80=9CAuxiliary=E2=80=9D, =
=E2=80=9CAncillary=E2=80=9D, and a number of other words mostly beginning w=
ith =E2=80=98A=E2=80=99. Individual AEntities are created only when reques=
ted by applications, either using mesh modification functions or by request=
ing adjacencies with a special =E2=80=9Ccreate if missing=E2=80=9D flag pas=
sed as =E2=80=9Ctrue=E2=80=9D. This reduces the overall memory usage when =
representing large meshes. Note entities must be explicitly represented be=
fore they can be assigned tag values or added to entity sets (described in =
following Sections).
+
+\ref contents
+
+ \subsection twothree 2.3. Entity Sets
+Entity sets are used to store arbitrary collections of entities and other =
sets. Sets are used for a variety of things in mesh-based applications, fr=
om the set of entities discretizing a given geometric model entity to the e=
ntities partitioned to a specific processor in a parallel finite element ap=
plication. MOAB entity sets can also store parent/child relations with oth=
er entity sets, with these relations distinct from contains relations. Par=
ent/child relations are useful for building directed graphs with graph node=
s representing collections of mesh entities; this construct can be used, fo=
r example, to represent an interface of mesh faces shared by two distinct c=
ollections of mesh regions. MOAB also defines one special set, the =E2=80=
=9Croot set=E2=80=9D or the interface itself; all entities are part of this=
set by definition. Defining a root set allows the use of a single set of =
MOAB API functions to query entities in the overall mesh as well as its sub=
sets.
+
+MOAB entity sets can be one of two distinct types: list-type entity sets p=
reserve the order in which entities are added to the set, and can store a g=
iven entity handle multiple times in the same set; set-type sets are always=
ordered by handle, regardless of the order of addition to the set, and can=
store a given entity handle only once. This characteristic is assigned wh=
en the set is created, and cannot be changed during the set=E2=80=99s lifet=
ime.
+
+MOAB provides the option to track or not track entities in a set. When en=
tities (and sets) are deleted by other operations in MOAB, they will also b=
e removed from containing sets for which tracking has been enabled. This b=
ehavior is assigned when the set is created, and cannot be changed during t=
he set=E2=80=99s lifetime. The cost of turning tracking on for a given set=
is sizeof(EntityHandle) for each entity added to the set; MOAB stores cont=
aining sets in the same list which stores adjacencies to other entities.
+
+Using an entity set looks like the following:
+\code
+using namespace moab;
+// load a file using MOAB, putting the loaded mesh into a file set
+EntityHandle file_set;
+ErrorCode rval =3D moab->create_meshset(MESHSET_SET, file_set);
+rval =3D moab->load_file(=E2=80=9Cfname.vtk=E2=80=9D, &file_set);
+Range set_ents;
+// get all the 3D entities in the set
+rval =3D moab->get_entities_by_dimension(file_set, 3, set_ents);
+\endcode
+
+Entity sets are often used in conjunction with tags (described in the next=
section), and provide a powerful mechanism to store a variety of meta-data=
with meshes.
+
+\ref contents
+
+ \subsection twofour 2.4. Tags=20
+
+Applications of a mesh database often need to attach data to mesh entities=
. The types of attached data are often not known at compile time, and can =
vary across individual entities and entity types. MOAB refers to this atta=
ched data as a =E2=80=9Ctag=E2=80=9D. Tags can be thought of loosely as a =
variable, which can be given a distinct value for individual entities, enti=
ty sets, or for the interface itself. A tag is referenced using a handle, =
similarly to how entities are referenced in MOAB. Each MOAB tag has the fo=
llowing characteristics, which can be queried through the MOAB interface:
+- Name
+- Size (in bytes)
+- Storage type
+- Data type (integer, double, opaque, entity handle)
+- Handle
+.
+
+The storage type determines how tag values are stored on entities. =20
+
+- Dense: Dense tag values are stored in arrays which match arrays of conti=
guous entity handles. Dense tags are more efficient in both storage and me=
mory if large numbers of entities are assigned the same tag. Storage for a=
given dense tag is not allocated until a tag value is set on an entity; me=
mory for a given dense tag is allocated for all entities in a given sequenc=
e at the same time.
+- Sparse: Sparse tags are stored as a list of (entity handle, tag value) t=
uples, one list per sparse tag, sorted by entity handle.
+- Bit: Bit tags are stored similarly to dense tags, but with special handl=
ing to allow allocation in bit-size amounts per entity.
+.
+
+MOAB also supports variable-length tags, which can have a different length=
for each entity they are assigned to. Variable length tags are stored sim=
ilarly to sparse tags.
+
+The data type of a tag can either be one understood at compile time (integ=
er, double, entity handle), in which case the tag value can be saved and re=
stored properly to/from files and between computers of different architectu=
re (MOAB provides a native HDF5-based save/restore format for this purpose;=
see Section 4.6). The opaque data type is used for character strings, or =
for allocating =E2=80=9Craw memory=E2=80=9D for use by applications (e.g. f=
or storage application-defined structures or other abstract data types). T=
hese tags are saved and restored as raw memory, with no special handling fo=
r endian or precision differences.
+
+An application would use the following code to attach a double-precision t=
ag to vertices in a mesh, e.g. to assign a temperature field to those verti=
ces:
+
+\code
+using namespace moab;
+// load a file using MOAB and get the vertices
+ErrorCode rval =3D moab->load_file(=E2=80=9Cfname.vtk=E2=80=9D);
+Range verts;
+rval =3D moab->get_entities_by_dimension(0, 0, verts);
+// create a tag called =E2=80=9CTEMPERATURE=E2=80=9D
+Tag temperature;
+double def_val =3D -1.0d-300, new_val =3D 273.0;
+rval =3D moab->tag_create(=E2=80=9CTEMPERATURE=E2=80=9D, sizeof(double), M=
B_TAG_DENSE,=20
+ MB_TYPE_DOUBLE, temperature, &def_val);
+// assign a value to vertices
+for (Range::iterator vit =3D verts.begin();=20
+ vit !=3D verts.end(); vit++)=20
+ rval =3D moab->tag_set_data(temperature, &(*rit), 1, &new_val);
+
+\endcode
+
+The semantic meaning of a tag is determined by applications using it. How=
ever, to promote interoperability between applications, there are a number =
of tag names reserved by MOAB which are intended to be used by convention. =
Mesh readers and writers in MOAB use these tag conventions, and applicatio=
ns can use them as well to access the same data. Ref. [4] maintains an up-t=
o-date list of conventions for meta-data usage in MOAB.
+
+ \ref contents
+
+ \section api 3. MOAB API Design Philosophy and Summary
+
+This section describes the design philosophy behind MOAB, and summarizes t=
he functions, data types and enumerated variables in the MOAB API. A compl=
ete description of the MOAB API is available in online documentation in the=
MOAB distribution [8].
+
+MOAB is designed to operate efficiently on collections of entities. Entit=
ies are often created or referenced in groups (e.g. the mesh faces discreti=
zing a given geometric face, the 3D elements read from a file), with those =
groups having some form of temporal or spatial locality. The interface pro=
vides special mechanisms for reading data directly into the native storage =
used in MOAB, and for writing large collections of entities directly from t=
hat storage, to avoid data copies. MOAB applications structured to take ad=
vantage of that locality will typically operate more efficiently.
+
+MOAB has been designed to maximize the flexibility of mesh data which can =
be represented. There is no explicit constraint on the geometric structure=
of meshes represented in MOAB, or on the connectivity between elements. I=
n particular, MOAB allows the representation of multiple entities with the =
same exact connectivity; however, in these cases, explicit adjacencies must=
be used to distinguish adjacencies with AEntities bounding such entities.
+
+The number of vertices used to represent a given topological entity can va=
ry, depending on analysis needs; this is often the case in FEA. For exampl=
e, applications often use =E2=80=9Cquadratic=E2=80=9D or 10-vertex tetrahed=
ral, with vertices at edge midpoints as well as corners. MOAB does not dis=
tinguish these variants by entity type, referring to all variants as =E2=80=
=9Ctetrahedra=E2=80=9D. The number of vertices for a given entity is used =
to distinguish the variants, with canonical numbering conventions used to d=
etermine placement of the vertices [9]. This is similar to how such variat=
ions are represented in the Exodus [10] and Patran [11] file formats. In p=
ractice, we find that this simplifies coding in applications, since in many=
cases the handling of entities depends only on the number of corner vertic=
es in the element. Some MOAB API functions provide a flag which determines=
whether corner or all vertices are requested.
+
+The MOAB API is designed to balance complexity and ease of use. This bala=
nce is evident in the following general design characteristics:
+
+- Entity lists: Lists of entities are passed to and from MOAB in a variety=
of forms. Lists output from MOAB are passed as either STL vector or Range=
data types. Either of these constructs may be more efficient in both time=
and memory, depending on the semantics of the data being requested. Input=
lists are passed as either Range=E2=80=99s, or as a pointer to EntityHandl=
e and a size. The latter allows the same function to be used when passing =
individual entities, without requiring construction of an otherwise unneede=
d STL vector.
+- Entity sets: Most query functions accept an entity set as input. Applic=
ations can pass zero to indicate a request for the whole interface. Note t=
hat this convention applies only to query functions; attempts to add or sub=
tract entities to/from the interface using set-based modification functions=
, or to add parents or children to the interface set, will fail. Allowing =
specification of the interface set in this manner avoids the need for a sep=
arate set of API functions to query the database as a whole.
+- Implicit Booleans in output lists: A number of query functions in MOAB a=
llow specification of a Boolean operation (Interface::INTERSECT or Interfac=
e::UNION). This operation is applied to the results of the query, often el=
iminating the need for code the application would need to otherwise impleme=
nt. For example, to find the set of vertices shared by a collection of qua=
drilaterals, the application would pass that list of quadrilaterals to a re=
quest for vertex adjacencies, with Interface::INTERSECT passed for the Bool=
ean flag. The list of vertices returned would be the same as if the applic=
ation called that function for each individual entity, and computed the int=
ersection of the results over all the quadrilaterals. Applications may als=
o input non-empty lists to store the results, in which case the intersectio=
n is also performed with entities already in the list. In many cases, this=
allows optimizations in both time and memory inside the MOAB implementatio=
n.=20
+.
+
+Since these objectives are at odds with each other, tradeoffs had to be ma=
de between them. Some specific issues that came up are:
+
+- Using ranges: Where possible, entities can be referenced using either ra=
nges (which allow efficient storage of long lists) or STL vectors (which al=
low list order to be preserved), in both input and output arguments.
+- Entities in sets: Accessing the entities in a set is done using the same=
functions which access entities in the entire mesh. The whole mesh is ref=
erenced by specifying a set handle of zero<sup>1</sup>.
+- Entity vectors on input: Functions which could normally take a single en=
tity as input are specified to take a vector of handles instead. Single en=
tities are specified by taking the address of that entity handle and specif=
ying a list length of one. This minimizes the number of functions, while p=
reserving the ability to input single entities.<sup>2</sup>
+.
+
+Table 2 lists basic data types and enumerated variables defined and used b=
y MOAB. Values of the ErrorCode enumeration are returned from most MOAB fu=
nctions, and can be compared to those listed in Appendix [ref-appendix].
+
+MOAB uses several pre-defined tag names to define data commonly found in v=
arious mesh-based analyses. Ref. [4] describes these meta-data conventions=
in more detail. These conventions will be added to as new conventions eme=
rge for using sets and tags in MOAB applications.
+
+ \subsection tabletwo Table 2: Basic data types and enums defined in MOAB.
+
+<table border=3D"1">
+<tr>
+<th>Enum / Type</th>
+<th>Description</th>
+</tr>
+<tr>
+<td>ErrorCode</td>
+<td>Specific error codes returned from MOAB</td>
+</tr>
+<tr>
+<td>EntityHandle</td>
+<td>Type used to represent entity handles</td>
+</tr>
+<tr>
+<td>Tag</td>
+<td>Type used to represent tag handles</td>
+</tr>
+<tr>
+<td>TagType</td>
+<td>Type used to represent tag storage type</td>
+</tr>
+<tr>
+<td>DataType</td>
+<td>Type used to represent tag data type</td>
+</tr>
+</table>
+
+Table 3 lists the various groups of functions that comprise the MOAB API. =
This is listed here strictly as a reference to the various types of functi=
onality supported by MOAB; for a more detailed description of the scope and=
syntax of the MOAB API, see the online documentation [8].
+
+ \subsection tablethree Table 3: Groups of functions in MOAB API. See Re=
f. [8] for more details.
+
+<table border=3D"1">
+<tr>
+<th>Function group</th>
+<th>Examples</th>
+<th>Description</th>
+</tr>
+<tr>
+<td>Constructor, destructor, interface</td>
+<td>Interface, ~Core, query_interface</td>
+<td>Construct/destroy interface; get pointer to read/write interface</td>
+</tr>
+<tr>
+<td>Entity query</td>
+<td>get_entities_by_dimension, get_entities_by_handle</td>
+<td>Get entities by dimension, type, etc.</td>
+</tr>
+<tr>
+<td>Adjacencies</td>
+<td>get_adjacencies, set_adjacencies, add_adjacencies</td>
+<td>Get topologically adjacent entities; set or add explicit adjacencies</=
td>
+</tr>
+<tr>
+<td>Vertex coordinates</td>
+<td>get_coords, set_coords</td>
+<td>Get/set vertex coordinates</td>
+</tr>
+<tr>
+<td>Connectivity</td>
+<td>get_connectivity, set_connectivity</td>
+<td>Get/set connectivity of non-vertex entities</td>
+</tr>
+<tr>
+<td>Sets</td>
+<td>create_meshset, add_entities, add_parent_child</td>
+<td>Create and work with entity sets</td>
+</tr>
+<tr>
+<td>Tags</td>
+<td>tag_get_data, tag_create</td>
+<td>Create, read, write tag data</td>
+</tr>
+<tr>
+<td>Handles</td>
+<td>type_from_handle, id_from_handle</td>
+<td>Go between handles and types/ids</td>
+</tr>
+<tr>
+<td>File handling</td>
+<td>load_mesh, save_mesh</td>
+<td>Read/write mesh files</td>
+</tr>
+<tr>
+<td>Geometric dimension</td>
+<td>get_dimension, set_dimension</td>
+<td>Get/set geometric dimension of mesh</td>
+</tr>
+<tr>
+<td>Mesh modification</td>
+<td>create_vertex, delete_entity</td>
+<td>Create or delete mesh entities</td>
+</tr>
+<tr>
+<td>Information</td>
+<td>list_entities, get_last_error</td>
+<td>Get or print certain information</td>
+</tr>
+<tr>
+<td>High-order nodes</td>
+<td>high_order_node</td>
+<td>Get information on high-order nodes</td>
+</tr>
+<tr>
+<td>Canonical numbering</td>
+<td>side_number</td>
+<td>Get canonical numbering information</td>
+</tr>
+</table>
+
+<sup>1</sup>In iMesh, the whole mesh is specified by a special entity set =
handle, referred to as the =E2=80=9Croot set=E2=80=9D.
+
+<sup>2</sup>Note that STL vectors of entity handles can be input in this m=
anner by using &vector[0] and vector.size() for the 1d vector address and s=
ize, respectively.
+
+ \ref contents
+
+ \section services 4. Related Mesh Services
+
+A number of mesh-based services are often used in conjunction with a mesh =
library. For example, parallel applications often need to visualize the me=
sh and associated data. Other services, like spatial interpolation or find=
ing the faces on the =E2=80=9Cskin=E2=80=9D of a 3D mesh, can be implemente=
d more efficiently using knowledge of specific data structures in MOAB. Se=
veral of these services provided with MOAB are described in this chapter.
+
+ \ref contents
+
+ \subsection fourone 4.1. Visualization
+
+Visualization is one of the most common needs associated with meshes. The=
primary tool used to visualize MOAB meshes is VisIt [12]. Users can speci=
fy that VisIt read mesh directly out of the MOAB instance, by specifying th=
e ITAPS-MOABC mesh format and a file readable by MOAB (see xxx).
+
+There are some initial capabilities in VisIt for limited viewing and manip=
ulation of tag data and some types of entity sets. Tag data is visualized =
using the same mechanisms used to view other field data in VisIt, e.g. usin=
g a pseudocolor plot; sets are viewed using VisIt=E2=80=99s SIL window, acc=
essed by selecting the SIL icon in the data selection window. xxx shows a =
vertex-based radiation temperature field computed by the Cooper rad-hydro c=
ode [1] for a subset of geometric volumes in a mesh. =20
+
+Reorganization of VisIt=E2=80=99s set handling is also underway, to increa=
se versatility and flexibility of this important mechanism.
+
+ \ref contents
+
+ \subsection fourtwo 4.2. Parallel Decomposition
+
+To support parallel simulation, applications often need to partition a mes=
h into parts, designed to balance the load and minimize communication betwe=
en sets. MOAB includes the MBZoltan tool for this purpose, constructed on =
the well-known Zoltan partitioning library [13]. After computing the parti=
tion using Zoltan, MBZoltan stores the partition as either tags on individu=
al entities in the partition, or as tagged sets, one set per part. Since a=
partition often exhibits locality similar to how the entities were created=
, storing it as sets (based on Range=E2=80=99s) is often more memory-effici=
ent than an entity tag-based representation. Xxx shows a partition compute=
d with MBZoltan (and visualized in VisIt).=20
+
+ \ref contents
+
+ \subsection fourthree 4.3. Skinner
+
+An operation commonly applied to mesh is to compute the outermost =E2=80=
=9Cskin=E2=80=9D bounding a contiguous block of elements. This skin consis=
ts of elements of one fewer topological dimension, arranged in one or more =
topological balls on the boundary of the elements. The Skinner tool comput=
es the skin of a mesh in a memory-efficient manner. Skinner uses knowledge=
about whether vertex-entity adjacencies and AEntities exist to minimize me=
mory requirements and searching time required during the skinning process. =
This skin can be provided as a single collection of entities, or as sets o=
f entities distinguished by forward and reverse orientation with respect to=
higher-dimensional entities in the set being skinned.
+
+The following code fragment shows how Skinner can be used to compute the s=
kin of a range of hex elements:
+
+ \code
+using namespace moab;
+Range hexes, faces;
+ErrorCode rval =3D moab->get_entities_by_dimension(0, 3, hexes);
+Skinner myskinner(moab);
+bool verts_too =3D false;
+ErrorCode rval =3D myskinner.find_skin(hexes, verts_too, faces);
+\endcode
+
+Skinner can also skin a mesh based on geometric topology groupings importe=
d with the mesh. The geometric topology groupings contain information abou=
t the mesh =E2=80=9Cowned=E2=80=9D by each of the entities in the geometric=
model, e.g. the model vertices, edges, etc. Links between the mesh sets c=
orresponding to those entities can be inferred directly from the mesh. Ski=
nning a mesh this way will typically be much faster than doing so on the ac=
tual mesh elements, because there is no need to create and destroy interior=
faces on the mesh.
+
+ \ref contents
+
+ \subsection fourfour 4.4. Tree Decompositions
+
+MOAB provides several mechanisms for spatial decomposition and searching i=
n a mesh:
+
+- AdaptiveKDTree: Adaptive KD tree, a space-filling decomposition with axi=
s-aligned splitting planes, enabling fast searching.
+- BSPTree: Binary Space Partition tree, with non-axis-aligned partitions, =
for fast spatial searches with slightly better memory efficiency than KD tr=
ees.
+- OrientedBoxTreeTool: Oriented Bounding Box tree hierarchy, useful for fa=
st ray-tracing on collections of mesh facets.
+.
+
+These trees have various space and time searching efficiencies. All are i=
mplemented based on entity sets and parent/child relations between those se=
ts, allowing storage of a tree decomposition using MOAB=E2=80=99s native fi=
le storage mechanism (see Section 4.6.1). MOAB=E2=80=99s entity set implem=
entation is specialized for memory efficiency when representing binary tree=
s. Tree decompositions in MOAB have been used to implement fast ray tracin=
g to support radiation transport [14], solution coupling between meshes [2]=
, and embedded boundary mesh generation [15]. MOAB also includes the DAGMC=
tool, supporting Monte Carlo radiation transport.
+
+The following code fragment shows very basic use of AdaptiveKDTree. A ran=
ge of entities is put in the tree; the leaf containing a given point is fou=
nd, and the entities in that leaf are returned.
+\code
+using namespace moab;
+// create the adaptive kd tree from a range of tets
+EntityHandle tree_root
+AdaptiveKDTree myTree(moab);
+ErrorCode rval =3D myTree.build_tree(tets, tree_root);
+
+// get the overall bounding box corners
+double boxmax[3], boxmin;
+rval =3D myTree.get_tree_box(tree_root, boxmax, boxmin);
+
+// get the tree leaf containing point xyz, and the tets in that leaf
+AdaptiveKDTreeIter treeiter;
+rval =3D myTree.leaf_containing_point(tree_root, xyz, treeiter);
+Range leaf_tets;
+rval =3D moab->get_entities_by_dimension(treeiter.handle(), 3,=20
+ leaf_tets, false);
+\endcode
+
+More detailed examples of using the various tree decompositions in MOAB ca=
n be found in [ref-treeexamples].
+
+ \ref contents
+
+ \subsection fourfive 4.5. File Reader/Writer Interfaces
+
+Mesh readers and writers communicate mesh into/out of MOAB from/to disk fi=
les. Reading a mesh often involves importing large sets of data, for examp=
le coordinates of all the nodes in the mesh. Normally, this process would =
involve reading data from the file into a temporary data buffer, then copyi=
ng data from there into its destination in MOAB. To avoid the expense of c=
opying data, MOAB has implemented a reader/writer interface that provides d=
irect access to blocks of memory used to represent mesh.
+
+The reader interface, declared in ReadUtilIface, is used to request blocks=
of memory for storing coordinate positions and element connectivity. The =
pointers returned from these functions point to the actual memory used to r=
epresent those data in MOAB. Once data is written to that memory, no furth=
er copying is done. This not only saves time, but it also eliminates the n=
eed to allocate a large memory buffer for intermediate storage of these dat=
a.=20
+
+MOAB allocates memory for nodes and elements (and their corresponding dens=
e tags) in chunks, to avoid frequent allocation/de-allocation of small chun=
ks of memory. The chunk size used depends on from where the mesh is being =
created, and can strongly affect the performance (and memory layout) of MOA=
B. Since dense tags are allocated at the chunk size, this can also affect =
overall memory usage in cases where the mesh size is small but the number o=
f dense tags or dense tag size is large. When creating vertices and elemen=
ts through the normal MOAB API, default chunk sizes defined in the Sequence=
Manager class are used. However, most of the file readers in MOAB allocate=
the exact amount of space necessary to represent the mesh being read. The=
re are also a few exceptions to this:
+
+- When compiled in parallel, this space is increased by a factor of 1.5, t=
o allow subsequent creation of ghost vertices/elements in the same chunk as=
the original mesh.
+- The .cub file reader, which creates nodes and elements for individual ge=
ometric entities in separate calls, allocates using the default vertex/elem=
ent sequence sizes, which are defined in the SequenceManager class in MOAB.
+.
+
+Applications calling the reader interface functions directly can specify t=
he allocation chunk size as an optional parameter.
+
+The reader interface consists of the following functions:
+
+- get_node_arrays: Given the number of vertices requested, the number of g=
eometric dimensions, and a requested start id, allocates a block of vertex =
handles and returns pointers to coordinate arrays in memory, along with the=
actual start handle for that block of vertices.
+- get_element_array: Given the number of elements requested, the number of=
vertices per element, the element type and the requested start id, allocat=
es the block of elements, and returns a pointer to the connectivity array f=
or those elements and the actual start handle for that block. The number o=
f vertices per element is necessary because those elements may include high=
er-order nodes, and MOAB stores these as part of the normal connectivity ar=
ray.
+- update_adjacencies: This function takes the start handle for a block of =
elements and the connectivity of those elements, and updates adjacencies fo=
r those elements. Which adjacencies are updated depends on the options set=
in AEntityFactory.
+.
+
+The following code fragment illustrates the use of ReadUtilIface to read a=
mesh directly into MOAB=E2=80=99s native representation. This code assume=
s that connectivity is specified in terms of vertex indices, with vertex in=
dices starting from 1.
+
+\code
+// get the read iface from moab
+ReadUtilIface *iface;
+ErrorCode rval =3D moab->get_interface("ReadUtilIface", &iface);
+
+// allocate a block of vertex handles and read xyz=E2=80=99s into them
+std::vector<double*> arrays;
+EntityHandle startv, *starth;
+rval =3D iface->get_node_arrays(3, num_nodes, 0, startv, arrays);
+for (int i =3D 0; i < num_nodes; i++)
+ infile >> arrays[0][i] >> arrays[1][i] >> arrays[2][i];
+
+// allocate block of hex handles and read connectivity into them
+rval =3D iface->get_element_array(num_hexes, 8, MBHEX, 0, starth);
+for (int i =3D 0; i < 8*num_hexes; i++)
+ infile >> starth[i];
+
+// change connectivity indices to vertex handles
+for (int i =3D 0; i < 8*num_hexes; i++)
+ starth[i] +=3D startv-1;
+\endcode
+
+The writer interface, declared in WriteUtilIface, provides functions that =
support writing vertex coordinates and element connectivity to storage loca=
tions input by the application. Assembling these data is a common task for=
writing mesh, and can be non-trivial when exporting only subsets of a mesh=
. The writer interface declares the following functions:
+
+- get_node_arrays: Given already-allocated memory and the number of vertic=
es and dimensions, and a range of vertices, this function writes vertex coo=
rdinates to that memory. If a tag is input, that tag is also written with =
integer vertex ids, starting with 1, corresponding to the order the vertice=
s appear in that sequence (these ids are used to write the connectivity arr=
ay in the form of vertex indices).
+- get_element_array: Given a range of elements and the tag holding vertex =
ids, and a pointer to memory, the connectivity of the specified elements ar=
e written to that memory, in terms of the indices referenced by the specifi=
ed tag. Again, the number of vertices per element is input, to allow the d=
irect output of higher-order vertices.
+- gather_nodes_from_elements: Given a range of elements, this function ret=
urns the range of vertices used by those elements. If a bit-type tag is in=
put, vertices returned are also marked with 0x1 using that tag. If no tag =
is input, the implementation of this function uses its own bit tag for mark=
ing, to avoid using an n2 algorithm for gathering vertices.
+- reorder: Given a permutation vector, this function reorders the connecti=
vity for entities with specified type and number of vertices per entity to =
match that permutation. This function is needed for writing connectivity i=
nto numbering systems other than that used internally in MOAB.
+.
+
+The following code fragment shows how to use WriteUtilIface to write the v=
ertex coordinates and connectivity indices for a subset of entities.
+
+\code
+using namespace moab;
+// get the write iface from moab
+WriteUtilIface *iface;
+ErrorCode rval =3D moab->get_interface("WriteUtilIface", &iface);
+
+// get all hexes the model, and choose the first 10 of those
+Range tmp_hexes, hexes, verts;
+rval =3D moab->get_entities_by_type(0, MBHEX, tmp_hexes);
+for (int i =3D 0; i < 10; i++) hexes.insert(tmp_hexes[i]);
+rval =3D iface->gather_nodes_from_elements(hexes, 0, verts);
+
+// assign vertex ids
+iface->assign_ids(verts, 0, 1);
+
+// allocate space for coordinates & write them
+std::vector<double*> arrays(3);
+for (int i =3D 0; i < 3; i++) arrays[i] =3D new double[verts.size()];
+iface->get_node_arrays(3, verts.size(), verts, 0, 1, arrays);
+
+// put connect=E2=80=99y in array, in the form of indices into vertex array
+std::vector<int> conn(8*hexes.size());
+iface->get_element_array(hexes.size(), 8, 0, hexes, 0, 1, &conn[0]);
+\endcode
+
+ \ref contents
+
+ \subsection foursix 4.6. File Readers/Writers Packaged With MOAB
+
+MOAB has been designed to efficiently represent data and metadata commonly=
found in finite element mesh files. Readers and writers are included with=
MOAB which import/export specific types of metadata in terms of MOAB sets =
and tags, as described earlier in this document. The number of readers and=
writers in MOAB will probably grow over time, and so they are not enumerat=
ed here. See the src/io/README file in the MOAB source distribution for a =
current list of supported formats.
+
+Because of its generic support for readers and writers, described in the p=
revious section, MOAB is also a good environment for constructing new mesh =
readers and writers. The ReadTemplate and WriteTemplate classes in src/io =
are useful starting points for constructing new file readers/writers; appli=
cations are encouraged to submit their own readers/writers for inclusion in=
MOAB=E2=80=99s contrib/io directory in the MOAB source.=20
+
+The usefulness of a file reader/writer is determined not only by its abili=
ty to read and write nodes and elements, but also in its ability to store t=
he various types of meta-data included with the typical mesh. MOAB readers=
and writers are distinguished by their ability to preserve meta-data in me=
shes that they read and write. For example, MOAB=E2=80=99s CUB reader impo=
rts not only the mesh saved from CUBIT, but also the grouping of mesh entit=
ies into sets which reflect the geometric topology of the model used to gen=
erate the mesh. See [4] for a more detailed description of meta-data conve=
ntions used in MOAB=E2=80=99s file readers and writers, and the individual =
file reader/writer header files in src/io for details about the specific re=
aders and writers.
+
+Three specific file readers in MOAB bear further discussion: MOAB=E2=80=99=
s native HDF5-based file reader/writer; the CUB reader, used to import mesh=
and meta-data represented in CUBIT; and the CGM reader, which imports geom=
etric models. These are described next.
+
+ \ref contents
+
+ \subsection native 4.6.1. Native HD5-Based Reader/Writer
+
+A mesh database must be able to save and restore the data stored in its da=
ta model, at least to the extent to which it understands the semantics of t=
hat data. MOAB defines an HDF5-based file format that can store data embed=
ded in MOAB. By convention, these files are given an =E2=80=9C.h5m=E2=80=
=9D file extension. When reading or writing large amounts of data, it is r=
ecommended to use this file format, as it is the most complete and also the=
most efficient of the file readers/writers in MOAB.=20
+
+ \subsection cub 4.6.2. CUB Reader
+
+CUBIT is a toolkit for generating tetrahedral and hexahedral finite elemen=
t meshes from solid model geometry [16]. This tool saves and restores data=
in a custom =E2=80=9C.cub=E2=80=9D file, which stores both mesh and geomet=
ry (and data relating the two). The CUB reader in MOAB can import and inte=
rpret much of the meta-data information saved in .cub files. Ref. [4] desc=
ribes the conventions used to store this meta-data in the MOAB data model. =
The information read from .cub files, and stored in the MOAB data model, i=
ncludes:
+
+- Geometric model entities and topology
+- Model entity names and ids
+- Groups, element blocks, nodesets, and sidesets, including model entities=
stored in them
+- Mesh scheme and interval size information assigned to model entities
+.
+
+Note that although information about model entities is recovered, MOAB by =
default does not depend on a solid modeling engine; this information is sto=
red in the form of entity sets and parent/child relations between them. Se=
e Ref. [4] for more information.
+
+ \ref contents
+
+ \subsection cgm 4.6.3. CGM Reader
+
+The Common Geometry Module (CGM) [17] is a library for representing solid =
model and other types of solid geometry data. The CUBIT mesh generation to=
olkit uses CGM for its geometric modeling support, and CGM can restore geom=
etric models in the exact state in which they were represented in CUBIT. M=
OAB contains a CGM reader, which can be enabled with a configure option. U=
sing this reader, MOAB can read geometric models, and represent their model=
topology using entity sets linked by parent/child relations. The mesh in =
these models comes directly from the modeling engine faceting routines; the=
se are the same facets used to visualize solid models in other graphics eng=
ines. When used in conjunction with the VisIt visualization tool (see Sect=
ion ), this provides a solution for visualizing geometric models. Xxx sho=
ws a model imported using MOAB=E2=80=99s CGM reader and visualized with Vis=
It.=20
+
+ \ref contents
+
+ \section parallel 5. Parallel Mesh Representation and Query
+
+A parallel mesh representation must strike a careful balance between provi=
ding an interface to mesh similar to that of a serial mesh, while also allo=
wing the discovery of parallel aspects of the mesh and performance of paral=
lel mesh-based operations efficiently. MOAB supports a spatial domain-deco=
mposed view of a parallel mesh, where each subdomain is assigned to a proce=
ssor, lower-dimensional entities on interfaces between subdomains are share=
d between processors, and ghost entities can be exchanged with neighboring =
processors. Locally, each subdomain, along with any locally-represented gh=
ost entities, are accessed through a local MOAB instance. Parallel aspects=
of the mesh, e.g. whether entities are shared, on an interface, or ghost e=
ntities, are embedded in the same data model (entities, sets, tags, interfa=
ce) used in the rest of MOAB. MOAB provides a suite of parallel functions =
for initializing and communicating with a parallel mesh, along with functio=
ns to query the parallel aspects of the mesh.
+
+ \ref contents
+
+ \subsection fiveone 5.1. Nomenclature & Representation
+
+Before discussing how to access parallel aspects of a mesh, several terms =
need to be defined: =20
+
+<B>Shared entity:</B> An entity shared by one or several other processors.
+
+<B>Multi-shared entity:</B> An entity shared by more than two processors.
+
+<B>Owning Processor:</B> Each shared entity is owned by exactly one proces=
sor. This processor has the right to set tag values on the entity and have=
those values propagated to any sharing processors. =20
+
+<B>Part:</B> The collection of entities assigned to a given processor. Wh=
en reading mesh in parallel, the entities in a Part, along with any lower-d=
imensional entities adjacent to those, will be read onto the assigned proce=
ssor.
+
+<B>Partition:</B> A collection of Parts which take part in parallel collec=
tive communication, usually associated with an MPI communicator.
+
+<B>Interface:</B> A collection of mesh entities bounding entities in multi=
ple parts. Interface entities are owned by a single processor, but are rep=
resented on all parts/processors they bound.
+
+<B>Ghost entity:</B> A shared, non-interface, non-owned entity.
+
+<B>Parallel status:</B> A characteristic of entities and sets represented =
in parallel. The parallel status, or =E2=80=9Cpstatus=E2=80=9D, is represen=
ted by a bit field stored in an unsigned character, with bit values as desc=
ribed in Table 4.
+
+ \subsection tablefour Table 4: Bits representing various parallel charac=
teristics of a mesh. Also listed are enumerated values that can be used in=
bitmask expressions; these enumerated variables are declared in MBParallel=
Conventions.h.
+
+<table border=3D"1">
+<tr>
+<th>Bit</th>
+<th>Name</th>
+<th>Represents</th>
+</tr>
+<tr>
+<td>0x1</td>
+<td>PSTATUS_NOT_OWNED</td>
+<td>Not owned by the local processor</td>
+</tr>
+<tr>
+<td>0x2</td>
+<td>PSTATUS_SHARED</td>
+<td>Shared by exactly two processorstd>
+</tr>
+<tr>
+<td>0x4</td>
+<td>PSTATUS_MULTISHARED</td>
+<td>Shared by three or more processors</td>
+</tr>
+<tr>
+<td>0x8</td>
+<td>PSTATUS_INTERFACE</td>
+<td>Part of lower-dimensional interface shared by multiple processors</td>
+</tr>
+<tr>
+<td>0x10</td>
+<td>PSTATUS_GHOST</td>
+<td>Non-owned, non-interface entities represented locally</td>
+</tr>
+</table>
+
+Parallel functionality is described in the following sections. First, met=
hods to load a mesh into a parallel representation are described; next, fun=
ctions for accessing parallel aspects of a mesh are described; functions fo=
r communicating mesh and tag data are described.
+
+ \ref contents
+
+ \subsection fivetwo 5.2. Parallel Mesh Initialization
+
+Parallel mesh is initialized in MOAB in several steps:
+
+-# Establish a local mesh on each processor, either by reading the mesh i=
nto that representation from disk, or by creating mesh locally through the =
normal MOAB interface. =20
+-# Find vertices, then other entities, shared between processors, based o=
n a globally-consistent vertex numbering stored on the GLOBAL_ID tag. =20
+-# Exchange ghost or halo elements within a certain depth of processor in=
terfaces with neighboring processors. =20
+.
+
+These steps can be executed by a single call to MOAB=E2=80=99s load_file f=
unction, using the procedure described in Section 5.2.1. Or, they can be e=
xecuted in smaller increments calling specific functions in MOAB=E2=80=99s =
ParallelComm class, as described in Section 5.2.2. Closely related to the =
latter method is the handling of communicators, described in more detail in=
Section.
+
+ \ref contents
+
+ \subsection initialization 5.2.1. Parallel Mesh Initialization by Loadin=
g a File
+
+In the file reading approach, a mesh must contain some definition of the p=
artition (the assignment of mesh, usually regions, to processors). Partiti=
ons can be derived from other set structures already on the mesh, or can be=
computed explicitly for that purpose by tools like mbzoltan (see Section 4=
.2). For example, geometric volumes used to generate the mesh, and region-=
based material type assignments, are both acceptable partitions (see Ref. [=
4] for information about this and other meta-data often accompanying mesh).=
In addition to defining the groupings of regions into parts, the assignme=
nt of specific parts to processors can be done implicitly or using addition=
al data stored with the partition.
+
+MOAB implements several specific methods for loading mesh into a parallel =
representation:
+
+- READ_PART: each processor reads only the mesh used by its part(s).
+
+- READ_DELETE: each processor reads the entire mesh, then deletes the mesh=
not used by its part(s).
+
+- BCAST_DELETE: the root processor reads and broadcasts the mesh; each pro=
cessor then deletes the mesh not used by its part(s).
+
+The READ_DELETE and BCAST_DELETE methods are supported for all file types =
MOAB is able to read, while READ_PART is only implemented for MOAB=E2=80=99=
s native HDF5-based file format.
+
+Various other options control the selection of part sets or other details =
of the parallel read process. For example, the application can specify the=
tags, and optionally tag values, which identify parts, and whether those p=
arts should be distributed according to tag value or using round-robin assi=
gnment.
+
+The options used to specify loading method, the data used to identify part=
s, and other parameters controlling the parallel read process, are shown in=
Table 5. =20
+ \subsection tablefive Table 5: Options passed to MOAB=E2=80=99s load_fil=
e function identifying the partition and other parameters controlling the p=
arallel read of mesh data. Options and values should appear in option stri=
ng as =E2=80=9Coption=3Dval=E2=80=9D, with a delimiter (usually =E2=80=9C;=
=E2=80=9D) between options.
+
+<table border=3D"1">
+<tr>
+<th>Option</th>
+<th>Value</th>
+<th>Description</th>
+</tr>
+<tr>
+<td>PARTITION</td>
+<td><tag_name></td>
+<td>Sets with the specified tag name should be used as part sets</td>
+</tr>
+<tr>
+<td>PARTITION_VAL</td>
+<td><val1, val2-val3, ...></td>
+<td>Integer values to be combined with tag name, with ranges input using v=
al1, val2-val3. Not meaningful unless PARTITION option is also given.</td>
+</tr>
+<tr>
+<td>PARTITION_DISTRIBUTE</td>
+<td>(none)</td>
+<td>If present, or values are not input using PARTITION_VAL, sets with tag=
indicated in PARTITION option are partitioned across processors in round-r=
obin fashion.</td>
+</tr>
+<tr>
+<td>PARALLEL_RESOLVE_SHARED_ENTS</td>
+<td><pd.sd></td>
+<td>Resolve entities shared between processors, where partition is made up=
of pd- dimensional entities, and entities of dimension sd and lower should=
be resolved.</td>
+</tr>
+<tr>
+<td>PARALLEL_GHOSTS</td>
+<td><gd.bd.nl[.ad]></td>
+<td>Exchange ghost elements at shared inter-processor interfaces. Ghost e=
lements of dimension gd will be exchanged. Ghost elements are chosen going=
through bd-dimensional interface entities. Number of layers of ghost elem=
ents is specified in nl. If ad is present, lower-dimensional entities boun=
ding exchanged ghost entities will also be exchanged; allowed values for ad=
are 1 (exchange bounding edges), 2 (faces), or 3 (edges and faces).</td>
+</tr>
+<tr>
+<td>CPUTIME</td>
+<td>(none)</td>
+<td>Print cpu time required for each step of parallel read & initializatio=
n.</td>
+</tr>
+<tr>
+<td>MPI_IO_RANK</td>
+<td><r></td>
+<td>If read method requires reading mesh onto a single processor, processo=
r with rank r is used to do that read.</td>
+</tr>
+</table>
+
+Several example option strings controlling parallel reading and initializa=
tion are:
+
+<B>=E2=80=9CPARALLEL=3DREAD_DELETE; PARTITION=3DMATERIAL_SET; PARTITION_VA=
L=3D100, 200, 600-700=E2=80=9D:</B> The whole mesh is read by every process=
or; this processor keeps mesh in sets assigned the tag whose name is =E2=80=
=9CMATERIAL_SET=E2=80=9D and whose value is any one of 100, 200, and 600-70=
0 inclusive.
+
+<B>=E2=80=9CPARALLEL=3DREAD_PART; PARTITION=3DPARALLEL_PARTITION, PARTITIO=
N_VAL=3D2=E2=80=9D:</B> Each processor reads only its mesh; this processor,=
whose rank is 2, is responsible for elements in a set with the PARALLEL_PA=
RTITION tag whose value is 2. This would by typical input for a mesh which=
had already been partitioned with e.g. Zoltan or Parmetis.
+
+<B>=E2=80=9CPARALLEL=3DBCAST_DELETE; PARTITION=3DGEOM_DIMENSION, PARTITION=
_VAL=3D3, PARTITION_DISTRIBUTE=E2=80=9D:</B> The root processor reads the f=
ile and broadcasts the whole mesh to all processors. If a list is construc=
ted with entity sets whose GEOM_DIMENSION tag is 3, i.e. sets corresponding=
to geometric volumes in the original geometric model, this processor is re=
sponsible for all elements with index R+iP, i >=3D 0 (i.e. a round-robin di=
stribution).
+
+ \ref contents
+
+ \subsection functions 5.2.2. Parallel Mesh Initialization Using Functions
+
+After creating the local mesh on each processor, an application can call t=
he following functions in ParallelComm to establish information on shared m=
esh entities. See the [ref-directparmesh example] in the MOAB source tree =
for a complete example of how this is done from an application.
+
+- ParallelComm::resolve_shared_entities (collective): Resolves shared enti=
ties between processors, based on GLOBAL_ID tag values of vertices. Variou=
s forms are available, based on entities to be evaluated and maximum dimens=
ion for which entity sharing should be found.
+
+- ParallelComm::exchange_ghost_cells (collective): Exchange ghost entities=
with processors sharing an interface with this processor, based on specifi=
ed ghost dimension (dimension of ghost entities exchanged), bridge dimensio=
n, number of layers, and type of adjacencies to ghost entities. An entity =
is sent as a ghost if it is within that number of layers, across entities o=
f the bridge dimension, with entities owned by the receiving processor, or =
if it is a lower-dimensional entity adjacent to a ghost entity and that opt=
ion is requested.
+.
+
+ \ref contents
+
+ \subsection communicator 5.2.3. Communicator Handling
+
+The ParallelComm constructor takes arguments for an MPI communicator and a=
MOAB instance. The ParallelComm instance stores the MPI communicator, and=
registers itself with the MOAB instance. Applications can specify the Par=
allelComm index to be used for a given file operation, thereby specifying t=
he MPI communicator for that parallel operation. For example:
+
+\code
+using namespace moab;
+// pass a communicator to the constructor, getting back the index
+MPI_Comm my_mpicomm;
+int pcomm_index;
+ParallelComm my_pcomm(moab, my_mpicomm, &pcomm_index);
+
+// write the pcomm index into a string option
+char load_opt[32];
+sprintf(load_opt, "PARALLEL=3DBCAST_DELETE;PARALLEL_COMM=3D%d",=20
+ pcomm_index);
+
+// specify that option in a parallel read operation
+ErrorCode rval =3D moab->load_file(load_opt, fname, ...)
+\endcode
+
+In the above code fragment, the ParallelComm instance with index pcomm_ind=
ex will be used in the parallel file read, so that the operation executes o=
ver the specified MPI communicator. If no ParallelComm instance is specifi=
ed for a parallel file operation, a default instance will be defined, using=
MPI_COMM_WORLD.
+
+Applications needing to retrieve a ParallelComm instance created previousl=
y and stored with the MOAB instance, e.g. by a different code component, ca=
n do so using a static function on ParallelComm:
+
+\code
+ParallelComm *my_pcomm =3D ParallelComm::get_pcomm(moab, pcomm_index);
+\endcode
+
+ParallelComm also provides the ParallelComm::get_all_pcomm function, for r=
etrieving all ParallelComm instances stored with a MOAB instance. For synt=
ax and usage of this function, see the MOAB online documentation for Parall=
elComm.hpp [8].
+
+ \ref contents
+
+ \subsection fivethree 5.3. Parallel Mesh Query Functions
+
+Various functions are commonly used in parallel mesh-based applications. =
Functions marked as being collective must be called collectively for all pr=
ocessors that are members of the communicator associated with the ParallelC=
omm instance used for the call.
+
+<B>ParallelComm::get_pstatus:</B> Get the parallel status for the entity.
+
+<B>ParallelComm::get_pstatus_entities:</B> Get all entities whose pstatus =
satisfies (pstatus & val).
+
+<B>ParallelComm::get_owner:</B> Get the rank of the owning processor for t=
he specified entity.
+
+<B>ParallelComm::get_owner_handle:</B> Get the rank of the owning processo=
r for the specified entity, and the entity's handle on the owning processor.
+
+<B>ParallelComm::get_sharing_data:</B> Get the sharing processor(s) and ha=
ndle(s) for an entity or entities. Various overloaded versions are availab=
le, some with an optional =E2=80=9Coperation=E2=80=9D argument, where Inter=
face::INTERSECT or Interface::UNION can be specified. This is similar to t=
he operation arguments to Interface::get_adjacencies.
+
+<B>ParallelComm::get_shared_entities:</B> Get entities shared with the giv=
en processor, or with all processors. This function has optional arguments=
for specifying dimension, whether interface entities are requested, and wh=
ether to return just owned entities.
+
+<B>ParallelComm::get_interface_procs:</B> Return all processors with whom =
this processor shares an interface.
+
+<B>ParallelComm::get_comm_procs:</B> Return all processors with whom this =
processor communicates.
+
+ \ref contents
+
+ \subsection fivefour 5.4. Parallel Mesh Communication
+
+Once a parallel mesh has been initialized, applications can call the Paral=
lelComm::exchange_tags function for exchanging tag values between processor=
s. This function causes the owning processor to send the specified tag val=
ues for all shared, owned entities to other processors sharing those entiti=
es. Asynchronous communication is used to hide latency, and only point-to-=
point communication is used in these calls.
+
+ \ref contents
+
+ \section applications 6. Building MOAB-Based Applications
+
+There are two primary mechanisms supported by MOAB for building applicatio=
ns, one based on MOAB-defined make variables, and the other based on the us=
e of libtool and autoconf. Both assume the use of a =E2=80=9Cmake=E2=80=9D=
-based build system. =20
+
+The easiest way to incorporate MOAB into an application=E2=80=99s build pr=
ocess is to include the =E2=80=9Cmoab.make=E2=80=9D file into the applicati=
on=E2=80=99s Makefile, adding the make variables MOAB_INCLUDES and MOAB_LIB=
S_LINK to application=E2=80=99s compile and link commands, respectively. M=
OAB_INCLUDES contains compiler options specifying the location of MOAB incl=
ude files, and any preprocessor definitions required by MOAB. MOAB_LIBS_LI=
NK contains both the options telling where libraries can be found, and the =
link options which incorporate those libraries into the application. Any l=
ibraries depended on by the particular configuration of MOAB are included i=
n that definition, e.g. the HDF5 library. Using this method to incorporate=
MOAB is the most straightforward; for example, the following Makefile is u=
sed to build one of the example problems packaged with the MOAB source:
+\code
+include ${MOAB_LIB_DIR}/moab.make
+
+GetEntities : GetEntities.o
+ ${CXX} $< ${MOAB_LIBS_LINK} -o $@
+
+.cpp.o :=20
+ ${CXX} ${MOAB_INCLUDES} -c $<
+\endcode
+
+Here, the MOAB_LIB_DIR environment variable or make argument definition sp=
ecifies where the MOAB library is installed; this is also the location of t=
he moab.make file. Once that file has been included, MOAB_INCLUDES and MOA=
B_LIBS_LINK can be used, as shown.
+
+Other make variables are defined in the moab.make file which simplify buil=
ding applications:
+
+- MOAB_LIBDIR, MOAB_INCLUDEDIR: the directories into which MOAB libraries =
and include files will be installed, respectively. Note that some include =
files are put in a subdirectory named =E2=80=9Cmoab=E2=80=9D below that dir=
ectory, to reflect namespace naming conventions used in MOAB.
+
+- MOAB_CXXFLAGS, MOAB_CFLAGS, MOAB_LDFLAGS: Options passed to the C++ and =
C compilers and the linker, respectively.
+
+- MOAB_CXX, MOAB_CC, MOAB_FC: C++, C, and Fortran compilers specified to M=
OAB at configure time, respectively.
+.
+
+The second method for incorporating MOAB into an application=E2=80=99s bui=
ld system is to use autoconf and libtool. MOAB is configured using these t=
ools, and generates the =E2=80=9C.la=E2=80=9D files that hold information o=
n library dependencies that can be used in application build systems also b=
ased on autoconf and libtool. Further information on this subject is beyon=
d the scope of this User=E2=80=99s Guide; see the =E2=80=9C.la=E2=80=9D fil=
es as installed by MOAB, and contact the MOAB developer=E2=80=99s mailing l=
ist [6] for more details.
+
+ \ref contents
+
+ \section implementation 7. iMesh (ITAPS Mesh Interface) Implementation =
in MOAB
+
+iMesh is a common API to mesh data developed as part of the Interoperable =
Tools for Advanced Petascale Simulations (ITAPS) project [18]. Application=
s using the iMesh interface can operate on any implementation of that inter=
face, including MOAB. MOAB-based applications can take advantage of other =
services implemented on top of iMesh, including the MESQUITE mesh improveme=
nt toolkit [19] and the GRUMMP mesh improvement library [20].
+
+MOAB=E2=80=99s native interface is accessed through the Interface abstract=
C++ base class. Wrappers are not provided in other languages; rather, app=
lications wanting to access MOAB from those languages should do so through =
iMesh. In most cases, the data models and functionality available through =
MOAB and iMesh are identical. However, there are a few differences, subtle=
and not-so-subtle, between the two:
+
+<B>SPARSE tags used by default:</B> MOAB=E2=80=99s iMesh implementation cr=
eates SPARSE tags by default, because of semantic requirements of other tag=
-related functions in iMesh. To create DENSE tags through iMesh, use the i=
Mesh_createTagWithOptions extension function (see below).
+
+<B>Higher-order elements:</B> ITAPS currently handles higher-order element=
s (e.g. a 10-node tetrahedron) usi[21]<sup>1</sup>. As described in [sec-e=
ntities], MOAB supports higher-order entities by allowing various numbers o=
f vertices to define topological entities like quadrilateral or tetrahedron=
. Applications can specify flags to the connectivity and adjacency functio=
ns specifying whether corner or all vertices are requested.
+
+<B>Self-adjacencies:</B> In MOAB=E2=80=99s native interface, entities are =
always self-adjacent<sup>2</sup>; that is, adjacencies of equal dimension r=
equested from an entity will always include that entity, while from iMesh w=
ill not include that entity.
+
+<B>Option strings:</B> The iMesh specification requires that options in th=
e options string passed to various functions (e.g. iMesh_load) be prepended=
with the implementation name required to parse them, and delimited with sp=
aces. Thus, a MOAB-targeted option would appear as =E2=80=9Cmoab:PARALLEL=
=3DREAD_PART moab:PARTITION=3DMATERIAL_SET=E2=80=9D.
+
+To provide complete MOAB support from other languages through iMesh, a col=
lection of iMesh extension functions are also available. A general descrip=
tion of these extensions appears below; for a complete description, see the=
online documentation for iMesh-extensions.h [8].
+
+- Recursive get_entities functions: There are many cases where sets includ=
e other sets (see [4] for more information). MOAB provides iMesh_getEntiti=
esRec, and other recursive-supporting functions, to get all non-set entitie=
s of a given type or topology accessible from input set(s). Similar functi=
ons are available for number of entities of a given type/topology.
+
+- Get entities by tag, and optionally tag value: It is common to search fo=
r entities with a given tag, and possibly tag value(s); functions like iMes=
h_getEntitiesByTag are provided for this purpose.
+
+- Options to createTag: To provide more control over the tag type, the iMe=
sh_createTagWithOptions is provided. The storage type is controlled with t=
he =E2=80=9C
+
+- MBCNType: Canonical numbering evaluations are commonly needed by applica=
tions, e.g. to apply boundary conditions locally. The MBCN package provide=
s these evaluations in terms of entity types defined in MOAB [9]; the getMB=
CNType is required to translate between iMesh_Topology and MBCN type.
+
+- Iterator step: Step an iterator a specified number of entities; allows a=
dvancement of an iterator without needing to allocate memory to hold the en=
tity handles stepped over.
+
+- Direct access to tag storage: The Interface::tag_iterate function allows=
an application get a pointer to the memory used to store a given tag. For=
dense tags on contiguous ranges of entities, this provides more efficient =
access to tags. The iMesh functionn iMesh_tagIterate provides access to th=
is functionlity. See examples/TagIterateC.c and examples/TagIterateF.F for=
examples of how to use this from C and Fortran, respectively.=20
+.
+
+As required by the iMesh specification, MOAB generates the =E2=80=9CiMesh-=
Defs.inc=E2=80=9D file and installs it with the iMesh and MOAB libraries. =
This file defines make variables which can be used to build iMesh-based app=
lications. The method used here is quite similar to that used for MOAB its=
elf (see Section 6). In particular, the IMESH_INCLUDES and IMESH_LIBS make=
variables can be used with application compile and link commands, respecti=
vely, with other make variables similar to those provided in moab.make also=
available.
+
+Note that using the iMesh interface from Fortran-based applications requir=
es a compiler that supports Cray pointers, along with the pass-by-value (%V=
AL) extension. Almost all compilers support those extensions; however, if =
using the gcc series of compilers, you must use gfortran 4.3 or later.
+
+<sup>1</sup>There are currently no implementations of this interface.
+
+<sup>2</sup>iMesh and MOAB both define adjacencies using the topological c=
oncept of closure. Since the closure of an entity includes the entity itse=
lf, the d-dimensional entities on the closure of a given entity should incl=
ude the entity itself.
+
+ \ref contents
+
+ \section representation 8. Structured Mesh Representation
+
+A structured mesh is defined as a D-dimensional mesh whose interior vertic=
es have 2D connected edges. Structured mesh can be stored without connect=
ivity, if certain information is kept about the parametric space of each st=
ructured block of mesh. MOAB can represent structured mesh with implicit c=
onnectivity, saving approximately 57% of the storage cost compared to an un=
structured representation<sup>1</sup>. Since connectivity must be computed=
on the fly, these queries execute a bit slower than those for unstructured=
mesh. More information on the theory behind MOAB's structured mesh repres=
entation can be found in=20
+
+Currently, MOAB's structured mesh representation can only be used by creat=
ing structured mesh at runtime; that is, structured mesh is saved/restored =
in an unstructured format in MOAB's HDF5-based native save format. For mor=
e details on how to use MOAB's structured mesh representation, see the scds=
eq_test.cpp source file in the test/ directory.
+
+<sup>1</sup> This assumes vertex coordinates are represented explicitly, a=
nd that there are approximately the same number of vertices and hexahedra i=
n a structured hex mesh.
+
+ \ref contents
+
+ \section element 9. Spectral Element Meshes
+
+The Spectral Element Method (SEM) is a high-order method, using a polynomi=
al Legendre interpolation basis with Gauss-Lobatto quadrature points, in co=
ntrast to the Lagrange basis used in (linear) finite elements [20]. SEM ob=
tains exponential convergence with decreasing mesh characteristic sizes, an=
d codes implementing this method typically have high floating-point intensi=
ty, making the method highly efficient on modern CPUs. Most Nth-order SEM =
codes require tensor product cuboid (quad/hex) meshes, with each d-dimensio=
nal element containing (N+1)d degrees of freedom (DOFs). There are various=
methods for representing SEM meshes and solution fields on them; this docu=
ment discusses these methods and the tradeoffs between them. The mesh part=
s of this discussion are given in terms of the iMesh mesh interface and its=
implementation by the MOAB mesh library [21].
+
+The figure above shows a two-dimensional 3rd-order SEM mesh consisting of =
four quadrilaterals. For this mesh, each quadrilateral has (N+1)^2=3D16 DO=
Fs, with corner and edge degrees of freedom shared between neighboring quad=
rilaterals.
+
+ \ref contents
+
+ \subsection nineone 9.1. Representations
+
+There are various representations of this mesh in a mesh database like MOA=
B, depending on how DOFs are related to mesh entities and tags on those ent=
ities. We mention several possible representations:
+
+-# Corner vertices, element-based DOFs: Each quadrilateral is defined by f=
our vertices, ordered in CCW order typical of FE meshes. DOFs are stored a=
s tags on quadrilaterals, with size (N+1)^2 values, ordered lexicographical=
ly (i.e. as a 2D array tag(i,j) with i varying faster than j.) In the figu=
re above, the connectivity for face 1 would be (1, 4, 16, 13), and DOFs wou=
ld be ordered (1..16). Note that in this representation, tag values for DO=
Fs shared by neighboring elements must be set multiple times, since there a=
re as many copies of these DOFs as elements sharing them.
+-# High-order FE-like elements: Each DOF is represented by a mesh vertex. =
Quadrilaterals each have (N+1)^2 vertices, ordered as they would be for hig=
h-order finite elements (corner vertices first, then mid-edge and mid-face =
elements; see [22]). Mid -face, -edge, and -region vertices for a given ed=
ge/face/region would be ordered lexicographically, according to positive di=
rection in a corresponding reference element. In the figure above, the con=
nectivity array for face 1 would be (1, 4, 16, 13, 2, 3, 8, 12, 14, 15, 5, =
9, 6, 7, 10, 11). DOF values are stored as tags on vertices. Since DOFs a=
re uniquely associated with vertices and vertices are shared by neighboring=
elements, tag values only need to be set once. Full vertex-quadrilateral =
adjacencies are available, for all vertices.
+-# Linear FE-like elements, one vertex per DOF, array with DOF vertices: E=
ach quadrilateral is defined by four (corner) vertices, with additional ver=
tices representing mid-edge and mid-face DOFs. An additional =E2=80=9CDOF =
array=E2=80=9D tag is assigned to each quadrilateral, storing the array of =
vertices representing the (N+1)^2 DOFs for the quadrilateral, ordered lexic=
ographically. For the figure above, the connectivity array for face 1 woul=
d be (1, 4, 16, 13), and the DOF array would be (1..16), assuming that vert=
ex handles are integers as shown in the figure. DOF values are stored as t=
ags on vertices, and lexicographically-ordered arrays of DOFs can be retrie=
ved using the DOF array tag as input to the tag_get_data function in MOAB. =
Adjacency functions would only be meaningful for corner vertices, but tag =
values would only need to be set once per DOF.
+-# High-order FE-like elements, array with DOF vertices: This is a combina=
tion of options 2 and 3. The advantage would be full vertex-quad adjacency=
support and direct availability of lexicographically-ordered vertex arrays=
, at the expense of more memory.
+-# Convert to linear mesh: Since a spectral element is a cuboid with highe=
r-order vertices, it can always be converted to N^2 linear cuboids using th=
e high-order vertices as corners of the finer quads/hexes. This is how rea=
ders in ParaView and VisIt typically import spectral meshes (CAM-SE also ex=
ports connectivity in this form).
+
+As a convenience for applications, functions could also be provided for im=
portant tasks, like assembling the vertex handles for an entity in lexograp=
hic order (useful for option 2 above), and getting an array of tag values i=
n lexicographic order (for option 3 above).
+
+ \ref contents
+
+ \subsection ninetwo 9.2. Tradeoffs
+
+There are various competing tradeoffs in the various representation types.=
These include:
+
+- Adjacencies: being able to retrieve the element(s) using a given (corner=
or higher-order) vertex.
+- Connectivity list: being able to retrieve the connectivity of a given el=
ement, consisting of all (corner + higher-order) vertices in the element, u=
sually in lexicographical order. This is closely linked with being able to=
access the connectivity list as a const*, i.e. using the list straight fro=
m memory without needing to copy it.
+- Memory vs. time: There is a memory vs. execution time tradeoff between d=
uplicating interface vertex solution/tag variables in neighboring elements =
(more memory but more time-efficient and allows direct access to tag storag=
e by applications) versus using vertex-based tags (less memory but requires=
assembly of variables into lexicographically-ordered arrays, and prevents =
direct access from applications).
+.
+
+The lower-memory option (storing variables on vertices and assembling into=
lexicographically-ordered arrays for application use) usually ends up cost=
ing more in memory anyway, since applications must allocate their own stora=
ge for these arrays. On the other hand, certain applications will always c=
hoose to do that, instead of sharing storage with MOAB for these variables.=
In the case where applications do share memory with MOAB, other tools wou=
ld need to interpret the lexicographically-ordered field arrays specially, =
instead of simply treating the vertex tags as a point-based field.
+
+ \ref contents
+
+ \subsection ninethree 9.3. MOAB Representation
+In choosing the right MOAB representation for spectral meshes, we are tryi=
ng to balance a) minimal memory usage, b) access to properly-ordered and -a=
ligned tag storage, and c) maximal compatibility with tools likely to use M=
OAB. The solution we propose is to use a representation most like option 2=
) above, with a few optional behaviors based on application requirements. =20
+
+In brief, we propose to represent elements using the linear, FE-ordered co=
nnectivity list (containing only corner vertices from the spectral element)=
, with field variables written to either vertices, lexicographically-ordere=
d arrays on elements, or both, and with a lexicographically-ordered array (=
stored on tag SPECTRAL_VERTICES) of all (corner+higher-order) vertices stor=
ed on elements. In the either/or case, the choice will be evident from the=
tag size and the entities on which the tag is set. In the both case, the =
tag name will have a =E2=80=9C-LEX=E2=80=9D suffix for the element tags, an=
d the size of the element tag will be (N+1)^2 times that of the vertex-base=
d tag. Finally, the file set containing the spectral elements (or the root=
set, if no file set was input to the read) will contain a =E2=80=9CSPECTRA=
L_ORDER=E2=80=9D tag whose value is N. These conventions are described in =
the =E2=80=9CMetadata Information=E2=80=9D document distributed with the MO=
AB source code.
+
+ \ref contents
+
+ \section performance 10. Performance and Using MOAB Efficiently from App=
lications
+
+MOAB is designed to operate efficiently on groups of entities and for larg=
e meshes. Applications will be most efficient when they operate on entitie=
s in groups, especially groups which are close in their order of creation. =
The MOAB API is structured to encourage operations on groups of entities. =
Conversely, MOAB will not perform as well as other libraries if there are =
frequent deletion and creation of entities. For those types of application=
s, a mesh library using a C++ object-based representation is more appropria=
te. In this section, performance of MOAB when executing a variety of tasks=
is described, and compared to that of other representations. Of course, t=
hese metrics are based on the particular models and environments where they=
are run, and may or may not be representative of other application types.
+
+One useful measure of MOAB performance is in the representation and query =
of a large mesh. MOAB includes a performance test, located in the test/per=
f directory, in which a single rectangular region of hexahedral elements is=
created then queried; the following steps are performed:
+
+- Create the vertices and hexes in the mesh
+- For each vertex, get the set of connected hexahedra
+- For each hex, get the connected vertices, their coordinates, average the=
m, and assign them as a tag on the hexes
+.
+
+This test can be run on your system to determine the runtime and memory pe=
rformance for these queries in MOAB.
+
+ \ref contents
+
+ \section conclusions 11. Conclusions and Future Plans
+
+MOAB, a Mesh-Oriented datABase, provides a simple but powerful data abstra=
ction to structured and unstructured mesh, and makes that abstraction avail=
able through a function API. MOAB provides the mesh representation for the=
VERDE mesh verification tool, which demonstrates some of the powerful mesh=
metadata representation capabilities in MOAB. MOAB includes modules that =
import mesh in the ExodusII, CUBIT .cub and Vtk file formats, as well as th=
e capability to write mesh to ExodusII, all without licensing restrictions =
normally found in ExodusII-based applications. MOAB also has the capabilit=
y to represent and query structured mesh in a way that optimizes storage sp=
ace using the parametric space of a structured mesh; see Ref. [17] for deta=
ils.
+
+Initial results have demonstrated that the data abstraction provided by MO=
AB is powerful enough to represent many different kinds of mesh data found =
in real applications, including geometric topology groupings and relations,=
boundary condition groupings, and inter-processor interface representation=
. Our future plans are to further explore how these abstractions can be us=
ed in the design through analysis process.
+
+ \ref contents
+
+ \section references 12. References
+
+[1] M. Fatenejad and G.A. Moses, =E2=80=9CCooper radiation hydrodynamics c=
ode..=E2=80=9D
+
+[2] T.J. Tautges and A. Caceres, =E2=80=9CScalable parallel solution coupl=
ing for multiphysics reactor simulation,=E2=80=9D Journal of Physics: Confe=
rence Series, vol. 180, 2009.
+
+[3] T.J. Tautges, MOAB Meta-Data Information, 2010.
+
+[4] T.J. Tautges, =E2=80=9CMOAB - ITAPS =E2=80=93 Trac.=E2=80=9D, http://t=
rac.mcs.anl.gov/projects/ITAPS/wiki/MOAB
+
+[5] =E2=80=9CMOAB Developers Email List.=E2=80=9D, moab-dev at mcs.anl.gov.
+
+[6] =E2=80=9CMOAB Users Email List.=E2=80=9D, moab at mcs.anl.gov.
+
+[7] =E2=80=9CMOAB online documentation.=E2=80=9D, http://gnep.mcs.anl.gov:=
8010/moab-docs/
+
+[8] T.J. Tautges, =E2=80=9CCanonical numbering systems for finite-element =
codes,=E2=80=9D Communications in Numerical Methods in Engineering, vol. O=
nline, Mar. 2009.
+
+[9] L.A. Schoof and V.R. Yarberry, EXODUS II: A Finite Element Data Model,=
Albuquerque, NM: Sandia National Laboratories, 1994.
+
+[10] M. PATRAN, =E2=80=9CPATRAN User=E2=80=99s Manual,=E2=80=9D 2005.
+
+[11] VisIt User's Guide.
+
+[12] K. Devine, E. Boman, R. Heaphy, B. Hendrickson, and C. Vaughan, =E2=
=80=9CZoltan Data Management Services for Parallel Dynamic Applications,=E2=
=80=9D Computing in Science and Engineering, vol. 4, 2002, pp. 90=E2=80=93=
97.
+
+[13] T.J. Tautges, P.P.H. Wilson, J. Kraftcheck, B.F. Smith, and D.L. Hend=
erson, =E2=80=9CAcceleration Techniques for Direct Use of CAD-Based Geomet=
ries in Monte Carlo Radiation Transport,=E2=80=9D International Conference =
on Mathematics, Computational Methods & Reactor Physics (M&C 2009), Sarato=
ga Springs, NY: American Nuclear Society, 2009.
+
+[14] H. Kim and T. Tautges, =E2=80=9CEBMesh: An Embedded Boundary Meshing =
Tool,=E2=80=9D in preparation.
+
+[15] G.D. Sjaardema, T.J. Tautges, T.J. Wilson, S.J. Owen, T.D. Blacker, W=
.J. Bohnhoff, T.L. Edwards, J.R. Hipp, R.R. Lober, and S.A. Mitchell, CUBIT=
mesh generation environment Volume 1: Users manual, Sandia National Labora=
tories, May 1994, 1994.
+
+[16] T.J. Tautges, =E2=80=9CCGM: A geometry interface for mesh generation,=
analysis and other applications,=E2=80=9D Engineering with Computers, vol=
. 17, 2001, pp. 299-314.
+
+[17] T. J. Tautges, MOAB-SD: Integrated structured and unstructured mesh r=
epresentation, Engineering with Computers, vol. 20, no. 3, pp. 286-293, 200=
4.
+
+[18] =E2=80=9CInteroperable Technologies for Advanced Petascale Simulation=
s (ITAPS),=E2=80=9D Interoperable Technologies for Advanced Petascale Simul=
ations (ITAPS).
+
+[19] P. Knupp, =E2=80=9CMesh quality improvement for SciDAC applications,=
=E2=80=9D Journal of Physics: Conference Series, vol. 46, 2006, pp. 458-46=
2.
+
+[20] M. O. Deville, P. F. Fischer, and E. H. Mund, High-order methods for =
incompressible fluid flow. Cambridge, UK; New York: Cambridge University Pr=
ess, 2002.
+
+[21] T. J. Tautges, =E2=80=9CMOAB Wiki.=E2=80=9D [Online]. Available: http=
://trac.mcs.anl.gov/projects/ITAPS/wiki/MOAB. [Accessed: 30-Oct-2012].
+
+[22] T. J. Tautges, =E2=80=9CCanonical numbering systems for finite-elemen=
t codes,=E2=80=9D International Journal for Numerical Methods in Biomedical=
Engineering, vol. 26, no. 12, pp. 1559=E2=80=931572, 2010.
+
+
+ \ref contents
+
+ \page musthave Differences Between iMesh and MOAB
+
+ The data models used in MOAB and iMesh are quite similar, but not identi=
cal.The most significant differences are the following:
+
+- Tags: MOAB differentiates between DENSE, SPARSE, and BIT tags, using dif=
ferent storage models for each, while iMesh uses a single tag concept. iMe=
sh allows application to query whether an entity has been given a tag of a =
specified type; this query is incompatible with the concept of a DENSE tag =
with a default value. Thus, MOAB=E2=80=99s iMesh implementation creates SP=
ARSE tags by default, and tags created and accessed through this interface =
will use more memory than DENSE tags created through MOAB=E2=80=99s native =
interface. To mitigate this problem, MOAB implements an extension of the i=
Mesh_createTag function which allows specification of the tag type (DENSE, =
SPARSE, etc.) to be created. See later in this section for more informatio=
n.
+
+- Higher-order nodes: ITAPS currently handles higher-order elements (e.g. =
a 10-node tetrahedron) using a special =E2=80=9CShape=E2=80=9D interface. =
In this interface, higher-order nodes are only accessible through the AEnti=
ties which they resolve. MOAB=E2=80=99s iMesh implementation provides acce=
ss to higher-order nodes in the same manner described in Section , by vary=
ing the number of vertices defining each entity. As a result, if higher-or=
der entities are used in a model, the functions returning connectivity and =
vertex adjacencies always return all vertices, rather than providing an opt=
ion to return just corner vertices.
+
+- Self-adjacencies: iMesh specifies that entities are not self-adjacent; t=
hat is, requesting adjacencies of the same dimension/type results in an err=
or. MOAB does not consider this an error, returning the entity itself.
+
+- Adjacency table and AEntities: iMesh uses the concept of an =E2=80=9Cadj=
acency table=E2=80=9D to determine which AEntities are available and create=
d by default. MOAB uses input arguments to the get_adjacencies functions t=
o control whether AEntities are created. These flags provide finer-grained=
control over AEntities, but make it slightly less convenient to ensure tha=
t AEntities of a given dimension are always created.
+.
+
+ \page figures List of Figures
+=20
+ This page is intended to be empty.
+=20
+ \page tables List of Tables
+=20
+ \ref tableone
+
+ \ref tabletwo
+
+ \ref tablethree
+
+ \ref tablefour
+
+ \ref tablefive
+=20
+ \page building Building & Installing
+=20
+ MOAB uses an autoconf and libtool-based build process by default. The p=
rocedure used to build MOAB from scratch depends on whether the source code=
was obtained from a =E2=80=9Ctarball=E2=80=9D or directly from the Subvers=
ion repository. Assuming the latter, the following steps should be execute=
d for building and installing MOAB:
+ - Locate and build any required dependencies. MOAB can be built with no=
dependencies on other libraries; this may be useful for applications only =
needing basic mesh representation and not needing to export mesh to formats=
implemented in other libraries. MOAB=E2=80=99s native save/restore capabi=
lity is built on HDF5-based files; applications needing to save and restore=
files from MOAB reliably should use this library. MOAB also uses ExodusII=
, a netCDF-based file format developed at Sandia National Laboratories [10]=
. Applications needing to execute these tests should also build netCDF. N=
ote that MOAB uses netCDF=E2=80=99s C++ interface, which is not enabled by =
default in netCDF but can be enabled using the =E2=80=9C=E2=80=93enable-cxx=
=E2=80=9D option to netCDF=E2=80=99s configure script.
+ - Unpack source code into <moab>, and change current working directory t=
o that location.
+ - Execute =E2=80=9Cautoreconf =E2=80=93fi=E2=80=9D.
+ - Run configure script, by executing =E2=80=9C./configure <options>=E2=
=80=9D. Recommended options:
+ -# =E2=80=93prefix=3D<install_dir>: directory below which MOAB libr=
ary and include files will be installed; can either be the directory used f=
or MOAB source (<moab> from step 1), or a different directory.
+ -# =E2=80=93hdf5-dir=3D<hdf5_dir>: directory whose =E2=80=9Cinclude=
=E2=80=9D and =E2=80=9Clib=E2=80=9D subdirectories hold HDF5 include and li=
brary, respectively. MOAB uses HDF5 for its native save/restore format (se=
e Section 4.6.1).
+ -# =E2=80=93netcdf-dir=3D: directory whose =E2=80=9Cinclude=E2=80=
=9D and =E2=80=9Clib=E2=80=9D subdirectories hold netCDF include and librar=
y, respectively. MOAB uses netCDF-based files for many of its build tests.=
If the location of netCDF cannot be found, MOAB=E2=80=99s build tests wil=
l not function properly, but MOAB will still be usable.
+ .
+ - Run =E2=80=9Cmake check=E2=80=9D; this runs a series of build tests, t=
o verify that the MOAB build was successful. Note this check will fail if =
netCDF is not used, but MOAB itself will still be usable from applications.
+ - Run =E2=80=9Cmake install=E2=80=9D; this copies include files and libr=
aries to subdirectories of the directory specified in the =E2=80=9Cprefix=
=E2=80=9D option.
+ .
+
+These steps are sufficient for building MOAB against HDF5 and netCDF. By =
default, a small number of standard MOAB-based applications are also built,=
including mbconvert (a utility for reading and writing files), mbsize (for=
querying basic information about a mesh), and the iMesh interface (see Sec=
tion 7). Other utilities can be enabled using various other options to the=
configure script; for a complete list of build options, execute =E2=80=9C.=
/configure =E2=80=93help=E2=80=9D.
+=20
+ */
+
https://bitbucket.org/fathomteam/moab/commits/60772c2789a8/
Changeset: 60772c2789a8
Branch: master
User: jhu
Date: 2013-06-06 18:23:43
Summary: Merge branch 'master' of https://bitbucket.org/fathomteam/moab
Affected #: 14 files
diff --git a/MeshFiles/unittest/io/Makefile.am b/MeshFiles/unittest/io/Make=
file.am
index ece20cf..8257e5f 100644
--- a/MeshFiles/unittest/io/Makefile.am
+++ b/MeshFiles/unittest/io/Makefile.am
@@ -3,6 +3,7 @@ EXTRA_DIST =3D HommeMapping.nc \
brick_cubit10.2.cub \
brick_cubit10.cub \
camEul26x48x96.t3.nc \
+ fv26x46x72.t.3.nc \
cubtest12.cub \
cubtest.jou \
dum.sat \
diff --git a/config/compiler.m4 b/config/compiler.m4
index 81596a6..6a7e2ed 100644
--- a/config/compiler.m4
+++ b/config/compiler.m4
@@ -102,15 +102,27 @@ if test "xno" !=3D "x$WITH_MPI"; then
DISTCHECK_CONFIGURE_FLAGS=3D"$DISTCHECK_CONFIGURE_FLAGS --with-mpi=3D\"$=
{withval}\""
=20
if test "xyes" =3D=3D "x$WITH_MPI"; then
+ if test "xno" !=3D "x$CHECK_CC"; then =20
FATHOM_SET_MPI_COMPILER([CC], [$CC_LIST])
+ fi
+ if test "xno" !=3D "x$CHECK_CXX"; then =20
FATHOM_SET_MPI_COMPILER([CXX],[$CXX_LIST])
- FATHOM_SET_MPI_COMPILER([FC], [$FC_LIST])
- FATHOM_SET_MPI_COMPILER([F77],[$F77_LIST])
+ fi
+ if test "xno" !=3D "x$CHECK_FC"; then =20
+ FATHOM_SET_MPI_COMPILER([FC], [$FC_LIST])
+ FATHOM_SET_MPI_COMPILER([F77],[$F77_LIST])
+ fi
else
+ if test "xno" !=3D "x$CHECK_CC"; then =20
FATHOM_SET_MPI_COMPILER([CC], [$CC_LIST],[${WITH_MPI}/bin])
+ fi
+ if test "xno" !=3D "x$CHECK_CXX"; then =20
FATHOM_SET_MPI_COMPILER([CXX],[$CXX_LIST],[${WITH_MPI}/bin])
- FATHOM_SET_MPI_COMPILER([FC], [$FC_LIST],[${WITH_MPI}/bin])
- FATHOM_SET_MPI_COMPILER([F77],[$F77_LIST],[${WITH_MPI}/bin])
+ fi
+ if test "xno" !=3D "x$CHECK_FC"; then
+ FATHOM_SET_MPI_COMPILER([FC], [$FC_LIST],[${WITH_MPI}/bin])
+ FATHOM_SET_MPI_COMPILER([F77],[$F77_LIST],[${WITH_MPI}/bin])
+ fi
WITH_MPI=3Dyes
fi
fi
diff --git a/src/io/Makefile.am b/src/io/Makefile.am
index 9c76524..b7eae6a 100644
--- a/src/io/Makefile.am
+++ b/src/io/Makefile.am
@@ -16,8 +16,12 @@ if NETCDF_FILE
MOAB_NETCDF_SRCS =3D ReadNCDF.cpp ReadNCDF.hpp \
WriteNCDF.cpp WriteNCDF.hpp \
WriteSLAC.cpp WriteSLAC.hpp \
- ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ ReadNC.cpp ReadNC.hpp \
+ NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
else
MOAB_NETCDF_SRCS =3D
endif
@@ -25,7 +29,11 @@ endif
if PNETCDF_FILE
if !NETCDF_FILE
MOAB_NETCDF_SRCS +=3D ReadNC.cpp ReadNC.hpp \
- ReadGCRM.cpp ReadGCRM.hpp
+ NCHelper.cpp NCHelper.hpp \
+ NCHelperEuler.cpp NCHelperEuler.hpp \
+ NCHelperFV.cpp NCHelperFV.hpp \
+ NCHelperHOMME.cpp NCHelperHOMME.hpp \
+ ReadGCRM.cpp ReadGCRM.hpp
endif
endif
=20
diff --git a/src/io/NCHelper.cpp b/src/io/NCHelper.cpp
new file mode 100644
index 0000000..e0f021e
--- /dev/null
+++ b/src/io/NCHelper.cpp
@@ -0,0 +1,21 @@
+#include "NCHelper.hpp"
+#include "NCHelperEuler.hpp"
+#include "NCHelperFV.hpp"
+#include "NCHelperHOMME.hpp"
+#include "moab/ReadUtilIface.hpp"
+
+namespace moab {
+
+NCHelper* NCHelper::get_nc_helper(ReadNC* readNC, int fileId, const FileOp=
tions& opts)
+{
+ if (NCHelperEuler::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperEuler(readNC, fileId);
+ else if (NCHelperFV::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperFV(readNC, fileId);
+ else if (NCHelperHOMME::can_read_file(readNC, fileId))
+ return new (std::nothrow) NCHelperHOMME(readNC, fileId, opts);
+ else // Unknown NetCDF grid (will fill this in later for POP, CICE and C=
LM)
+ return NULL;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelper.hpp b/src/io/NCHelper.hpp
new file mode 100644
index 0000000..852fc61
--- /dev/null
+++ b/src/io/NCHelper.hpp
@@ -0,0 +1,36 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelper.hpp
+//
+// Purpose : Climate NC file helper
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPER_HPP
+#define NCHELPER_HPP
+
+#include "ReadNC.hpp"
+
+namespace moab {
+
+//! Helper class to isolate reading of several different nc file formats
+class NCHelper
+{
+public:
+ NCHelper(ReadNC* readNC, int fileId) : _readNC(readNC), _fileId(fileId) =
{}
+
+ static NCHelper* get_nc_helper(ReadNC* readNC, int fileId, const FileOpt=
ions& opts);
+
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle f=
ile_set) =3D 0;
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptio=
ns& opts, EntityHandle file_set, Range& quads) =3D 0;
+ virtual std::string get_mesh_type_name() =3D 0;
+ virtual bool is_scd_mesh() =3D 0;
+
+protected:
+ ReadNC* _readNC;
+ int _fileId;
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
new file mode 100644
index 0000000..a01d4b8
--- /dev/null
+++ b/src/io/NCHelperEuler.cpp
@@ -0,0 +1,481 @@
+#include "NCHelperEuler.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS !=3D rval) {_readNC->readMeshIface->report_error("%s", =
str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_F=
AILURE;}
+
+namespace moab {
+
+bool NCHelperEuler::can_read_file(ReadNC* readNC, int fileId)
+{
+ std::vector<std::string>& dimNames =3D readNC->dimNames;
+
+ // If dimension names "lon" AND "lat' exist then it could be either the =
Eulerian Spectral grid or the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) !=
=3D dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) !=3D dimNames.end())) {
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then=
it should be the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("slon")) =
!=3D dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("slat")) !=3D dimNames.end()))
+ return false;
+
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt =3D readNC->glo=
balAtts.find("source");
+ if (attIt =3D=3D readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source=
global attribute.");
+ return false;
+ }
+ unsigned int sz =3D attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] =3D '\000';
+ int success =3D NCFUNC(get_att_text)(fileId, attIt->second.attVarId, a=
ttIt->second.attName.c_str(), &att_data[0]);
+ if (success !=3D 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source glo=
bal attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") =3D=3D std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHan=
dle file_set)
+{
+ Interface*& mbImpl =3D _readNC->mbImpl;
+ std::vector<std::string>& dimNames =3D _readNC->dimNames;
+ std::vector<int>& dimVals =3D _readNC->dimVals;
+ std::string& iName =3D _readNC->iName;
+ std::string& jName =3D _readNC->jName;
+ std::string& tName =3D _readNC->tName;
+ std::string& iCName =3D _readNC->iCName;
+ std::string& jCName =3D _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo =3D _readNC->varInfo;
+ int& tMin =3D _readNC->tMin;
+ int& tMax =3D _readNC->tMax;
+ int (&gDims)[6] =3D _readNC->gDims;
+ int (&lDims)[6] =3D _readNC->lDims;
+ int (&lCDims)[6] =3D _readNC->lCDims;
+ int (&gCDims)[6] =3D _readNC->gCDims;
+ std::vector<double>& ilVals =3D _readNC->ilVals;
+ std::vector<double>& jlVals =3D _readNC->jlVals;
+ std::vector<double>& tVals =3D _readNC->tVals;
+ std::vector<double>& ilCVals =3D _readNC->ilCVals;
+ std::vector<double>& jlCVals =3D _readNC->jlCVals;
+ int& tDim =3D _readNC->tDim;
+ int& iCDim =3D _readNC->iCDim;
+ int& jCDim =3D _readNC->jCDim;
+ DebugOutput& dbgOut =3D _readNC->dbgOut;
+ bool& isParallel =3D _readNC->isParallel;
+ int& partMethod =3D _readNC->partMethod;
+ int (&locallyPeriodic)[2] =3D _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] =3D _readNC->globallyPeriodic;
+ ScdParData& parData =3D _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm =3D _readNC->myPcomm;
+#endif
+
+ // look for names of center i/j dimensions
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ iCName =3D std::string("lon");
+ iName =3D std::string("slon");
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), iCName.c_str())=
) !=3D dimNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center i variable.");
+ }
+ iCDim =3D idx;
+
+ // decide on i periodicity using math for now
+ std::vector<double> tilVals(dimVals[idx]);
+ ErrorCode rval =3D _readNC->read_coordinate(iCName.c_str(), 0, dimVals[i=
dx] - 1, tilVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ if (std::fabs(2 * (*(tilVals.rbegin())) - *(tilVals.rbegin() + 1) - 360)=
< 0.001)
+ globallyPeriodic[0] =3D 1;
+
+ // now we can set gCDims and gDims for i
+ gCDims[0] =3D 0;
+ gDims[0] =3D 0;
+ gCDims[3] =3D dimVals[idx] - 1; // these are stored directly in file
+ gDims[3] =3D gCDims[3] + (globallyPeriodic[0] ? 0 : 1); // only if not p=
eriodic is vertex param max > elem param max
+
+ // now j
+ jCName =3D std::string("lat");
+ jName =3D std::string("slat");
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), jCName.c_str())=
) !=3D dimNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find center j variable.");
+ }
+ jCDim =3D idx;
+
+ // for Eul models, will always be non-periodic in j
+ gCDims[1] =3D 0;
+ gDims[1] =3D 0;
+ gCDims[4] =3D dimVals[idx] - 1;
+ gDims[4] =3D gCDims[4] + 1;
+
+ // try a truly 2d mesh
+ gDims[2] =3D -1;
+ gDims[5] =3D -1;
+
+ // look for time dimensions
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "time")) !=3D d=
imNames.end())
+ idx =3D vit - dimNames.begin();
+ else if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "t")) !=3D=
dimNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim =3D idx;
+ tMax =3D dimVals[idx] - 1;
+ tMin =3D 0;
+ tName =3D dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i =3D 0; i < 6; i++)
+ parData.gDims[i] =3D gDims[i];
+ for (int i =3D 0; i < 3; i++)
+ parData.gPeriodic[i] =3D globallyPeriodic[i];
+ parData.partMethod =3D partMethod;
+ int pdims[3];
+
+ rval =3D ScdInterface::compute_partition(myPcomm->proc_config().proc_s=
ize(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS !=3D rval)
+ return rval;
+ for (int i =3D 0; i < 3; i++)
+ parData.pDims[i] =3D pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() =3D=3D 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDim=
s[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i =3D 0; i < 6; i++)
+ lDims[i] =3D gDims[i];
+ locallyPeriodic[0] =3D globallyPeriodic[0];
+ }
+
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first=
resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # v=
ertex coords =3D # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] =3D lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3])=
{
+ // globally periodic and I'm the last proc, get fewer vertex coords =
than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] =3D lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] =3D lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] =3D lDims[0];
+ lCDims[4] =3D lDims[4] - 1;
+ lCDims[1] =3D lDims[1];
+
+ if (-1 !=3D lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 !=3D tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // now read coord values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (!ilCVals.empty()) {
+ if ((vmit =3D varInfo.find(iCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(iCName.c_str(), lDims[0], lDims[0]=
+ ilCVals.size() - 1, ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (!jlCVals.empty()) {
+ if ((vmit =3D varInfo.find(jCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(jCName.c_str(), lDims[1], lDims[1]=
+ jlCVals.size() - 1, jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] !=3D -1) {
+ if ((vmit =3D varInfo.find(iCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ double dif =3D (ilCVals[1] - ilCVals[0]) / 2;
+ std::size_t i;
+ for (i =3D 0; i !=3D ilCVals.size(); i++)
+ ilVals[i] =3D ilCVals[i] - dif;
+ // the last one is needed only if not periodic
+ if (!locallyPeriodic[0])
+ ilVals[i] =3D ilCVals[i - 1] + dif;
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] !=3D -1) {
+ if ((vmit =3D varInfo.find(jCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) =3D=3D (lDims[4] - lDims[1=
]))) {
+ std::string gwName("gw");
+ std::vector<double> gwVals(lDims[4] - lDims[1] - 1);
+ rval =3D _readNC->read_coordinate(gwName.c_str(), lDims[1], lDims[=
4] - 2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] =3D -(M_PI / 2) * 180 / M_PI;
+ unsigned int i =3D 0;
+ double gwSum =3D -1;
+ for (i =3D 1; i !=3D gwVals.size() + 1; i++) {
+ gwSum =3D gwSum + gwVals[i - 1];
+ jlVals[i] =3D std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] =3D 90.0; // using value of i after loop exits.
+ }
+ else {
+ std::string gwName("gw");
+ double gwSum =3D 0;
+
+ // If this is the first row
+ if (lDims[1] =3D=3D gDims[1]) {
+ std::vector<double> gwVals(lDims[4]);
+ rval =3D _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - =
1, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ jlVals[0] =3D -(M_PI / 2) * 180 / M_PI;
+ gwSum =3D -1;
+ for (std::size_t i =3D 1; i !=3D jlVals.size(); i++) {
+ gwSum =3D gwSum + gwVals[i - 1];
+ jlVals[i] =3D std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ // or if it's the last row
+ else if (lDims[4] =3D=3D gDims[4]) {
+ std::vector<double> gwVals(lDims[4] - 1);
+ rval =3D _readNC->read_coordinate(gwName.c_str(), 0, lDims[4] - =
2, gwVals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ // copy the correct piece
+ gwSum =3D -1;
+ for (int j =3D 0; j !=3D lDims[1] - 1; j++) {
+ gwSum =3D gwSum + gwVals[j];
+ }
+ std::size_t i =3D 0;
+ for (; i !=3D jlVals.size() - 1; i++) {
+ gwSum =3D gwSum + gwVals[lDims[1] - 1 + i];
+ jlVals[i] =3D std::asin(gwSum) * 180 / M_PI;
+ }
+ jlVals[i] =3D 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ int start =3D lDims[1] - 1;
+ int end =3D lDims[4] - 1;
+ std::vector<double> gwVals(end);
+ rval =3D _readNC->read_coordinate(gwName.c_str(), 0, end - 1, gw=
Vals);
+ ERRORR(rval, "Trouble reading gw variable.");
+ gwSum =3D -1;
+ for (int j =3D 0; j !=3D start - 1; j++) {
+ gwSum =3D gwSum + gwVals[j];
+ }
+ std::size_t i =3D 0;
+ for (; i !=3D jlVals.size(); i++) {
+ gwSum =3D gwSum + gwVals[start - 1 + i];
+ jlVals[i] =3D std::asin(gwSum) * 180 / M_PI;
+ }
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin !=3D -1) {
+ if ((vmit =3D varInfo.find(tName)) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=3D%d-%d, J=3D%d-%d\n", lDims[0], lDims[3], lDims[1]=
, lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * =
(lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit =3D varInfo.begin(); mit !=3D varInfo.end(); ++mit) {
+ ReadNC::VarData& vd =3D (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) !=3D vd.va=
rDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) !=3D vd.varDims.end()))
+ vd.entLoc =3D ReadNC::ENTLOCQUAD;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] =3D "__slon";
+ ijdimNames[1] =3D "__slat";
+ ijdimNames[2] =3D "__lon";
+ ijdimNames[3] =3D "__lat";
+
+ std::string tag_name;
+ int val_len =3D 0;
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ tag_name =3D ijdimNames[i];
+ void * val =3D NULL;
+ if (tag_name =3D=3D "__slon") {
+ val =3D &ilVals[0];
+ val_len =3D ilVals.size();
+ }
+ else if (tag_name =3D=3D "__slat") {
+ val =3D &jlVals[0];
+ val_len =3D jlVals.size();
+ }
+ else if (tag_name =3D=3D "__lon") {
+ val =3D &ilCVals[0];
+ val_len =3D ilCVals.size();
+ }
+ else if (tag_name =3D=3D "__lat") {
+ val =3D &jlCVals[0];
+ val_len =3D jlCVals.size();
+ }
+ Tag tagh =3D 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type =3D MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type =3D MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type =3D MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std:=
:endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, =
MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval =3D mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name =3D ss_tag_name.str();
+ Tag tagh =3D 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] =3D=3D "__slon") {
+ val[0] =3D lDims[0];
+ val[1] =3D lDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__slat") {
+ val[0] =3D lDims[1];
+ val[1] =3D lDims[4];
+ }
+ else if (ijdimNames[i] =3D=3D "__lon") {
+ val[0] =3D lCDims[0];
+ val[1] =3D lCDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__lat") {
+ val[0] =3D lCDims[1];
+ val[1] =3D lCDims[4];
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, =
tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.=
");
+ rval =3D mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MIN=
MAX tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name =3D ss_tag_name.str();
+ Tag tagh =3D 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] =3D=3D "__slon") {
+ val[0] =3D gDims[0];
+ val[1] =3D gDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__slat") {
+ val[0] =3D gDims[1];
+ val[1] =3D gDims[4];
+ }
+ else if (ijdimNames[i] =3D=3D "__lon") {
+ val[0] =3D gCDims[0];
+ val[1] =3D gCDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__lat") {
+ val[0] =3D gCDims[1];
+ val[1] =3D gCDims[4];
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, =
tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX t=
ag.");
+ rval =3D mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_=
MINMAX tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for variables like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperEuler::create_verts_quads(ScdInterface* scdi, const File=
Options& opts, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperEuler.hpp b/src/io/NCHelperEuler.hpp
new file mode 100644
index 0000000..b1c4c21
--- /dev/null
+++ b/src/io/NCHelperEuler.hpp
@@ -0,0 +1,36 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperEuler.hpp
+//
+// Purpose : Climate NC file helper for Eulerian Spectral grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPEREULER_HPP
+#define NCHELPEREULER_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Eulerian Spectral grid (CAM_EUL)
+class NCHelperEuler : public NCHelper
+{
+public:
+ NCHelperEuler(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle f=
ile_set);
+
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptio=
ns& opts, EntityHandle file_set, Range& quads);
+
+ virtual std::string get_mesh_type_name() { return "CAM_EUL"; }
+
+ virtual bool is_scd_mesh() { return true; }
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
new file mode 100644
index 0000000..566b719
--- /dev/null
+++ b/src/io/NCHelperFV.cpp
@@ -0,0 +1,477 @@
+#include "NCHelperFV.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+
+#include <cmath>
+#include <sstream>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS !=3D rval) {_readNC->readMeshIface->report_error("%s", =
str); return rval;}
+
+namespace moab {
+
+bool NCHelperFV::can_read_file(ReadNC* readNC, int fileId)
+{
+ std::vector<std::string>& dimNames =3D readNC->dimNames;
+
+ // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then i=
t should be the FV grid
+ if ((std::find(dimNames.begin(), dimNames.end(), std::string("lon")) !=
=3D dimNames.end()) && (std::find(dimNames.begin(),
+ dimNames.end(), std::string("lat")) !=3D dimNames.end()) && (std::fi=
nd(dimNames.begin(), dimNames.end(), std::string("slon"))
+ !=3D dimNames.end()) && (std::find(dimNames.begin(), dimNames.end(),=
std::string("slat")) !=3D dimNames.end())) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt =3D readNC->glo=
balAtts.find("source");
+ if (attIt =3D=3D readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source=
global attribute.");
+ return false;
+ }
+ unsigned int sz =3D attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] =3D '\000';
+ int success =3D NCFUNC(get_att_text)(fileId, attIt->second.attVarId, a=
ttIt->second.attName.c_str(), &att_data[0]);
+ if (success !=3D 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source glo=
bal attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") =3D=3D std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle=
file_set)
+{
+ Interface*& mbImpl =3D _readNC->mbImpl;
+ std::vector<std::string>& dimNames =3D _readNC->dimNames;
+ std::vector<int>& dimVals =3D _readNC->dimVals;
+ std::string& iName =3D _readNC->iName;
+ std::string& jName =3D _readNC->jName;
+ std::string& tName =3D _readNC->tName;
+ std::string& iCName =3D _readNC->iCName;
+ std::string& jCName =3D _readNC->jCName;
+ std::map<std::string, ReadNC::VarData>& varInfo =3D _readNC->varInfo;
+ int& tMin =3D _readNC->tMin;
+ int& tMax =3D _readNC->tMax;
+ int (&gDims)[6] =3D _readNC->gDims;
+ int (&lDims)[6] =3D _readNC->lDims;
+ int (&gCDims)[6] =3D _readNC->gCDims;
+ int (&lCDims)[6] =3D _readNC->lCDims;
+ std::vector<double>& ilVals =3D _readNC->ilVals;
+ std::vector<double>& jlVals =3D _readNC->jlVals;
+ std::vector<double>& tVals =3D _readNC->tVals;
+ std::vector<double>& ilCVals =3D _readNC->ilCVals;
+ std::vector<double>& jlCVals =3D _readNC->jlCVals;
+ int& iDim =3D _readNC->iDim;
+ int& jDim =3D _readNC->jDim;
+ int& tDim =3D _readNC->tDim;
+ int& iCDim =3D _readNC->iCDim;
+ int& jCDim =3D _readNC->jCDim;
+ DebugOutput& dbgOut =3D _readNC->dbgOut;
+ bool& isParallel =3D _readNC->isParallel;
+ int& partMethod =3D _readNC->partMethod;
+ int (&locallyPeriodic)[2] =3D _readNC->locallyPeriodic;
+ int (&globallyPeriodic)[2] =3D _readNC->globallyPeriodic;
+ ScdParData& parData =3D _readNC->parData;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm =3D _readNC->myPcomm;
+#endif
+
+ std::vector<std::string>::iterator vit;
+ unsigned int idx;
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "slon")) !=3D d=
imNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slon variable.");
+ }
+
+ iDim =3D idx;
+ gDims[3] =3D dimVals[idx] - 1;
+ gDims[0] =3D 0;
+ iName =3D dimNames[idx];
+
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "slat")) !=3D d=
imNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find slat variable.");
+ }
+ jDim =3D idx;
+ gDims[4] =3D dimVals[idx] - 1 + 2; // add 2 for the pole points
+ gDims[1] =3D 0;
+ jName =3D dimNames[idx];
+
+ // look for names of center i/j dimensions
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "lon")) !=3D di=
mNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon variable.");
+ }
+ iCDim =3D idx;
+ gCDims[3] =3D dimVals[idx] - 1;
+ gCDims[0] =3D 0;
+ iCName =3D dimNames[idx];
+
+ // check and set globallyPeriodic[0]
+ std::vector<double> til_vals(2);
+ ErrorCode rval =3D _readNC->read_coordinate(iCName.c_str(), dimVals[idx]=
- 2, dimVals[idx] - 1, til_vals);
+ ERRORR(rval, "Trouble reading slon variable.");
+ if (std::fabs(2 * til_vals[1] - til_vals[0] - 360) < 0.001)
+ globallyPeriodic[0] =3D 1;
+ if (globallyPeriodic[0])
+ assert("Number of vertices and edges should be same" && gDims[3] =3D=
=3D gCDims[3]);
+ else
+ assert("Number of vertices should equal to number of edges plus one" &=
& gDims[3] =3D=3D gCDims[3] + 1);
+
+#ifdef USE_MPI
+ // if serial, use a locally-periodic representation only if local mesh i=
s periodic, otherwise don't
+ if ((isParallel && myPcomm->proc_config().proc_size() =3D=3D 1) && globa=
llyPeriodic[0])
+ locallyPeriodic[0] =3D 1;
+#else
+ if (globallyPeriodic[0])
+ locallyPeriodic[0] =3D 1;
+#endif
+
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "lat")) !=3D di=
mNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat variable.");
+ }
+ jCDim =3D idx;
+ gCDims[4] =3D dimVals[idx] - 1;
+ gCDims[1] =3D 0;
+ jCName =3D dimNames[idx];
+
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "time")) !=3D d=
imNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim =3D idx;
+ tMax =3D dimVals[idx] - 1;
+ tMin =3D 0;
+ tName =3D dimNames[idx];
+
+ // parse options to get subset
+ if (isParallel) {
+#ifdef USE_MPI
+ for (int i =3D 0; i < 6; i++)
+ parData.gDims[i] =3D gDims[i];
+ for (int i =3D 0; i < 3; i++)
+ parData.gPeriodic[i] =3D globallyPeriodic[i];
+ parData.partMethod =3D partMethod;
+ int pdims[3];
+
+ rval =3D ScdInterface::compute_partition(myPcomm->proc_config().proc_s=
ize(),
+ myPcomm->proc_config().proc_rank(),
+ parData, lDims, locallyPeriodic, pdims);
+ if (MB_SUCCESS !=3D rval)
+ return rval;
+ for (int i =3D 0; i < 3; i++)
+ parData.pDims[i] =3D pdims[i];
+
+ dbgOut.tprintf(1, "Partition: %dx%d (out of %dx%d)\n",
+ lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
+ gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1);
+ if (myPcomm->proc_config().proc_rank() =3D=3D 0)
+ dbgOut.tprintf(1, "Contiguous chunks of size %d bytes.\n", 8 * (lDim=
s[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1));
+#endif
+ }
+ else {
+ for (int i =3D 0; i < 6; i++)
+ lDims[i] =3D gDims[i];
+ locallyPeriodic[0] =3D globallyPeriodic[0];
+ }
+ opts.get_int_option("IMIN", lDims[0]);
+ opts.get_int_option("IMAX", lDims[3]);
+ opts.get_int_option("JMIN", lDims[1]);
+ opts.get_int_option("JMAX", lDims[4]);
+
+ // now get actual coordinate values for vertices and cell centers; first=
resize
+ if (locallyPeriodic[0]) {
+ // if locally periodic, doesn't matter what global periodicity is, # v=
ertex coords =3D # elem coords
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0] + 1);
+ lCDims[3] =3D lDims[3];
+ }
+ else {
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3])=
{
+ // globally periodic and I'm the last proc, get fewer vertex coords =
than vertices in i
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] =3D lDims[3] - 1;
+ }
+ else {
+ ilVals.resize(lDims[3] - lDims[0] + 1);
+ ilCVals.resize(lDims[3] - lDims[0]);
+ lCDims[3] =3D lDims[3] - 1;
+ }
+ }
+
+ lCDims[0] =3D lDims[0];
+ lCDims[4] =3D lDims[4] - 1;
+ lCDims[1] =3D lDims[1];
+
+ if (-1 !=3D lDims[1]) {
+ jlVals.resize(lDims[4] - lDims[1] + 1);
+ jlCVals.resize(lCDims[4] - lCDims[1] + 1);
+ }
+
+ if (-1 !=3D tMin)
+ tVals.resize(tMax - tMin + 1);
+
+ // ... then read actual values
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (lCDims[0] !=3D -1) {
+ if ((vmit =3D varInfo.find(iCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(iCName.c_str(), lCDims[0], lCDims[=
3], ilCVals);
+ ERRORR(rval, "Trouble reading lon variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lon coordinate.");
+ }
+ }
+
+ if (lCDims[1] !=3D -1) {
+ if ((vmit =3D varInfo.find(jCName)) !=3D varInfo.end() && (*vmit).seco=
nd.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(jCName.c_str(), lCDims[1], lCDims[=
4], jlCVals);
+ ERRORR(rval, "Trouble reading lat variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find lat coordinate.");
+ }
+ }
+
+ if (lDims[0] !=3D -1) {
+ if ((vmit =3D varInfo.find(iName)) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ // last column
+ if (!locallyPeriodic[0] && globallyPeriodic[0] && lDims[3] > gDims[3=
]) {
+ til_vals.resize(ilVals.size() - 1, 0.0);
+ rval =3D _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3=
] - 1, til_vals);
+ double dif =3D til_vals[1] - til_vals[0];
+ std::size_t i;
+ for (i =3D 0; i !=3D til_vals.size(); i++)
+ ilVals[i] =3D til_vals[i];
+ ilVals[i] =3D ilVals[i - 1] + dif;
+ }
+ else {
+ rval =3D _readNC->read_coordinate(iName.c_str(), lDims[0], lDims[3=
], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ if (lDims[1] !=3D -1) {
+ if ((vmit =3D varInfo.find(jName)) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ if (!isParallel || ((gDims[4] - gDims[1]) =3D=3D (lDims[4] - lDims[1=
]))) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1] - 1);
+ rval =3D _readNC->read_coordinate(jName.c_str(), lDims[1], lDims[4=
] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] =3D -90.0;
+ unsigned int i =3D 0;
+ for (i =3D 1; i !=3D dummyVar.size() + 1; i++)
+ jlVals[i] =3D dummyVar[i - 1];
+ jlVals[i] =3D 90.0; // using value of i after loop exits.
+ }
+ else {
+ // If this is the first row
+ // need to read one less then available and read it into a dummy v=
ar
+ if (lDims[1] =3D=3D gDims[1]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval =3D _readNC->read_coordinate(jName.c_str(), lDims[1], lDims=
[4] - 1, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ jlVals[0] =3D -90.0;
+ for (int i =3D 1; i < lDims[4] + 1; i++)
+ jlVals[i] =3D dummyVar[i - 1];
+ }
+ // or if it's the last row
+ else if (lDims[4] =3D=3D gDims[4]) {
+ std::vector<double> dummyVar(lDims[4] - lDims[1]);
+ rval =3D _readNC->read_coordinate(jName.c_str(), lDims[1] - 1, l=
Dims[4] - 2, dummyVar);
+ ERRORR(rval, "Trouble reading y variable.");
+ // copy the correct piece
+ std::size_t i =3D 0;
+ for (i =3D 0; i !=3D dummyVar.size(); i++)
+ jlVals[i] =3D dummyVar[i];
+ jlVals[i] =3D 90.0; // using value of i after loop exits.
+ }
+ // it's in the middle
+ else {
+ rval =3D _readNC->read_coordinate(jCName.c_str(), lDims[1] - 1, =
lDims[4] - 1, jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (tMin !=3D -1) {
+ if ((vmit =3D varInfo.find(tName)) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ dbgOut.tprintf(1, "I=3D%d-%d, J=3D%d-%d\n", lDims[0], lDims[3], lDims[1]=
, lDims[4]);
+ dbgOut.tprintf(1, "%d elements, %d vertices\n", (lDims[3] - lDims[0]) * =
(lDims[4] - lDims[1]), (lDims[3] - lDims[0] + 1)
+ * (lDims[4] - lDims[1] + 1));
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit =3D varInfo.begin(); mit !=3D varInfo.end(); ++mit) {
+ ReadNC::VarData& vd =3D (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iCDim) !=3D vd.va=
rDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), jCDim) !=3D vd.varDims.end()))
+ vd.entLoc =3D ReadNC::ENTLOCQUAD;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jDim) !=3D v=
d.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iCDim) !=3D vd.varDims.end()))
+ vd.entLoc =3D ReadNC::ENTLOCNSEDGE;
+ else if ((std::find(vd.varDims.begin(), vd.varDims.end(), jCDim) !=3D =
vd.varDims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), iDim) !=3D vd.varDims.end()))
+ vd.entLoc =3D ReadNC::ENTLOCEWEDGE;
+ }
+
+ // <coordinate_dim_name>
+ std::vector<std::string> ijdimNames(4);
+ ijdimNames[0] =3D "__slon";
+ ijdimNames[1] =3D "__slat";
+ ijdimNames[2] =3D "__lon";
+ ijdimNames[3] =3D "__lat";
+
+ std::string tag_name;
+ int val_len =3D 0;
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ tag_name =3D ijdimNames[i];
+ void * val =3D NULL;
+ if (tag_name =3D=3D "__slon") {
+ val =3D &ilVals[0];
+ val_len =3D ilVals.size();
+ }
+ else if (tag_name =3D=3D "__slat") {
+ val =3D &jlVals[0];
+ val_len =3D jlVals.size();
+ }
+ else if (tag_name =3D=3D "__lon") {
+ val =3D &ilCVals[0];
+ val_len =3D ilCVals.size();
+ }
+ else if (tag_name =3D=3D "__lat") {
+ val =3D &jlCVals[0];
+ val_len =3D jlCVals.size();
+ }
+ Tag tagh =3D 0;
+ DataType data_type;
+
+ // assume all has same data type as lon
+ switch (varInfo["lon"].varDataType) {
+ case NC_BYTE:
+ case NC_CHAR:
+ case NC_DOUBLE:
+ data_type =3D MB_TYPE_DOUBLE;
+ break;
+ case NC_FLOAT:
+ data_type =3D MB_TYPE_DOUBLE;
+ break;
+ case NC_INT:
+ data_type =3D MB_TYPE_INTEGER;
+ break;
+ case NC_SHORT:
+ default:
+ std::cerr << "Unrecognized data type for tag " << tag_name << std:=
:endl;
+ ERRORR(MB_FAILURE, "Unrecognized data type");
+ break;
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 0, data_type, tagh, =
MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN);
+ ERRORR(rval, "Trouble creating <coordinate_dim_name> tag.");
+ rval =3D mbImpl->tag_set_by_ptr(tagh, &file_set, 1, &val, &val_len);
+ ERRORR(rval, "Trouble setting data for <coordinate_dim_name> tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_LOC_MINMAX
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
+ tag_name =3D ss_tag_name.str();
+ Tag tagh =3D 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] =3D=3D "__slon") {
+ val[0] =3D lDims[0];
+ val[1] =3D lDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__slat") {
+ val[0] =3D lDims[1];
+ val[1] =3D lDims[4];
+ }
+ else if (ijdimNames[i] =3D=3D "__lon") {
+ val[0] =3D lCDims[0];
+ val[1] =3D lCDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__lat") {
+ val[0] =3D lCDims[1];
+ val[1] =3D lCDims[4];
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, =
tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_LOC_MINMAX tag.=
");
+ rval =3D mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_LOC_MIN=
MAX tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // __<coordinate_dim_name>_GLOBAL_MINMAX
+ for (unsigned int i =3D 0; i !=3D ijdimNames.size(); i++) {
+ std::stringstream ss_tag_name;
+ ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
+ tag_name =3D ss_tag_name.str();
+ Tag tagh =3D 0;
+ std::vector<int> val(2, 0);
+ if (ijdimNames[i] =3D=3D "__slon") {
+ val[0] =3D gDims[0];
+ val[1] =3D gDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__slat") {
+ val[0] =3D gDims[1];
+ val[1] =3D gDims[4];
+ }
+ else if (ijdimNames[i] =3D=3D "__lon") {
+ val[0] =3D gCDims[0];
+ val[1] =3D gCDims[3];
+ }
+ else if (ijdimNames[i] =3D=3D "__lat") {
+ val[0] =3D gCDims[1];
+ val[1] =3D gCDims[4];
+ }
+ rval =3D mbImpl->tag_get_handle(tag_name.c_str(), 2, MB_TYPE_INTEGER, =
tagh, MB_TAG_SPARSE | MB_TAG_CREAT);
+ ERRORR(rval, "Trouble creating __<coordinate_dim_name>_GLOBAL_MINMAX t=
ag.");
+ rval =3D mbImpl->tag_set_data(tagh, &file_set, 1, &val[0]);
+ ERRORR(rval, "Trouble setting data for __<coordinate_dim_name>_GLOBAL_=
MINMAX tag.");
+ if (MB_SUCCESS =3D=3D rval)
+ dbgOut.tprintf(2, "Tag created for variable %s\n", tag_name.c_str());
+ }
+
+ // hack: create dummy tags, if needed, for dimensions like nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperFV::create_verts_quads(ScdInterface* scdi, const FileOpt=
ions& opts, EntityHandle file_set, Range& quads)
+{
+ return _readNC->create_scd_verts_quads(scdi, file_set, quads);
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperFV.hpp b/src/io/NCHelperFV.hpp
new file mode 100644
index 0000000..afcd9b8
--- /dev/null
+++ b/src/io/NCHelperFV.hpp
@@ -0,0 +1,32 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperFV.hpp
+//
+// Purpose : Climate NC file helper for Finite Volume grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERFV_HPP
+#define NCHELPERFV_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for Finite Volume grid (CAM_FV)
+class NCHelperFV : public NCHelper
+{
+public:
+ NCHelperFV(ReadNC* readNC, int fileId) : NCHelper(readNC, fileId) {}
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle f=
ile_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptio=
ns& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_FV"; }
+ virtual bool is_scd_mesh() { return true; }
+};
+
+} // namespace moab
+
+#endif
diff --git a/src/io/NCHelperHOMME.cpp b/src/io/NCHelperHOMME.cpp
new file mode 100644
index 0000000..32faa58
--- /dev/null
+++ b/src/io/NCHelperHOMME.cpp
@@ -0,0 +1,506 @@
+#include "NCHelperHOMME.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "FileOptions.hpp"
+#include "moab/SpectralMeshTool.hpp"
+
+#include <cmath>
+
+#define ERRORR(rval, str) \
+ if (MB_SUCCESS !=3D rval) {_readNC->readMeshIface->report_error("%s", =
str); return rval;}
+
+#define ERRORS(err, str) \
+ if (err) {_readNC->readMeshIface->report_error("%s", str); return MB_F=
AILURE;}
+
+namespace moab {
+
+NCHelperHOMME::NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions=
& opts) : NCHelper(readNC, fileId), _spectralOrder(-1)
+{
+ // Calculate spectral order
+ std::map<std::string, ReadNC::AttData>::iterator attIt =3D readNC->globa=
lAtts.find("np");
+ if (attIt !=3D readNC->globalAtts.end()) {
+ int success =3D NCFUNC(get_att_int)(readNC->fileId, attIt->second.attV=
arId, attIt->second.attName.c_str(), &_spectralOrder);
+ if (success !=3D 0)
+ readNC->readMeshIface->report_error("%s", "Failed to read np global =
attribute int data.");
+ else
+ _spectralOrder--; // Spectral order is one less than np
+
+ if (MB_SUCCESS =3D=3D opts.match_option("PARTITION_METHOD", "NODAL_PAR=
TITION"))
+ readNC->partMethod =3D -1;
+ }
+}
+
+bool NCHelperHOMME::can_read_file(ReadNC* readNC, int fileId)
+{
+ // If global attribute "np" exists then it should be the HOMME grid
+ if (readNC->globalAtts.find("np") !=3D readNC->globalAtts.end()) {
+ // Make sure it is CAM grid
+ std::map<std::string, ReadNC::AttData>::iterator attIt =3D readNC->glo=
balAtts.find("source");
+ if (attIt =3D=3D readNC->globalAtts.end()) {
+ readNC->readMeshIface->report_error("%s", "File does not have source=
global attribute.");
+ return false;
+ }
+ unsigned int sz =3D attIt->second.attLen;
+ std::string att_data;
+ att_data.resize(sz + 1);
+ att_data[sz] =3D '\000';
+ int success =3D NCFUNC(get_att_text)(fileId, attIt->second.attVarId, a=
ttIt->second.attName.c_str(), &att_data[0]);
+ if (success !=3D 0) {
+ readNC->readMeshIface->report_error("%s", "Failed to read source glo=
bal attribute char data.");
+ return false;
+ }
+ if (att_data.find("CAM") =3D=3D std::string::npos)
+ return false;
+
+ return true;
+ }
+
+ return false;
+}
+
+ErrorCode NCHelperHOMME::init_mesh_vals(const FileOptions& opts, EntityHan=
dle file_set)
+{
+ std::vector<std::string>& dimNames =3D _readNC->dimNames;
+ std::vector<int>& dimVals =3D _readNC->dimVals;
+ std::string& kName =3D _readNC->kName;
+ std::string& tName =3D _readNC->tName;
+ std::map<std::string, ReadNC::VarData>& varInfo =3D _readNC->varInfo;
+ int& tMin =3D _readNC->tMin;
+ int& tMax =3D _readNC->tMax;
+ int (&gDims)[6] =3D _readNC->gDims;
+ int (&lDims)[6] =3D _readNC->lDims;
+ int& iDim =3D _readNC->iDim;
+ int& kDim =3D _readNC->kDim;
+ int& tDim =3D _readNC->tDim;
+ std::vector<double>& ilVals =3D _readNC->ilVals;
+ std::vector<double>& jlVals =3D _readNC->jlVals;
+ std::vector<double>& klVals =3D _readNC->klVals;
+ std::vector<double>& tVals =3D _readNC->tVals;
+
+ ErrorCode rval;
+ unsigned int idx;
+ std::vector<std::string>::iterator vit;
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "time")) !=3D d=
imNames.end())
+ idx =3D vit - dimNames.begin();
+ else if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "t")) !=3D=
dimNames.end())
+ idx =3D vit - dimNames.begin();
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time variable.");
+ }
+ tDim =3D idx;
+ tMax =3D dimVals[idx] - 1;
+ tMin =3D 0;
+ tName =3D dimNames[idx];
+
+ // get number of vertices (labeled as number of columns) and levels
+ gDims[0] =3D gDims[3] =3D -1;
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "ncol")) !=3D d=
imNames.end()) {
+ idx =3D vit - dimNames.begin();
+ gDims[3] =3D dimVals[idx] - 1;
+ gDims[0] =3D 0;
+ iDim =3D idx;
+ }
+ if (-1 =3D=3D gDims[0])
+ return MB_FAILURE;
+
+ // set j coordinate to the number of quads
+ gDims[1] =3D gDims[0];
+ gDims[4] =3D gDims[3] - 2;
+
+ gDims[2] =3D gDims[5] =3D -1;
+ if ((vit =3D std::find(dimNames.begin(), dimNames.end(), "lev")) !=3D di=
mNames.end()) {
+ idx =3D vit - dimNames.begin();
+ gDims[5] =3D dimVals[idx] - 1, gDims[2] =3D 0, kName =3D std::string("=
lev");
+ kDim =3D idx;
+ }
+ if (-1 =3D=3D gDims[2])
+ return MB_FAILURE;
+
+ // read coordinate data
+ std::map<std::string, ReadNC::VarData>::iterator vmit;
+ if (gDims[0] !=3D -1) {
+ if ((vmit =3D varInfo.find("lon")) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate("lon", gDims[0], gDims[3], ilVals);
+ ERRORR(rval, "Trouble reading x variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find x coordinate.");
+ }
+ }
+
+ // store lat values in jlVals parameterized by j
+ if (gDims[1] !=3D -1) {
+ if ((vmit =3D varInfo.find("lat")) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate("lat", gDims[0], gDims[3], jlVals);
+ ERRORR(rval, "Trouble reading y variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
+ }
+ }
+
+ if (gDims[2] !=3D -1) {
+ if ((vmit =3D varInfo.find("lev")) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate("lev", gDims[2], gDims[5], klVals);
+ ERRORR(rval, "Trouble reading z variable.");
+
+ // decide whether down is positive
+ char posval[10];
+ int success =3D NCFUNC(get_att_text)(_readNC->fileId, (*vmit).second=
.varId, "positive", posval);
+ if (0 =3D=3D success && !strcmp(posval, "down")) {
+ for (std::vector<double>::iterator dvit =3D klVals.begin(); dvit !=
=3D klVals.end(); ++dvit)
+ (*dvit) *=3D -1.0;
+ }
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find z coordinate.");
+ }
+ }
+
+ if (tMin !=3D -1) {
+ if ((vmit =3D varInfo.find(tName)) !=3D varInfo.end() && (*vmit).secon=
d.varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+ }
+
+ if ((vmit =3D varInfo.find(tName)) !=3D varInfo.end() && (*vmit).second.=
varDims.size() =3D=3D 1) {
+ rval =3D _readNC->read_coordinate(tName.c_str(), tMin, tMax, tVals);
+ ERRORR(rval, "Trouble reading time variable.");
+ }
+ else {
+ ERRORR(MB_FAILURE, "Couldn't find time coordinate.");
+ }
+
+ // determine the entity location type of a variable
+ std::map<std::string, ReadNC::VarData>::iterator mit;
+ for (mit =3D varInfo.begin(); mit !=3D varInfo.end(); ++mit) {
+ ReadNC::VarData& vd =3D (*mit).second;
+ if ((std::find(vd.varDims.begin(), vd.varDims.end(), iDim) !=3D vd.var=
Dims.end()) && (std::find(vd.varDims.begin(),
+ vd.varDims.end(), kDim) !=3D vd.varDims.end()))
+ vd.entLoc =3D ReadNC::ENTLOCNODE;
+ }
+
+ std::copy(gDims, gDims + 6, lDims);
+
+ // don't read coordinates of columns until we actually create the mesh
+
+ // hack: create dummy tags, if needed, for variables like ncol and nbnd
+ // with no corresponding variables
+ _readNC->init_dims_with_no_cvars_info();
+
+ // This check is for HOMME and other ucd mesh. When ReadNC class instance
+ // gets out of scope in a script (and deleted), the localGid will be lost
+ rval =3D _readNC->check_ucd_localGid(file_set);
+ ERRORR(rval, "Trouble checking local Gid for ucd mesh.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode NCHelperHOMME::create_verts_quads(ScdInterface* scdi, const File=
Options& opts, EntityHandle file_set, Range& quads)
+{
+ Interface*& mbImpl =3D _readNC->mbImpl;
+ std::string& fileName =3D _readNC->fileName;
+ int& connectId =3D _readNC->connectId;
+ int (&gDims)[6] =3D _readNC->gDims;
+ int (&lDims)[6] =3D _readNC->lDims;
+ std::vector<double>& ilVals =3D _readNC->ilVals;
+ std::vector<double>& jlVals =3D _readNC->jlVals;
+ std::vector<double>& klVals =3D _readNC->klVals;
+ Tag& mGlobalIdTag =3D _readNC->mGlobalIdTag;
+ const Tag*& mpFileIdTag =3D _readNC->mpFileIdTag;
+ DebugOutput& dbgOut =3D _readNC->dbgOut;
+ bool& isParallel =3D _readNC->isParallel;
+ Range& localGid =3D _readNC->localGid;
+#ifdef USE_MPI
+ ParallelComm*& myPcomm =3D _readNC->myPcomm;
+#endif
+ bool& spectralMesh =3D _readNC->spectralMesh;
+
+ // need to get/read connectivity data before creating elements
+ std::string conn_fname;
+
+ // try to open the connectivity file through CONN option, if used
+ ErrorCode rval =3D opts.get_str_option("CONN", conn_fname);
+ if (MB_SUCCESS !=3D rval) {
+ // default convention for reading HOMME is a file HommeMapping.nc in s=
ame dir as data file
+ conn_fname =3D std::string(fileName);
+ size_t idx =3D conn_fname.find_last_of("/");
+ if (idx !=3D std::string::npos)
+ conn_fname =3D conn_fname.substr(0, idx).append("/HommeMapping.nc");
+ else
+ conn_fname =3D "HommeMapping.nc";
+ }
+
+ int success;
+
+ int rank, procs;
+#ifdef PNETCDF_FILE
+ if (isParallel) {
+ success =3D NCFUNC(open)(myPcomm->proc_config().proc_comm(), conn_fnam=
e.c_str(), 0, MPI_INFO_NULL, &connectId);
+ rank =3D myPcomm->proc_config().proc_rank();
+ procs =3D myPcomm->proc_config().proc_size();
+ }
+ else {
+ success =3D NCFUNC(open)(MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INF=
O_NULL, &connectId);
+ rank =3D 0;
+ procs =3D 1;
+ }
+#else
+ success =3D NCFUNC(open)(conn_fname.c_str(), 0, &connectId);
+ rank =3D 0;
+ procs =3D 1;
+#endif
+ ERRORS(success, "Failed on open.");
+
+ std::vector<std::string> conn_names;
+ std::vector<int> conn_vals;
+ rval =3D _readNC->get_dimensions(connectId, conn_names, conn_vals);
+ ERRORR(rval, "Failed to get dimensions for connectivity.");
+
+ if (conn_vals[0] !=3D gDims[3] - gDims[0] + 1 - 2) {
+ dbgOut.tprintf(1, "Warning: number of quads from %s and vertices from =
%s are inconsistent; nverts =3D %d, nquads =3D %d.\n",
+ conn_fname.c_str(), fileName.c_str(), gDims[3] - gDims[0] + 1, con=
n_vals[0]);
+ }
+
+ // read connectivity into temporary variable
+ int num_fine_quads, num_coarse_quads, start_idx;
+ std::vector<std::string>::iterator vit;
+ int idx;
+ if ((vit =3D std::find(conn_names.begin(), conn_names.end(), "ncells")) =
!=3D conn_names.end())
+ idx =3D vit - conn_names.begin();
+ else if ((vit =3D std::find(conn_names.begin(), conn_names.end(), "ncent=
ers")) !=3D conn_names.end())
+ idx =3D vit - conn_names.begin();
+ else {
+ ERRORR(MB_FAILURE, "Failed to get number of quads.");
+ }
+ int num_quads =3D conn_vals[idx];
+
+ // get the connectivity into tmp_conn2 and permute into tmp_conn
+ int cornerVarId;
+ success =3D NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId=
);
+ ERRORS(success, "Failed to get variable id.");
+ NCDF_SIZE tmp_dims[2] =3D {0, 0}, tmp_counts[2] =3D {4, static_cast<size=
_t>(num_quads)};
+ std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads);
+ success =3D NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_dims, tmp_co=
unts, &tmp_conn2[0] NCREQ);
+ ERRORS(success, "Failed to get temporary connectivity.");
+ success =3D NCFUNC(close)(connectId);
+ ERRORS(success, "Failed on close.");
+ // permute the connectivity
+ for (int i =3D 0; i < num_quads; i++) {
+ tmp_conn[4 * i] =3D tmp_conn2[i];
+ tmp_conn[4 * i + 1] =3D tmp_conn2[i + 1 * num_quads];
+ tmp_conn[4 * i + 2] =3D tmp_conn2[i + 2 * num_quads];
+ tmp_conn[4 * i + 3] =3D tmp_conn2[i + 3 * num_quads];
+ }
+
+ // need to know whether we'll be creating gather mesh later, to make sur=
e we allocate enough space
+ // in one shot
+ bool create_gathers =3D true;
+#ifdef USE_MPI
+ if (isParallel)
+ if (myPcomm->proc_config().proc_rank() !=3D 0)
+ create_gathers =3D false;
+#endif
+
+ // compute the number of local quads, accounting for coarse or fine repr=
esentation
+ // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
+ int spectral_unit =3D (spectralMesh ? _spectralOrder * _spectralOrder : =
1);
+ // num_coarse_quads is the number of quads instantiated in MOAB; if !spe=
ctralMesh, num_coarse_quads =3D num_fine_quads
+ num_coarse_quads =3D int(std::floor(1.0 * num_quads / (spectral_unit * p=
rocs)));
+ // start_idx is the starting index in the HommeMapping connectivity list=
for this proc, before converting to coarse quad representation
+ start_idx =3D 4 * rank * num_coarse_quads * spectral_unit;
+ // iextra =3D # coarse quads extra after equal split over procs
+ int iextra =3D num_quads % (procs * spectral_unit);
+ if (rank < iextra)
+ num_coarse_quads++;
+ start_idx +=3D 4 * spectral_unit * std::min(rank, iextra);
+ // num_fine_quads is the number of quads in the connectivity list in Hom=
meMapping file assigned to this proc
+ num_fine_quads =3D spectral_unit * num_coarse_quads;
+
+ // now create num_coarse_quads
+ EntityHandle *conn_arr;
+ EntityHandle start_vertex;
+ Range tmp_range;
+
+ // read connectivity into that space
+ EntityHandle *sv_ptr =3D NULL, start_quad;
+ SpectralMeshTool smt(mbImpl, _spectralOrder);
+ if (!spectralMesh) {
+ rval =3D _readNC->readMeshIface->get_element_connect(num_coarse_quads,=
4,
+ MBQUAD, 0, start_quad, conn_=
arr,
+ // might have to create ga=
ther mesh later
+ (create_gathers ? num_coarse=
_quads + num_quads : num_coarse_quads));
+ ERRORR(rval, "Failed to create quads.");
+ tmp_range.insert(start_quad, start_quad + num_coarse_quads - 1);
+ std::copy(&tmp_conn[start_idx], &tmp_conn[start_idx + 4 * num_fine_qua=
ds], conn_arr);
+ std::copy(conn_arr, conn_arr + 4 * num_fine_quads, range_inserter(loca=
lGid));
+ }
+ else {
+ rval =3D smt.create_spectral_elems(&tmp_conn[0], num_fine_quads, 2, tm=
p_range, start_idx, &localGid);
+ ERRORR(rval, "Failed to create spectral elements.");
+ int count, v_per_e;
+ rval =3D mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), c=
onn_arr, v_per_e, count);
+ ERRORR(rval, "Failed to get connectivity of spectral elements.");
+ rval =3D mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_rang=
e.begin(), tmp_range.end(),
+ count, (void*&)sv_ptr);
+ ERRORR(rval, "Failed to get fine connectivity of spectral elements.");
+ }
+
+ // on this proc, I get columns ldims[1]..ldims[4], inclusive; need to fi=
nd which vertices those correpond to
+ unsigned int num_local_verts =3D localGid.size();
+ unsigned int num_total_verts =3D gDims[3] - gDims[0] + 1;
+
+ // create vertices
+ std::vector<double*> arrays;
+ rval =3D _readNC->readMeshIface->get_node_coords(3, num_local_verts, 0, =
start_vertex, arrays,
+ // might have to create gather m=
esh later
+ (create_gathers ? num_local_verts+=
num_total_verts : num_local_verts));
+ ERRORR(rval, "Couldn't create vertices in ucd mesh.");
+
+ // set vertex coordinates
+ Range::iterator rit;
+ double *xptr =3D arrays[0], *yptr =3D arrays[1], *zptr =3D arrays[2];
+ int i;
+ for (i =3D 0, rit =3D localGid.begin(); i < (int)num_local_verts; i++, +=
+rit) {
+ assert(*rit < ilVals.size() + 1);
+ xptr[i] =3D ilVals[(*rit) - 1];
+ yptr[i] =3D jlVals[(*rit) - 1];
+ zptr[i] =3D klVals[lDims[2]];
+ }
+
+ const double pideg =3D acos(-1.0) / 180.0;
+ for (i =3D 0; i < (int)num_local_verts; i++) {
+ double cosphi =3D cos(pideg * yptr[i]);
+ double zmult =3D sin(pideg * yptr[i]), xmult =3D cosphi * cos(xptr[i] =
* pideg), ymult =3D cosphi * sin(xptr[i] * pideg);
+ double rad =3D 8.0e3 + klVals[lDims[2]];
+ xptr[i] =3D rad * xmult, yptr[i] =3D rad * ymult, zptr[i] =3D rad * zm=
ult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range vert_range(start_vertex, start_vertex + num_local_verts - 1);
+ void *data;
+ int count;
+ rval =3D mbImpl->tag_iterate(mGlobalIdTag, vert_range.begin(), vert_rang=
e.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count =3D=3D (int) num_local_verts);
+ int *gid_data =3D (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ // duplicate global id data, which will be used to resolve sharing
+ if (mpFileIdTag) {
+ rval =3D mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_ra=
nge.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator on file id tag.");
+ assert(count =3D=3D (int) num_local_verts);
+ gid_data =3D (int*) data;
+ std::copy(localGid.begin(), localGid.end(), gid_data);
+ }
+
+ // create map from file ids to vertex handles, used later to set connect=
ivity
+ std::map<EntityHandle, EntityHandle> vert_handles;
+ for (rit =3D localGid.begin(), i =3D 0; rit !=3D localGid.end(); ++rit, =
i++) {
+ vert_handles[*rit] =3D start_vertex + i;
+ }
+
+ // compute proper handles in connectivity using offset
+ for (int q =3D 0; q < 4 * num_coarse_quads; q++) {
+ conn_arr[q] =3D vert_handles[conn_arr[q]];
+ assert(conn_arr[q]);
+ }
+ if (spectralMesh) {
+ int verts_per_quad =3D (_spectralOrder + 1) * (_spectralOrder + 1);
+ for (int q =3D 0; q < verts_per_quad * num_coarse_quads; q++) {
+ sv_ptr[q] =3D vert_handles[sv_ptr[q]];
+ assert(sv_ptr[q]);
+ }
+ }
+
+ // add new vertices and elements to the set
+ quads.merge(tmp_range);
+ tmp_range.insert(start_vertex, start_vertex + num_local_verts - 1);
+ rval =3D mbImpl->add_entities(file_set, tmp_range);
+ ERRORR(rval, "Couldn't add new vertices and quads/hexes to file set.");
+
+ // mark the set with the spectral order
+ Tag sporder;
+ rval =3D mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sp=
order, MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create spectral order tag.");
+ rval =3D mbImpl->tag_set_data(sporder, &file_set, 1, &_spectralOrder);
+ ERRORR(rval, "Couldn't set value for spectral order tag.");
+
+#ifdef USE_MPI
+ if (isParallel && myPcomm->proc_config().proc_rank() =3D=3D 0) {
+#endif
+ EntityHandle gather_set;
+ rval =3D mbImpl->create_meshset(MESHSET_SET, gather_set);
+ ERRORR(rval, "Trouble creating gather set.");
+
+ // create vertices
+ arrays.clear();
+ // don't need to specify allocation number here, because we know enoug=
h verts were created before
+ rval =3D _readNC->readMeshIface->get_node_coords(3, num_total_verts, 0=
, start_vertex, arrays);
+ ERRORR(rval, "Couldn't create vertices in ucd mesh for gather set.");
+
+ xptr =3D arrays[0], yptr =3D arrays[1], zptr =3D arrays[2];
+ for (i =3D 0; i < (int)num_total_verts; i++) {
+ double cosphi =3D cos(pideg * jlVals[i]);
+ double zmult =3D sin(pideg * jlVals[i]);
+ double xmult =3D cosphi * cos(ilVals[i] * pideg);
+ double ymult =3D cosphi * sin(ilVals[i] * pideg);
+ double rad =3D 8.0e3 + klVals[lDims[2]];
+ xptr[i] =3D rad * xmult;
+ yptr[i] =3D rad * ymult;
+ zptr[i] =3D rad * zmult;
+ }
+
+ // get ptr to gid memory for vertices
+ Range gather_verts(start_vertex, start_vertex + num_total_verts - 1);
+ rval =3D mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gathe=
r_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator.");
+ assert(count =3D=3D (int) num_total_verts);
+ gid_data =3D (int*) data;
+ for (int j =3D 1; j <=3D (int) num_total_verts; j++)
+ gid_data[j - 1] =3D j;
+ // set the file id tag too, it should be bigger something not interfer=
ing with global id
+ if (mpFileIdTag) {
+ rval =3D mbImpl->tag_iterate(*mpFileIdTag, gather_verts.begin(), gat=
her_verts.end(), count, data);
+ ERRORR(rval, "Failed to get tag iterator in file id tag.");
+ assert(count =3D=3D (int) num_total_verts);
+ gid_data =3D (int*) data;
+ for (int j =3D 1; j <=3D (int) num_total_verts; j++)
+ gid_data[j - 1] =3D num_total_verts + j; // bigger than global id =
tag
+ }
+
+ rval =3D mbImpl->add_entities(gather_set, gather_verts);
+ ERRORR(rval, "Couldn't add vertices to gather set.");
+
+ // create quads
+ Range gather_quads;
+ // don't need to specify allocation number here, because we know enoug=
h quads were created before
+ rval =3D _readNC->readMeshIface->get_element_connect(num_quads, 4,
+ MBQUAD, 0, start_quad, conn_=
arr);
+ ERRORR(rval, "Failed to create quads.");
+ gather_quads.insert(start_quad, start_quad + num_quads - 1);
+ std::copy(&tmp_conn[0], &tmp_conn[4 * num_quads], conn_arr);
+ for (i =3D 0; i !=3D 4 * num_quads; i++)
+ conn_arr[i] +=3D start_vertex - 1; // connectivity array is shifted =
by where the gather verts start
+ rval =3D mbImpl->add_entities(gather_set, gather_quads);
+ ERRORR(rval, "Couldn't add quads to gather set.");
+
+ Tag gathersettag;
+ rval =3D mbImpl->tag_get_handle("GATHER_SET", 1, MB_TYPE_INTEGER, gath=
ersettag,
+ MB_TAG_CREAT | MB_TAG_SPARSE);
+ ERRORR(rval, "Couldn't create gather set tag.");
+ int gatherval =3D 1;
+ rval =3D mbImpl->tag_set_data(gathersettag, &gather_set, 1, &gatherval=
);
+ ERRORR(rval, "Couldn't set value for gather set tag.");
+
+#ifdef USE_MPI
+ }
+#endif
+
+ return MB_SUCCESS;
+}
+
+} // namespace moab
diff --git a/src/io/NCHelperHOMME.hpp b/src/io/NCHelperHOMME.hpp
new file mode 100644
index 0000000..feb3cad
--- /dev/null
+++ b/src/io/NCHelperHOMME.hpp
@@ -0,0 +1,35 @@
+//-------------------------------------------------------------------------
+// Filename : NCHelperHOMME.hpp
+//
+// Purpose : Climate NC file helper for HOMME grid
+//
+// Creator : Danqing Wu
+//-------------------------------------------------------------------------
+
+#ifndef NCHELPERHOMME_HPP
+#define NCHELPERHOMME_HPP
+
+#include "NCHelper.hpp"
+
+namespace moab {
+
+//! Child helper class for HOMME grid (CAM_SE)
+class NCHelperHOMME : public NCHelper
+{
+public:
+ NCHelperHOMME(ReadNC* readNC, int fileId, const FileOptions& opts);
+ static bool can_read_file(ReadNC* readNC, int fileId);
+
+private:
+ virtual ErrorCode init_mesh_vals(const FileOptions& opts, EntityHandle f=
ile_set);
+ virtual ErrorCode create_verts_quads(ScdInterface* scdi, const FileOptio=
ns& opts, EntityHandle file_set, Range& quads);
+ virtual std::string get_mesh_type_name() { return "CAM_SE"; }
+ virtual bool is_scd_mesh() { return false; }
+
+private:
+ int _spectralOrder; // read from variable 'np'
+};
+
+} // namespace moab
+
+#endif
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list