[MOAB-dev] r3508 - in MOAB/trunk: . m4
tautges at mcs.anl.gov
tautges at mcs.anl.gov
Thu Jan 28 23:24:01 CST 2010
Author: tautges
Date: 2010-01-28 23:24:00 -0600 (Thu, 28 Jan 2010)
New Revision: 3508
Added:
MOAB/trunk/ReadCCMIO.cpp
MOAB/trunk/ReadCCMIO.hpp
MOAB/trunk/WriteCCMIO.cpp
MOAB/trunk/WriteCCMIO.hpp
MOAB/trunk/m4/ccmio.m4
Modified:
MOAB/trunk/MBCN.cpp
MOAB/trunk/MBCN.hpp
MOAB/trunk/MBReaderWriterSet.cpp
MOAB/trunk/Makefile.am
MOAB/trunk/ReadSmf.cpp
MOAB/trunk/SMF_State.cpp
MOAB/trunk/configure.ac
Log:
Adding a CCMIO reader/writer. Basic functionality for reading hex meshes works.
No metadata processing yet.
Specific changes:
ReadSmf.cpp: commenting out open so make check will work
SMF_State.cpp: adding some includes for gcc4.4
configure.ac, m4/ccmio.m4, Makefile.am: processing options for ccmio reader
MBCN.cpp: un-commenting SubEntityConn function; why was this commented out in the first place?
MBReaderWriterSet, ReadCCMIO, WriteCCMIO: CCMIO reader/writer
Make check works now, again.
Modified: MOAB/trunk/MBCN.cpp
===================================================================
--- MOAB/trunk/MBCN.cpp 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/MBCN.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -122,21 +122,21 @@
//! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
//! ordering for parent_type
//! \param num_sub_vertices Number of vertices in sub-entity
-//void MBCN::SubEntityConn(const void *parent_conn, const MBEntityType parent_type,
-// const int sub_dimension,
-// const int sub_index,
-// void *sub_entity_conn, int &num_sub_vertices)
-//{
-// static int sub_indices[MB_MAX_SUB_ENTITY_VERTICES];
-//
-// SubEntityVertexIndices(parent_type, sub_dimension, sub_index, sub_indices);
-//
-// num_sub_vertices = VerticesPerEntity(SubEntityType(parent_type, sub_dimension, sub_index));
-// void **parent_conn_ptr = static_cast<void **>(const_cast<void *>(parent_conn));
-// void **sub_conn_ptr = static_cast<void **>(sub_entity_conn);
-// for (int i = 0; i < num_sub_vertices; i++)
-// sub_conn_ptr[i] = parent_conn_ptr[sub_indices[i]];
-//}
+void MBCN::SubEntityConn(const void *parent_conn, const MBEntityType parent_type,
+ const int sub_dimension,
+ const int sub_index,
+ void *sub_entity_conn, int &num_sub_vertices)
+{
+ static int sub_indices[MB_MAX_SUB_ENTITY_VERTICES];
+
+ SubEntityVertexIndices(parent_type, sub_dimension, sub_index, sub_indices);
+
+ num_sub_vertices = VerticesPerEntity(SubEntityType(parent_type, sub_dimension, sub_index));
+ void **parent_conn_ptr = static_cast<void **>(const_cast<void *>(parent_conn));
+ void **sub_conn_ptr = static_cast<void **>(sub_entity_conn);
+ for (int i = 0; i < num_sub_vertices; i++)
+ sub_conn_ptr[i] = parent_conn_ptr[sub_indices[i]];
+}
//! given an entity and a target dimension & side number, get that entity
short int MBCN::AdjacentSubEntities(const MBEntityType this_type,
Modified: MOAB/trunk/MBCN.hpp
===================================================================
--- MOAB/trunk/MBCN.hpp 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/MBCN.hpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -216,10 +216,10 @@
//! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
//! ordering for parent_type
//! \param num_sub_vertices Number of vertices in sub-entity
-// static void SubEntityConn(const void *parent_conn, const MBEntityType parent_type,
-// const int sub_dimension,
-// const int sub_index,
-// void *sub_entity_conn, int &num_sub_vertices);
+ static void SubEntityConn(const void *parent_conn, const MBEntityType parent_type,
+ const int sub_dimension,
+ const int sub_index,
+ void *sub_entity_conn, int &num_sub_vertices);
//! For a specified set of sides of given dimension, return the intersection
//! or union of all sides of specified target dimension adjacent to those sides.
Modified: MOAB/trunk/MBReaderWriterSet.cpp
===================================================================
--- MOAB/trunk/MBReaderWriterSet.cpp 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/MBReaderWriterSet.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -48,6 +48,11 @@
# include "WriteSLAC.hpp"
#endif
+#ifdef CCMIO_FILE
+# include "WriteCCMIO.hpp"
+# include "ReadCCMIO.hpp"
+#endif
+
#ifdef HDF5_FILE
# include "ReadHDF5.hpp"
# ifdef HDF5_PARALLEL
@@ -110,6 +115,12 @@
register_factory( NULL, WriteSLAC::factory, "SLAC", "slac", "SLAC" );
#endif
+#ifdef CCMIO_FILE
+ const char* ccmio_sufxs[] = { "ccm", "ccmg", NULL };
+ register_factory( NULL, WriteCCMIO::factory, "CCMIO files", ccmio_sufxs, "CCMIO");
+ register_factory( ReadCCMIO::factory, NULL, "CCMIO files", ccmio_sufxs, "CCMIO");
+#endif
+
register_factory( NULL, WriteGMV::factory, "GMV", "gmv", "GMV" );
register_factory( NULL, WriteAns::factory, "Ansys", "ans", "ANSYS" );
Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/Makefile.am 2010-01-29 05:24:00 UTC (rev 3508)
@@ -78,6 +78,10 @@
INCLUDES += -I$(srcdir)/mhdf/include
MOAB_EXTRA_SRCS += ReadHDF5.cpp ReadHDF5.hpp WriteHDF5.cpp WriteHDF5.hpp
endif
+if CCMIO_FILE
+ MOAB_EXTRA_SRCS += WriteCCMIO.cpp WriteCCMIO.hpp ReadCCMIO.cpp ReadCCMIO.hpp
+ libMOAB_la_LIBADD += $(CCMIO_LIBS)
+endif
if PARALLEL
libMOAB_la_LIBADD += $(top_builddir)/parallel/libMOABpar.la
Added: MOAB/trunk/ReadCCMIO.cpp
===================================================================
--- MOAB/trunk/ReadCCMIO.cpp (rev 0)
+++ MOAB/trunk/ReadCCMIO.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -0,0 +1,3372 @@
+#include <stdlib.h> // For exit()
+#include <vector>
+#include <map>
+#include <iostream>
+#include <string>
+#include <algorithm>
+
+#include "MBCN.hpp"
+#include "MBRange.hpp"
+#include "MBInterface.hpp"
+#include "MBTagConventions.hpp"
+#include "MBInternals.hpp"
+#include "MBReadUtilIface.hpp"
+#include "FileOptions.hpp"
+#include "ReadCCMIO.hpp"
+#include "MeshTopoUtil.hpp"
+
+#include "ccmio.h"
+
+/*
+ * CCMIO file structure
+ *
+ * Root
+ * State / Problem = "default"
+ * Processor
+ * Vertices
+ * ->CCMIOReadVerticesx
+ * Topology
+ * Cells
+ * Solution
+ * Phase
+ * Field
+ * FieldData
+ * Processor ...
+
+
+ */
+
+enum DataType { kScalar, kVector, kVertex, kCell, kInternalFace, kBoundaryFace,
+ kBoundaryData, kBoundaryFaceData, kCellType };
+
+static int const kNValues = 10; // Number of values of each element to print
+static char const kDefaultState[] = "default";
+static char const kUnitsName[] = "Units";
+static int const kVertOffset = 2;
+static int const kCellInc = 4;
+
+MBReaderIface* ReadCCMIO::factory( MBInterface* iface )
+{ return new ReadCCMIO( iface ); }
+
+ReadCCMIO::ReadCCMIO(MBInterface* impl)
+ : mbImpl(impl)
+{
+ assert(impl != NULL);
+
+ void* ptr = 0;
+ impl->query_interface( "MBReadUtilIface", &ptr );
+ readMeshIface = reinterpret_cast<MBReadUtilIface*>(ptr);
+
+ // initialize in case tag_get_handle fails below
+ mMaterialSetTag = 0;
+ mDirichletSetTag = 0;
+ mNeumannSetTag = 0;
+ mHasMidNodesTag = 0;
+ mGlobalIdTag = 0;
+
+ //! get and cache predefined tag handles
+ int dum_val = 0;
+ MBErrorCode result = impl->tag_get_handle(MATERIAL_SET_TAG_NAME, mMaterialSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(MATERIAL_SET_TAG_NAME,
+ sizeof(int),
+ MB_TAG_SPARSE,
+ MB_TYPE_INTEGER,
+ mMaterialSetTag,
+ &dum_val);
+
+ result = impl->tag_get_handle(DIRICHLET_SET_TAG_NAME, mDirichletSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(DIRICHLET_SET_TAG_NAME,
+ sizeof(int),
+ MB_TAG_SPARSE,
+ MB_TYPE_INTEGER,
+ mDirichletSetTag,
+ &dum_val);
+
+ result = impl->tag_get_handle(NEUMANN_SET_TAG_NAME, mNeumannSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(NEUMANN_SET_TAG_NAME,
+ sizeof(int),
+ MB_TAG_SPARSE,
+ MB_TYPE_INTEGER,
+ mNeumannSetTag,
+ &dum_val);
+
+ result = impl->tag_get_handle(HAS_MID_NODES_TAG_NAME, mHasMidNodesTag);
+ if (MB_TAG_NOT_FOUND == result) {
+ int dum_val_array[] = {0, 0, 0, 0};
+ result = impl->tag_create(HAS_MID_NODES_TAG_NAME,
+ 4*sizeof(int),
+ MB_TAG_SPARSE,
+ MB_TYPE_INTEGER,
+ mHasMidNodesTag,
+ dum_val_array);
+ }
+
+ result = impl->tag_get_handle(GLOBAL_ID_TAG_NAME, mGlobalIdTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_SPARSE,
+ MB_TYPE_INTEGER, mGlobalIdTag, &dum_val);
+
+}
+
+ReadCCMIO::~ReadCCMIO()
+{}
+
+MBErrorCode ReadCCMIO::load_file(const char *file_name,
+ const MBEntityHandle* file_set,
+ const FileOptions& opts,
+ const MBReaderIface::IDTag* subset_list,
+ int subset_list_length,
+ const MBTag* file_id_tag)
+{
+ CCMIOID rootID, problemID, stateID, processorID;
+ CCMIOError error = kCCMIONoErr;
+
+ CCMIOOpenFile(&error, file_name, kCCMIORead, &rootID);
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Problem opening file.");
+ return MB_FAILURE;
+ }
+
+ // get the file state
+ MBErrorCode rval = get_state(rootID, problemID, stateID);
+ if (MB_SUCCESS != rval) return MB_FAILURE;
+
+ // get processors
+ std::set<CCMIOSize_t> procs;
+ rval = get_processors(stateID, processorID, procs);
+ if (MB_SUCCESS != rval) return MB_FAILURE;
+
+ std::set<CCMIOSize_t>::iterator sit;
+ for (sit = procs.begin(); sit != procs.end(); sit++) {
+ rval = read_processor(stateID, processorID, *sit);
+ if (MB_SUCCESS != rval) return MB_FAILURE;
+ }
+
+ return rval;
+}
+
+MBErrorCode ReadCCMIO::read_processor(CCMIOID stateID, CCMIOID processorID, CCMIOSize_t proc)
+{
+ CCMIOError error = kCCMIONoErr;
+ MBErrorCode rval;
+ bool has_solution = true;
+ CCMIOID verticesID, topologyID, solutionID;
+
+ // read the vertices, topology, and solution ids
+ proc = CCMIOSIZEC(0);
+ CCMIONextEntity(&error, stateID, kCCMIOProcessor, &proc, &processorID);
+ CCMIOReadProcessor(&error, processorID, &verticesID, &topologyID, NULL, &solutionID);
+ if(kCCMIONoErr != error) {
+ // might not be a solution; try reading just verts & topology
+ error = kCCMIONoErr;
+ CCMIOReadProcessor(&error, processorID, &verticesID, &topologyID, NULL, NULL);
+ if(kCCMIONoErr == error)
+ hasSolution = false;
+ else {
+ readMeshIface->report_error("Couldn't get vertices and topology.");
+ return MB_FAILURE;
+ }
+ }
+
+ MBRange verts;
+ // vert_map fields: s: none, i: gid, ul: vert handle, r: none
+ //TupleList vert_map(0, 1, 1, 0, 0);
+ TupleList vert_map;
+ rval = read_vertices(proc, processorID, verticesID, topologyID,
+ solutionID, has_solution, verts, vert_map);
+ if (MB_SUCCESS != rval) return rval;
+
+ rval = read_cells(proc, processorID, verticesID, topologyID,
+ solutionID, has_solution, vert_map);
+
+ return rval;
+}
+
+MBErrorCode ReadCCMIO::read_cells(CCMIOSize_t proc, CCMIOID processorID,
+ CCMIOID verticesID, CCMIOID topologyID,
+ CCMIOID solutionID, bool has_solution,
+ TupleList &vert_map)
+{
+
+ CCMIOID cellsID, mapID;
+ CCMIOError error = kCCMIONoErr;
+
+ // get the cells entity and number of cells
+ CCMIOSize_t num_cells;
+ CCMIOGetEntity(&error, topologyID, kCCMIOCells, 0, &cellsID);
+ CCMIOEntitySize(&error, cellsID, &num_cells, NULL);
+
+ // read the cell types and gid map
+ std::vector<int> cell_gids(GETINT32(num_cells)), cell_types(GETINT32(num_cells));
+ CCMIOReadCells(&error, cellsID, &mapID, NULL,
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ // this reads the cellids from the map.
+ CCMIOReadMap(&error, mapID, &cell_gids[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Couldn't read cells or cell id map.");
+ return MB_FAILURE;
+ }
+
+ // read the faces.
+ // face_map fields: s:forward/reverse, i: cell id, ul: face handle, r: none
+#ifdef TUPLE_LIST
+ TupleList face_map(1, 1, 1, 0, 0);
+#else
+ TupleList face_map;
+ SenseList sense_map;
+#endif
+ MBErrorCode rval = read_all_faces(topologyID, vert_map, face_map
+#ifndef TUPLE_LIST
+ , sense_map
+#endif
+ );
+
+ // now construct the cells; sort the face map by cell ids first
+#ifdef TUPLE_LIST
+ rval = face_map.sort(1);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't sort face map by cell id.");
+ return MB_FAILURE;
+ }
+#endif
+ MBRange new_cells;
+ rval = construct_cells(face_map,
+#ifndef TUPLE_LIST
+ sense_map,
+#endif
+ vert_map, new_cells);
+ if (MB_SUCCESS != rval) return rval;
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::construct_cells(TupleList &face_map,
+#ifndef TUPLE_LIST
+ SenseList &sense_map,
+#endif
+ TupleList &vert_map,
+ MBRange &new_cells)
+{
+ std::vector<MBEntityHandle> facehs;
+ std::vector<int> senses;
+ MBEntityHandle cell;
+ MBErrorCode tmp_rval, rval = MB_SUCCESS;
+#ifdef TUPLE_LIST
+ unsigned int i = 0;
+ while (i < face_map.n) {
+ // pull out face handles bounding the same cell
+ facehs.clear();
+ int this_id = face_map.get_int(i);
+ unsigned int inext = i;
+ while (face_map.get_int(inext) == this_id && inext <= face_map.n) {
+ inext++;
+ MBEntityHandle face = face_map.get_ulong(inext);
+ facehs.push_back(face);
+ senses.push_back(face_map.get_short(inext));
+ }
+#else
+ std::map<int,std::vector<MBEntityHandle> >::iterator fmit;
+ std::map<int,std::vector<int> >::iterator smit;
+ for (fmit = face_map.begin(), smit = sense_map.begin();
+ fmit != face_map.end(); fmit++, smit++) {
+
+ // pull out face handles bounding the same cell
+ facehs.clear();
+ int this_id = (*fmit).first;
+ facehs = (*fmit).second;
+ senses.clear();
+ senses = (*smit).second;
+#endif
+ tmp_rval = create_cell_from_faces(facehs, senses, cell);
+ if (MB_SUCCESS != tmp_rval) rval = tmp_rval;
+ else {
+ new_cells.insert(cell);
+ // tag cell with global id
+ tmp_rval = mbImpl->tag_set_data(mGlobalIdTag, &cell, 1, &this_id);
+ if (MB_SUCCESS != tmp_rval) rval = tmp_rval;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::create_cell_from_faces(std::vector<MBEntityHandle> &facehs,
+ std::vector<int> &senses,
+ MBEntityHandle &cell)
+{
+ // test to see if they're one type
+ MBEntityType this_type = mbImpl->type_from_handle(facehs[0]);
+ bool same_type = true;
+ for (std::vector<MBEntityHandle>::iterator vit = facehs.begin(); vit != facehs.end(); vit++) {
+ if (this_type != mbImpl->type_from_handle(*vit)) {
+ same_type = false;
+ break;
+ }
+ }
+
+ // if different, we can quit here, we'll consider this a polyhedron
+ MBErrorCode rval = MB_SUCCESS;
+ if (!same_type ||
+ (MBTRI == this_type && facehs.size() != 4) ||
+ (MBQUAD == this_type && facehs.size() != 6) ||
+ (MBQUAD != this_type && MBTRI != this_type)) {
+ rval = mbImpl->create_element(MBPOLYHEDRON, &facehs[0], facehs.size(), cell);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't make polyhedron.");
+ return MB_FAILURE;
+ }
+ else return MB_SUCCESS;
+ }
+
+ // try tet and hex elements; get connectivity of first face
+ std::vector<MBEntityHandle> verts;
+ rval = mbImpl->get_connectivity(&facehs[0], 1, verts);
+ bool match = false;
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't get connectivity.");
+ return MB_FAILURE;
+ }
+
+ // reverse connectivity if sense is forward, since base face always points
+ // into entity
+ if (senses[0] > 0) std::reverse(verts.begin(), verts.end());
+
+ std::vector<MBEntityHandle> storage;
+ MeshTopoUtil mtu(mbImpl);
+ if (MBTRI == this_type) {
+ // get the 4th vertex through the next tri
+ const MBEntityHandle *conn; int conn_size;
+ rval = mbImpl->get_connectivity(facehs[1], conn, conn_size, true, &storage);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't get connectivity.");
+ return MB_FAILURE;
+ }
+ int i = 0;
+ while (std::find(verts.begin(), verts.end(), conn[i]) != verts.end() && i < conn_size) i++;
+ if (conn_size == i) return MB_FAILURE;
+ match = true;
+ this_type = MBTET;
+ }
+ else if (MBQUAD == this_type) {
+ // build hex from quads
+ // algorithm:
+ // - verts = vertices from 1st quad
+ // - find quad q1 sharing verts[0] and verts[1]
+ // - find quad q2 sharing other 2 verts in q1
+ // - find v1 = opposite vert from verts[1] in q1 , v2 = opposite from verts[0]
+ // - get i = offset of v1 in verts2 of q2, rotate verts2 by i
+ // - if verts2[i+1%4] != v2, flip verts2 by switching verts2[1] and verts2[3]
+ // - append verts2 to verts
+
+
+ // get the other vertices for this hex; need to find the quad with no common vertices
+ MBRange tmp_faces, tmp_verts;
+
+ // get q1, which shares 2 vertices with q0
+ std::copy(facehs.begin(), facehs.end(), mb_range_inserter(tmp_faces));
+ rval = mbImpl->get_adjacencies(&verts[0], 2, 2, false, tmp_faces);
+ if (MB_SUCCESS != rval || tmp_faces.size() != 2) {
+ readMeshIface->report_error("Couldn't get adj face.");
+ return MB_FAILURE;
+ }
+ tmp_faces.erase(facehs[0]);
+ MBEntityHandle q1 = *tmp_faces.begin();
+ // get other 2 verts of q1
+ rval = mbImpl->get_connectivity(&q1, 1, tmp_verts);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't get adj verts.");
+ return MB_FAILURE;
+ }
+ tmp_verts.erase(verts[0]); tmp_verts.erase(verts[1]);
+ // get q2
+ std::copy(facehs.begin(), facehs.end(), mb_range_inserter(tmp_faces));
+ rval = mbImpl->get_adjacencies(tmp_verts, 2, false, tmp_faces);
+ if (MB_SUCCESS != rval || tmp_faces.size() != 2) {
+ readMeshIface->report_error("Couldn't get adj face.");
+ return MB_FAILURE;
+ }
+ tmp_faces.erase(q1);
+ MBEntityHandle q2 = *tmp_faces.begin();
+ // get verts in q2
+ rval = mbImpl->get_connectivity(&q2, 1, storage);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't get adj vertices.");
+ return MB_FAILURE;
+ }
+ // get verts in q1 opposite from v[1] and v[0] in q0
+ MBEntityHandle v0 = 0, v1 = 0;
+ rval = mtu.opposite_entity(q1, verts[1], v0);
+ rval = mtu.opposite_entity(q1, verts[0], v1);
+ if (!v0 || !v1) {
+ readMeshIface->report_error("Trouble finding opposite vertices.");
+ return MB_FAILURE;
+ }
+ // offset of v0 in q2, then rotate and flip
+ unsigned int ioff = std::find(storage.begin(), storage.end(), v0) - storage.begin();
+ if (4 == ioff) {
+ readMeshIface->report_error("Trouble finding offset.");
+ return MB_FAILURE;
+ }
+ if (storage[(ioff+1)%4] != v1) {
+ std::reverse(storage.begin(), storage.end());
+ ioff = std::find(storage.begin(), storage.end(), v0) - storage.begin();
+ }
+ if (0 != ioff)
+ std::rotate(storage.begin(), storage.begin()+ioff, storage.end());
+
+ // copy into verts, and make hex
+ std::copy(storage.begin(), storage.end(), std::back_inserter(verts));
+ match = true;
+ this_type = MBHEX;
+ }
+ if (!match) return MB_FAILURE;
+
+ // now make the element
+ rval = mbImpl->create_element(this_type, &verts[0], verts.size(), cell);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("create_element failed.");
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::read_all_faces(CCMIOID topologyID, TupleList &vert_map,
+ TupleList &face_map
+#ifndef TUPLE_LIST
+ ,SenseList &sense_map
+#endif
+)
+{
+ CCMIOSize_t index = CCMIOSIZEC(0);
+ CCMIOID faceID;
+ MBErrorCode rval;
+
+ // get total # internal/bdy faces, size the face map accordingly
+ int nint_faces = 0, nbdy_faces = 0;
+ CCMIOSize_t nf;
+ CCMIOError error = kCCMIONoErr;
+ while (kCCMIONoErr == CCMIONextEntity(NULL, topologyID, kCCMIOBoundaryFaces, &index,
+ &faceID))
+ {
+ CCMIOEntitySize(&error, faceID, &nf, NULL);
+ nbdy_faces = nbdy_faces + nf;
+ }
+ CCMIOGetEntity(&error, topologyID, kCCMIOInternalFaces, 0, &faceID);
+ CCMIOEntitySize(&error, faceID, &nf, NULL);
+ nint_faces = nint_faces + nf;
+#ifdef TUPLE_LIST
+ face_map.resize(2*nint_faces + nbdy_faces);
+#endif
+
+ // get multiple blocks of bdy faces
+ index = CCMIOSIZEC(0);
+ while (kCCMIONoErr == CCMIONextEntity(NULL, topologyID, kCCMIOBoundaryFaces, &index,
+ &faceID))
+ {
+ rval = read_faces(faceID, kCCMIOBoundaryFaces, vert_map, face_map
+#ifndef TUPLE_LIST
+ , sense_map
+#endif
+ );
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Trouble reading boundary faces.");
+ return MB_FAILURE;
+ }
+ }
+
+ // now get internal faces
+ CCMIOGetEntity(&error, topologyID, kCCMIOInternalFaces, 0, &faceID);
+
+ rval = read_faces(faceID, kCCMIOInternalFaces, vert_map,face_map
+#ifndef TUPLE_LIST
+ , sense_map
+#endif
+ );
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Trouble reading internal faces.");
+ return MB_FAILURE;
+ }
+
+ return rval;
+}
+
+MBErrorCode ReadCCMIO::read_faces(CCMIOID faceID, CCMIOEntity bdy_or_int,
+ TupleList &vert_map,
+ TupleList &face_map
+#ifndef TUPLE_LIST
+ ,SenseList &sense_map
+#endif
+ )
+{
+ if (kCCMIOInternalFaces != bdy_or_int && kCCMIOBoundaryFaces != bdy_or_int) {
+ readMeshIface->report_error("Face type isn't boundary or internal.");
+ return MB_FAILURE;
+ }
+
+ CCMIOSize_t num_faces;
+ CCMIOError error = kCCMIONoErr;
+ CCMIOEntitySize(&error, faceID, &num_faces, NULL);
+
+ // get the size of the face connectivity array (not really a straight connect
+ // array, has n, connect(n), ...)
+ CCMIOSize_t farray_size = CCMIOSIZEC(0);
+ CCMIOReadFaces(&error, faceID, bdy_or_int, NULL, &farray_size, NULL,
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading face connectivity length.");
+ return MB_FAILURE;
+ }
+
+
+ // allocate vectors for holding farray and cells for each face; use new for finer
+ // control of de-allocation
+ int num_sides = (kCCMIOInternalFaces == bdy_or_int ? 2 : 1);
+ int *farray = new int[GETINT32(farray_size)];
+
+ // read farray and make the faces
+ CCMIOID mapID;
+ CCMIOReadFaces(&error, faceID, bdy_or_int, &mapID, NULL,
+ farray, CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading face connectivity.");
+ return MB_FAILURE;
+ }
+
+ std::vector<MBEntityHandle> face_handles(GETINT32(num_faces), 0);
+ MBErrorCode rval = make_faces(farray, vert_map, face_handles);
+ if (MB_SUCCESS != rval) return rval;
+
+ // read face cells and make tuples
+ int *face_cells;
+ if (num_sides*num_faces < farray_size) face_cells = new int[num_sides*GETINT32(num_faces)];
+ else face_cells = farray;
+ CCMIOReadFaceCells(&error, faceID, bdy_or_int, face_cells,
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading face cells.");
+ return MB_FAILURE;
+ }
+
+ int *tmp_ptr = face_cells;
+ for (int i = 0; i < num_faces; i++) {
+#ifdef TUPLE_LIST
+ short forward = 1, reverse = -1;
+ face_map.push_back(&forward, tmp_ptr++, &face_handles[i], NULL);
+ if (2 == num_sides)
+ face_map.push_back(&reverse, tmp_ptr++, &face_handles[i], NULL);
+#else
+ face_map[*tmp_ptr].push_back(face_handles[i]);
+ sense_map[*tmp_ptr++].push_back(1);
+ if (2 == num_sides) {
+ face_map[*tmp_ptr].push_back(face_handles[i]);
+ sense_map[*tmp_ptr++].push_back(-1);
+ }
+#endif
+ }
+
+ // now read & set face gids, reuse face_cells 'cuz we know it's big enough
+ CCMIOReadMap(&error, mapID, face_cells, CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading face gids.");
+ return MB_FAILURE;
+ }
+
+ rval = mbImpl->tag_set_data(mGlobalIdTag, &face_handles[0], face_handles.size(), face_cells);
+ if (MB_SUCCESS != rval) {
+ readMeshIface->report_error("Couldn't set face global ids.");
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::make_faces(int *farray,
+ TupleList &vert_map,
+ std::vector<MBEntityHandle> &new_faces)
+{
+ unsigned int num_faces = new_faces.size();
+ std::vector<MBEntityHandle> verts;
+ MBErrorCode tmp_rval = MB_SUCCESS, rval = MB_SUCCESS;
+
+ for (unsigned int i = 0; i < num_faces; i++) {
+ int num_verts = *farray++;
+ verts.resize(num_verts);
+
+ // fill in connectivity by looking up by gid in vert tuple_list
+ for (int j = 0; j < num_verts; j++) {
+#ifdef TUPLE_LIST
+ int tindex = vert_map.find(1, farray[j]);
+ if (-1 == tindex) {
+ tmp_rval = MB_FAILURE;
+ break;
+ }
+ verts[j] = vert_map.get_ulong(tindex, 0);
+#else
+ verts[j] = (vert_map[farray[j]])[0];
+#endif
+ }
+ farray += num_verts;
+
+ if (MB_SUCCESS == tmp_rval) {
+
+ // make face
+ MBEntityType ftype = (3 == num_verts ? MBTRI :
+ (4 == num_verts ? MBQUAD : MBPOLYGON));
+ MBEntityHandle faceh;
+ tmp_rval = mbImpl->create_element(ftype, &verts[0], num_verts, faceh);
+ if (faceh) new_faces[i] = faceh;
+ }
+
+ if (MB_SUCCESS != tmp_rval) rval = tmp_rval;
+ }
+
+ return rval;
+}
+
+MBErrorCode ReadCCMIO::read_vertices(CCMIOSize_t proc, CCMIOID processorID, CCMIOID verticesID,
+ CCMIOID topologyID, CCMIOID solutionID, bool has_solution,
+ MBRange &verts, TupleList &vert_map)
+{
+ CCMIOError error = kCCMIONoErr;
+
+ // pre-read the number of vertices, so we can pre-allocate & read directly in
+ CCMIOSize_t nverts = CCMIOSIZEC(0);
+ CCMIOEntitySize(&error, verticesID, &nverts, NULL);
+ if(kCCMIONoErr != error) {
+ readMeshIface->report_error("Couldn't get number of vertices.");
+ return MB_FAILURE;
+ }
+
+ // get # dimensions
+ CCMIOSize_t dims;
+ float scale;
+ CCMIOReadVerticesf(&error, verticesID, &dims, NULL, NULL, NULL, CCMIOINDEXC(0), CCMIOINDEXC(1));
+ if(kCCMIONoErr != error) {
+ readMeshIface->report_error("Couldn't get number of dimensions.");
+ return MB_FAILURE;
+ }
+
+ // allocate vertex space
+ MBEntityHandle node_handle = 0;
+ std::vector<double*> arrays;
+ readMeshIface->get_node_arrays(3, GETINT32(nverts), MB_START_ID, node_handle, arrays);
+
+ // read vertex coords
+ CCMIOID mapID;
+ std::vector<double> tmp_coords(GETINT32(dims)*GETINT32(nverts));
+ CCMIOReadVerticesd(&error, verticesID, &dims, &scale, &mapID, &tmp_coords[0],
+ CCMIOINDEXC(0), CCMIOINDEXC(0+nverts));
+ if(kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading vertex coordinates.");
+ return MB_FAILURE;
+ }
+
+ // copy interleaved coords into moab blocked coordinate space
+ int i = 0, threei = 0;
+ for (; i < nverts; i++) {
+ arrays[0][i] = tmp_coords[threei++];
+ arrays[1][i] = tmp_coords[threei++];
+ if (3 == GETINT32(dims)) arrays[2][i] = tmp_coords[threei++];
+ else arrays[2][i] = 0.0;
+ }
+
+ // scale, if necessary
+ if (1.0 != scale) {
+ for(i = 0; i < nverts; i++) {
+ arrays[0][i] *= scale;
+ arrays[1][i] *= scale;
+ if (3 == GETINT32(dims)) arrays[2][i] *= scale;
+ }
+ }
+
+ // put new vertex handles into range
+ verts.insert(node_handle, node_handle+nverts);
+
+ // pack vert_map with global ids and handles for these vertices
+ std::vector<int> gids(GETINT32(nverts));
+ CCMIOReadMap(&error, mapID, &gids[0], CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ if(kCCMIONoErr != error) {
+ readMeshIface->report_error("Trouble reading vertex global ids.");
+ return MB_FAILURE;
+ }
+#ifdef TUPLE_LIST
+ vert_map.resize(GETINT32(nverts));
+ for (i = 0; i < GETINT32(nverts); i++) {
+ vert_map.push_back(NULL, &gids[i], &node_handle, NULL);
+#else
+ for (i = 0; i < GETINT32(nverts); i++) {
+ (vert_map[gids[i]]).push_back(node_handle);
+#endif
+ node_handle += 1;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::get_processors(CCMIOID stateID, CCMIOID &processorID,
+ std::set<CCMIOSize_t> &procs)
+{
+ CCMIOSize_t proc = CCMIOSIZEC(0);
+ while (kCCMIONoErr ==
+ CCMIONextEntity(NULL, stateID, kCCMIOProcessor, &proc, &processorID))
+ procs.insert(proc);
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::get_state(CCMIOID rootID, CCMIOID &problemID, CCMIOID &stateID)
+{
+ CCMIOError error = kCCMIONoErr;
+
+ // first try default
+ CCMIOGetState(&error, rootID, "default", &problemID, &stateID);
+ if (kCCMIONoErr != error) {
+ CCMIOSize_t i = CCMIOSIZEC(0);
+ CCMIOError tmp_error = kCCMIONoErr;
+ CCMIONextEntity(&tmp_error, rootID, kCCMIOState, &i, &stateID);
+ if (kCCMIONoErr == tmp_error)
+ CCMIONextEntity(&error, rootID, kCCMIOProblemDescription,
+ &i, &problemID);
+ }
+ if (kCCMIONoErr != error) {
+ readMeshIface->report_error("Couldn't find state.");
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode ReadCCMIO::read_tag_values( const char* file_name,
+ const char* tag_name,
+ const FileOptions& opts,
+ std::vector<int>& tag_values_out,
+ const IDTag* subset_list,
+ int subset_list_length)
+{
+ return MB_FAILURE;
+}
+
+/*
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetIDsForDomain
+//
+// Purpose:
+// Gets nodes for state, processor, vertices, topology and solution that can
+// be used to query attributes for variables and meshes.
+//
+// Arguments:
+//
+// Returns:
+//
+// Note:
+//
+// Programmer: Brad Whitlock
+// Creation: Mon Aug 6 09:10:29 PDT 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+bool
+avtCCMFileFormat::GetIDsForDomain(int dom,
+ CCMIOID &processor, CCMIOID &vertices, CCMIOID &topology,
+ CCMIOID &solution, bool &hasSolution)
+{
+ const char *mName = "avtCCMFileFormat::GetIDsForDomain: ";
+
+ // Try and get the requested processor.
+ int proc = dom;
+ bool ret = (
+ CCMIONextEntity(NULL, GetState(), kCCMIOProcessor, &proc, &processor) ==
+ kCCMIONoErr);
+ if(ret)
+ {
+ hasSolution = true;
+ ccmErr = kCCMIONoErr;
+ // Try and read the vertices, topology, and solution ids for this
+ // processor.
+ CCMIOReadProcessor(&ccmErr, processor, &vertices, &topology, NULL,
+ &solution);
+ if(ccmErr != kCCMIONoErr)
+ {
+ // That didn't work. (Maybe no solution). See if we can at least
+ // get the vertices and processor.
+ ccmErr = kCCMIONoErr;
+ CCMIOReadProcessor(&ccmErr, processor, &vertices, &topology, NULL,
+ NULL);
+ if(ccmErr == kCCMIONoErr)
+ hasSolution = false;
+ else
+ ret = false;
+ }
+ }
+
+ return ret;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetFaces
+//
+// Purpose: Reads the face info.
+//
+// Arguments:
+// faceID The ID of the face entity.
+// faceType The type of faces (internal or boundary).
+// nFaces How many faces are in the entity.
+// cellIDMap Used to map cell IDs to indices
+// vertexIDMap Used to map vertex IDs to indices
+// minSize Min num verts in a face
+// maxSize Max num verts in a face
+// ci A place to store the face info, must be allocated by
+// calling method.
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 5, 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Removed unused variable.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+// ****************************************************************************
+
+void
+avtCCMFileFormat::GetFaces(CCMIOID faceID, CCMIOEntity faceType,
+ unsigned int nFaces, const IDMap &cellIDMap,
+ const IDMap &vertexIDMap,
+ int &minSize, int &maxSize, CellInfoVector &ci)
+{
+ if (faceType != kCCMIOInternalFaces && faceType != kCCMIOBoundaryFaces)
+ {
+ debug1 << "avtCCMFileFormat::GetFaces encountered an internal error"
+ << endl;
+ return;
+ }
+ int getFacesTimer = visitTimer->StartTimer();
+ CCMIOID mapID;
+ unsigned int nCells = 0, size = 0;
+ //intVector faces;
+ intVector faceNodes, faceCells;
+
+ // Determine the size of the faceNodes array, which is of the
+ // form n1, v1, v2, ...vn1, n2, v1, v2, ... vn2, )
+ CCMIOReadFaces(&ccmErr, faceID, faceType, NULL, &size, NULL,
+ kCCMIOStart, kCCMIOEnd);
+ faceNodes.resize(size);
+ //faces.resize(nFaces);
+ if (faceType == kCCMIOInternalFaces)
+ faceCells.resize(nFaces*2);
+ else
+ faceCells.resize(nFaces);
+ CCMIOReadFaces(&ccmErr, faceID, faceType, &mapID, NULL,
+ &faceNodes[0], kCCMIOStart, kCCMIOEnd);
+ CCMIOReadFaceCells(&ccmErr, faceID, faceType, &faceCells[0],
+ kCCMIOStart, kCCMIOEnd);
+ //CCMIOReadMap(&ccmErr, mapID, &faces[0], kCCMIOStart, kCCMIOEnd);
+
+ unsigned int pos = 0;
+ for (unsigned int i = 0; i < nFaces; ++i)
+ {
+ FaceInfo newFace;
+ //newFace.id = faces[i];
+ if (faceType == kCCMIOInternalFaces)
+ {
+ newFace.cells[0] = faceCells[i*2];
+ newFace.cells[1] = faceCells[i*2+1];
+ }
+ else
+ {
+ newFace.cells[0] = faceCells[i];
+ }
+ int nVerts = faceNodes[pos];
+
+ if (nVerts < minSize)
+ minSize = nVerts;
+ if (nVerts > maxSize)
+ maxSize = nVerts;
+
+ for (unsigned int j = 0; j < nVerts; ++j)
+ {
+ newFace.nodes.push_back( vertexIDMap.IDtoIndex(faceNodes[pos+1+j]) );
+ }
+ // cell ids are 1-origin, so must subract 1 to get the
+ // correct index into the CellInfoVector
+ if (faceType == kCCMIOInternalFaces)
+ {
+ int c0 = cellIDMap.IDtoIndex(newFace.cells[0]);
+ int c1 = cellIDMap.IDtoIndex(newFace.cells[1]);
+
+ ci[c0].faceTypes.push_back(1);
+ ci[c0].faces.push_back(newFace);
+ ci[c1].faceTypes.push_back(2);
+ ci[c1].faces.push_back(newFace);
+ }
+ else
+ {
+ int c0 = cellIDMap.IDtoIndex(newFace.cells[0]);
+
+ ci[c0].faceTypes.push_back(0);
+ ci[c0].faces.push_back(newFace);
+ }
+ pos += faceNodes[pos] +1;
+ }
+ visitTimer->StopTimer(getFacesTimer, "GetFaces");
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::PopulateDatabaseMetaData
+//
+// Purpose:
+// This database meta-data object is like a table of contents for the
+// file. By populating it, you are telling the rest of VisIt what
+// information it can request from you.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Aug 2 15:01:17 PST 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// Fixed determination of centering for Vector variables. Added code to
+// determine if any variables are defined only on parts of the mesh, by
+// reading the number of cells, then mapData size for each var.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
+{
+ const char *mName = "avtCCMFileFormat::PopulateDatabaseMetaData: ";
+
+ // Count the number of processors in the file.
+ // Use that for the number of domains.
+ CCMIOID processor;
+ int proc = 0;
+ int nblocks = 0;
+ while (CCMIONextEntity(NULL, GetState(), kCCMIOProcessor, &proc,
+ &processor) == kCCMIONoErr)
+ {
+ ++nblocks;
+ }
+ debug4 << mName << "Found " << nblocks << " domains in the file." << endl;
+
+#if 0
+ // this will be useful for subsetting by the cell type
+ int nCellTypes = 0;
+ if (CCMIOIsValidEntity(ccmProblem))
+ {
+ // Determine the number of kCCMIOCellTypes present.
+ int i = 0;
+ CCMIOID next;
+ while(CCMIONextEntity(NULL, ccmProblem, kCCMIOCellType, &i, &next)
+ == kCCMIONoErr)
+ ++nCellTypes;
+ }
+#endif
+
+
+#if 0
+ // Read the simulation title.
+ char *title = NULL;
+ ccmErr = CCMIOGetTitle(&ccmErr, GetRoot(), &title);
+ if(title != NULL)
+ {
+ md->SetDatabaseComment(title);
+ free(title);
+ }
+#endif
+
+ for (unsigned int i = 0; i < nblocks; ++i)
+ originalCells.push_back(NULL);
+ // Determine the spatial dimensions.
+ int dims = 3;
+ CCMIOID vertices, topology, solution;
+ bool hasSolution = true;
+ if(GetIDsForDomain(0, processor, vertices, topology, solution,
+ hasSolution))
+ {
+ dims = 1 ;
+ CCMIOReadVerticesf(&ccmErr, vertices, &dims, NULL, NULL, NULL, 0, 1);
+ }
+ else
+ {
+ EXCEPTION1(InvalidFilesException, filenames[0]);
+ }
+
+ // If there's just 1 block, read the mesh now and decompose it into
+ // more meshes, caching them for later.
+ int nDomains = nblocks;
+#ifdef ENABLE_SUBDIVISION
+ if(nblocks == 1)
+ {
+ subdividingSingleMesh = true;
+ md->SetFormatCanDoDomainDecomposition(true);
+ }
+#endif
+
+ // Create a mesh.
+ avtMeshMetaData *mmd = new avtMeshMetaData;
+ mmd->name = "Mesh";
+ mmd->spatialDimension = dims;
+ mmd->topologicalDimension = dims;
+ mmd->meshType = AVT_UNSTRUCTURED_MESH;
+ mmd->numBlocks = nDomains;
+ mmd->cellOrigin = 1;
+ mmd->nodeOrigin = 1;
+ md->Add(mmd);
+
+ // Find variable data
+ if (hasSolution)
+ {
+ // Determine number of cells, in order to find out
+ // which variables are NOT defined on entire mesh.
+ CCMIOError err = kCCMIONoErr;
+ CCMIOID cellsID;
+ CCMIOSize nCells;
+ CCMIOGetEntity(&err, topology, kCCMIOCells, 0, &cellsID);
+ CCMIOEntitySize(&err, cellsID, &nCells, NULL);
+
+ debug5 << mName << "Reading variable information" << endl;
+
+ CCMIOID field, phase = solution;
+ int h = 0;
+ int i = 0;
+ bool oldFile = CCMIONextEntity(NULL, solution, kCCMIOFieldPhase,
+ &h, &phase) != kCCMIONoErr;
+
+ h = 0;
+ while(oldFile ||
+ CCMIONextEntity(NULL, solution, kCCMIOFieldPhase, &h, &phase) ==
+ kCCMIONoErr)
+ {
+ debug5 << mName << "CCMIONextEntity for solution " << solution
+ << " returned h=" << h << ", phase = " << phase << endl;
+
+ if (oldFile)
+ phase = solution;
+ else
+ {
+ int phaseNum = 0;
+ CCMIOGetEntityIndex(NULL, phase, &phaseNum);
+ }
+
+ // This i needs to be set to 0 here because we're asking for the
+ // first field from the current phase. Since we're iterating
+ // over phases so far, we need to reset i to 0 here to get the
+ // first field.
+ while(CCMIONextEntity(NULL, phase, kCCMIOField, &i, &field) ==
+ kCCMIONoErr)
+ {
+ char name[kCCMIOMaxStringLength+1];
+ char sName[kCCMIOProstarShortNameLength+1];
+ char *units = NULL;
+ int usize;
+ CCMIODataType datatype;
+ CCMIODimensionality cdims;
+ CCMIOID fieldData, mapID;
+ CCMIODataLocation type;
+ CCMIOID varField = field;
+
+ debug5 << mName << "CCMIONextEntity for phase " << phase
+ << " returned i=" << i << ", field = " << field << endl;
+
+ // Get the field's name, dims, and data type.
+ CCMIOReadField(&ccmErr, field, name, sName, &cdims, &datatype);
+ debug5 << mName << "CCMIOReadField for field " << field
+ << " returned name=" << name << ", sName = "
+ << sName << endl;
+
+ if (cdims == kCCMIOVector)
+ {
+ CCMIOID scalar;
+ CCMIOError err = kCCMIONoErr;
+ CCMIOReadMultiDimensionalFieldData(&err, field,
+ kCCMIOVectorX, &scalar);
+ if (err == kCCMIONoErr)
+ {
+ // componentized vector, use the X component scalar
+ // in order to retrieve more info about the vector
+ field = scalar;
+ }
+ }
+
+ char *usename;
+ if (oldFile)
+ usename = name;
+ else
+ usename = sName;
+
+ // Read the field's units
+ if (CCMIOReadOptstr(NULL, field, "Units", &usize, NULL) ==
+ kCCMIONoErr)
+ {
+ units = new char[usize+1];
+ CCMIOReadOptstr(&ccmErr, field, "Units", NULL, units);
+ }
+
+ // Reset j to 0 to get the first field data. We read the
+ // field data to determine the variable centering.
+ int j = 0;
+ avtCentering centering = AVT_UNKNOWN_CENT;
+ CCMIOSize cellVarSize = 0;
+
+ while (CCMIONextEntity(NULL, field, kCCMIOFieldData,
+ &j, &fieldData) == kCCMIONoErr)
+ {
+ debug5 << mName << "CCMIONextEntity for field " << field
+ << " returned kCCMIONoErr, j=" << j
+ << ",fieldData=" << fieldData << endl;
+
+ // Read information about the field.
+ CCMIOReadFieldDataf(&ccmErr, fieldData, &mapID, &type,
+ NULL, kCCMIOStart, kCCMIOEnd);
+ if(ccmErr != kCCMIONoErr)
+ {
+ debug5 << mName << "CCMIOReadFieldDataf failed."
+ << endl;
+ continue;
+ }
+
+ // Determine the variable centering.
+ if (type == kCCMIOVertex)
+ {
+ centering = AVT_NODECENT;
+ debug4 << "Var " << usename << " is node centered"
+ << endl;
+ }
+ else if (type == kCCMIOCell)
+ {
+ centering = AVT_ZONECENT;
+ debug4 << "Var " << usename << " is zone centered"
+ << endl;
+ CCMIOSize n;
+ CCMIOIndex max;
+ CCMIOEntitySize(&ccmErr, fieldData, &n, &max);
+ cellVarSize += n;
+ }
+ else if (type == kCCMIOFace)
+ {
+ debug4 << "Var " << usename << " is face-centered, "
+ << "ignoring for now. " << endl;
+ }
+ }
+
+ if (cellVarSize != nCells)
+ {
+ varsOnSubmesh.push_back(usename);
+ }
+ // If we don't have metadata for the variable yet, add it.
+ if(varsToFields.find(usename) == varsToFields.end())
+ {
+ // Variables with an unsupported centering are tagged
+ // invalid.
+ bool validVariable = (centering != AVT_UNKNOWN_CENT);
+
+ if(cdims==kCCMIOScalar)
+ {
+ avtScalarMetaData *smd = new avtScalarMetaData(usename,
+ "Mesh", centering);
+ if (units != NULL)
+ {
+ smd->hasUnits = true;
+ smd->units = units;
+ }
+ smd->validVariable = validVariable;
+ md->Add(smd);
+ varsToFields[usename] = varField;
+ }
+ else if(cdims==kCCMIOVector)
+ {
+ avtVectorMetaData *vmd = new avtVectorMetaData(usename,
+ "Mesh", centering, 3);
+ if (units != NULL)
+ {
+ vmd->hasUnits = true;
+ vmd->units = units;
+ }
+ vmd->validVariable = validVariable;
+ md->Add(vmd);
+ varsToFields[usename] = varField;
+ }
+ else if(cdims==kCCMIOTensor)
+ {
+ avtTensorMetaData *tmd = new avtTensorMetaData(usename,
+ "Mesh", centering, 9);
+ if (units != NULL)
+ {
+ tmd->hasUnits = true;
+ tmd->units = units;
+ }
+#if 1
+ // Just make tensor vars display for now.
+ // Don't plot them.
+ tmd->validVariable = false;
+#else
+ tmd->validVariable = validVariable;
+#endif
+ md->Add(tmd);
+ varsToFields[usename] = varField;
+ }
+ }
+
+ if (units != NULL)
+ {
+ delete [] units;
+ units = NULL;
+ }
+ }
+ oldFile = false;
+ }
+ }
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetMesh
+//
+// Purpose:
+// Gets the mesh associated with this file. The mesh is returned as a
+// derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
+// vtkUnstructuredGrid, etc).
+//
+// Arguments:
+// domain The index of the domain. If there are NDomains, this
+// value is guaranteed to be between 0 and NDomains-1,
+// regardless of block origin.
+// meshname The name of the mesh of interest. This can be ignored if
+// there is only one mesh.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Aug 2 15:01:17 PST 2007
+//
+// Modifications:
+// Brad Whitlock, Tue Dec 18 12:57:28 PST 2007
+// Added support for 2D polygonal shapes.
+//
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// If the primary variable (activeVisItVar) is defined only on a portion of
+// the mesh, only retrieve (and tesselate) those cells it is defined upon.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// Brad Whitlock, Thu Oct 1 13:36:48 PDT 2009
+// I refactored this routine into helper routines and added support for
+// automatically subdividing a mesh on the fly.
+//
+// ****************************************************************************
+
+#ifndef MDSERVER
+#include <PolygonToTriangles.C>
+#endif
+
+vtkDataSet *
+avtCCMFileFormat::GetMesh(int domain, const char *meshname)
+{
+#ifdef MDSERVER
+ return 0;
+#else
+ // Override domain if we're automatically dividing the data
+ int dom = subdividingSingleMesh ? 0 : domain;
+
+ vtkUnstructuredGrid *ugrid = NULL;
+ vtkPoints *points = NULL;
+
+ TRY
+ {
+ // Read the points
+ points = ReadPoints(dom, meshname);
+
+ // Read the cell connectivity
+ CellInfoVector cellInfo;
+ int minFaceSize = VTK_LARGE_INTEGER;
+ int maxFaceSize = -1;
+ ReadCellInfo(dom, meshname,
+ cellInfo, minFaceSize, maxFaceSize);
+
+ //
+ // Convert cellInfo into vtkUnstructuredGrid
+ //
+ SelectCellsForThisProcessor(cellInfo, points);
+
+ ugrid = vtkUnstructuredGrid::New();
+
+ // Determine cell topology from face lists
+ if (minFaceSize == 2 && maxFaceSize == 2)
+ {
+ // 2D edges that we must assemble into polygons and tessellate into
+ // triangles that VisIt can digest.
+ TesselateCells2D(domain, cellInfo, points, ugrid);
+ }
+ else if (minFaceSize <= 4 && maxFaceSize <= 4)
+ {
+#ifdef ENABLE_SUBDIVISION
+ // If we're subdividing a single domain on the fly then we create
+ // original cell numbers so we can use them in the GetVar method
+ // to return only the cell values that we selected for this chunk
+ // of the mesh.
+ unsigned int oc[2] = {dom, 0};
+ vtkUnsignedIntArray *origCells = 0;
+ if(subdividingSingleMesh)
+ {
+ int useCount = 0;
+ for(unsigned int i = 0; i < cellInfo.size(); ++i)
+ useCount += (cellInfo[i].id != -1) ? 1 : 0;
+ origCells = vtkUnsignedIntArray::New();
+ origCells->SetName("avtOriginalCellNumbers");
+ origCells->SetNumberOfComponents(2);
+ origCells->Allocate(useCount * 3);
+ originalCells[dom] = origCells;
+ }
+#endif
+ ugrid->SetPoints(points);
+ // We have zoo elements that we can deal with.
+ vtkCellArray *cellArray = vtkCellArray::New();
+ intVector cellTypes;
+ bool unhandledCellType = false;
+ for (unsigned int i = 0; i < cellInfo.size(); i++)
+ {
+ const CellInfo &ci = cellInfo[i];
+ if(ci.id == -1)
+ continue;
+#ifdef ENABLE_SUBDIVISION
+ oc[1] = ci.id;
+#endif
+ switch(ci.faces.size())
+ {
+ case 4 :
+ BuildTet(ci, cellArray, cellTypes);
+#ifdef ENABLE_SUBDIVISION
+ if(subdividingSingleMesh)
+ origCells->InsertNextTupleValue(oc);
+#endif
+ break;
+ case 5 :
+ {
+ int nNodes = 0;
+ for (size_t j = 0; j < ci.faces.size(); j++)
+ {
+ nNodes += ci.faces[j].nodes.size();
+ }
+ if (nNodes == 16) // Pyramid
+ {
+ BuildPyramid(ci, cellArray, cellTypes);
+#ifdef ENABLE_SUBDIVISION
+ if(subdividingSingleMesh)
+ origCells->InsertNextTupleValue(oc);
+#endif
+ }
+ else if (nNodes == 18) // Wedge
+ {
+ BuildWedge(ci, cellArray, cellTypes);
+#ifdef ENABLE_SUBDIVISION
+ if(subdividingSingleMesh)
+ origCells->InsertNextTupleValue(oc);
+#endif
+ }
+ else
+ unhandledCellType = true;
+ break;
+ }
+ case 6 :
+ BuildHex(ci, cellArray, cellTypes);
+#ifdef ENABLE_SUBDIVISION
+ if(subdividingSingleMesh)
+ origCells->InsertNextTupleValue(oc);
+#endif
+ break;
+ default :
+ unhandledCellType = true;
+ break;
+ }
+ }
+ ugrid->SetCells(&cellTypes[0], cellArray);
+ cellArray->Delete();
+ }
+ else
+ {
+ TesselateCell(domain, cellInfo, points, ugrid);
+ }
+
+ points->Delete();
+ }
+ CATCHALL
+ {
+ if(points != NULL)
+ points->Delete();
+ if(ugrid != NULL)
+ ugrid->Delete();
+ }
+ ENDTRY
+
+ return ugrid;
+#endif
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::ReadPoints
+//
+// Purpose:
+// Reads the points associated with a specified domain/mesh.
+//
+// Arguments:
+// dom : The domain number.
+// meshname : The name of the mesh.(for error reporting)
+//
+// Returns: The vtkPoints object that contains the points.
+//
+// Note: This method was broken out from the original GetMesh routine.
+//
+// Programmer: Kathleen Bonnell, Brad Whitlock
+// Creation: Thu Oct 1 13:29:03 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+vtkPoints *
+avtCCMFileFormat::ReadPoints(int dom, const char *meshname)
+{
+ CCMIOID processor, vertices, topology, solution;
+ bool hasSolution = true;
+ if(!GetIDsForDomain(dom, processor, vertices, topology, solution,
+ hasSolution))
+ {
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Read the size of the vertices
+ CCMIOSize nnodes = 0;
+ CCMIOEntitySize(&ccmErr, vertices, &nnodes, NULL);
+ if(ccmErr != kCCMIONoErr)
+ {
+ debug4 << "CCMIOEntitySize for vertices failed with error " ;
+ debug4 << ccmErr << endl;
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Read the dimensions of the vertex.
+ int dims = 1;
+ float scale;
+ CCMIOReadVerticesf(&ccmErr, vertices, &dims, NULL, NULL, NULL, 0, 1);
+ if(ccmErr != kCCMIONoErr)
+ {
+ debug4 << "CCMIOReadVertices for first vertex dimensions ";
+ debug4 << "failed with error " << ccmErr << endl;
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Allocate VTK memory.
+ vtkPoints *points = vtkPoints::New();
+ points->SetNumberOfPoints(nnodes);
+ float *pts = (float *)points->GetVoidPointer(0);
+
+ // Read the data into the VTK points.
+ CCMIOID mapID;
+ if(dims == 2)
+ {
+ // Read 2D points and convert to 3D, storing into VTK.
+ float *pts2d = new float[2 * nnodes];
+ CCMIOReadVerticesf(&ccmErr, vertices, &dims, &scale, &mapID, pts2d,
+ 0, nnodes);
+ float *src = pts2d;
+ float *dest = pts;
+ for(int i = 0; i < nnodes; ++i)
+ {
+ *dest++ = *src++;
+ *dest++ = *src++;
+ *dest++ = 0.;
+ }
+ delete [] pts2d;
+ }
+ else
+ {
+ // Read the data directly into the VTK buffer.
+ CCMIOReadVerticesf(&ccmErr, vertices, &dims, &scale, &mapID, pts,
+ 0, nnodes);
+ }
+
+ // Scale the points, according to the scale factor read with the
+ // vertices.
+ for(int i = 0; i < nnodes; ++i)
+ {
+ pts[0] *= scale;
+ pts[1] *= scale;
+ pts[2] *= scale;
+ pts += 3;
+ }
+
+ return points;
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::ReadCellInfo
+//
+// Purpose:
+// Reads the cell info associated with a specified domain/mesh.
+//
+// Arguments:
+// dom : The domain number.
+// meshname : The name of the mesh.(for error reporting)
+// cellInfo : The cell information to populate.
+// minFaceSize : Return value of the min face size.
+// maxFaceSize : Return value of the max face size.
+//
+// Returns:
+//
+// Note: This method was broken out from the original GetMesh routine.
+// min/maxFaceSize are used to determine whether tesselation is
+// required.
+//
+// Programmer: Kathleen Bonnell, Brad Whitlock
+// Creation: Thu Oct 1 13:29:03 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::ReadCellInfo(int dom, const char *meshname,
+ CellInfoVector &cellInfo, int &minFaceSize, int &maxFaceSize)
+{
+ CCMIOID processor, vertices, topology, solution;
+ bool hasSolution = true;
+ if(!GetIDsForDomain(dom, processor, vertices, topology, solution,
+ hasSolution))
+ {
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Read the size of the vertices
+ CCMIOSize nnodes = 0;
+ CCMIOEntitySize(&ccmErr, vertices, &nnodes, NULL);
+ if(ccmErr != kCCMIONoErr)
+ {
+ debug4 << "CCMIOEntitySize for vertices failed with error " ;
+ debug4 << ccmErr << endl;
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Read the dimensions of the vertex and get the mapID
+ int dims = 1;
+ CCMIOID mapID;
+ CCMIOReadVerticesf(&ccmErr, vertices, &dims, NULL, &mapID, NULL, 0, 1);
+ if(ccmErr != kCCMIONoErr)
+ {
+ debug4 << "CCMIOReadVertices for first vertex dimensions ";
+ debug4 << "failed with error " << ccmErr << endl;
+ EXCEPTION1(InvalidVariableException, meshname);
+ }
+
+ // Read the vertex ids
+ intVector tmpVertexMap(nnodes);
+ CCMIOReadMap(&ccmErr, mapID, &tmpVertexMap[0], kCCMIOStart, kCCMIOEnd);
+ IDMap vertexIDMap;
+ vertexIDMap.SetIDs(tmpVertexMap);
+ tmpVertexMap.clear();
+
+ // Get the topology information
+ CCMIOID faceID, cellsID;
+ unsigned int nIFaces = 0, nCells = 0, size = 0;
+ intVector cells;
+ //intVector cellMatType;
+
+ // Read the cells entity
+ CCMIOGetEntity(&ccmErr, topology, kCCMIOCells, 0, &cellsID);
+ // Read the cells entity size (num cells)
+ CCMIOEntitySize(&ccmErr, cellsID, &nCells, NULL);
+ cells.resize(nCells);
+ cellInfo.resize(nCells);
+ //cellMatType.resize(nCells);
+ // this gets the cell types and the map that stores the cell ids
+ CCMIOReadCells(&ccmErr, cellsID, &mapID, NULL,
+ kCCMIOStart, kCCMIOEnd);
+ // this reads the cellids from the map.
+ CCMIOReadMap(&ccmErr, mapID, &cells[0], kCCMIOStart, kCCMIOEnd);
+ for (int i = 0; i < nCells; ++i)
+ cellInfo[i].id = cells[i];
+
+ IDMap cellIDMap;
+ cellIDMap.SetIDs(cells);
+
+ // Read the boundary faces.
+ int index = 0;
+ int count = 0;
+ int nBoundaries = 0;
+ while (CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index,
+ &faceID) == kCCMIONoErr)
+ {
+ nBoundaries++;
+ }
+
+ index = 0;
+ while (CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index,
+ &faceID) == kCCMIONoErr)
+ {
+ CCMIOSize nBFaces;
+ CCMIOEntitySize(&ccmErr, faceID, &nBFaces, NULL);
+ GetFaces(faceID, kCCMIOBoundaryFaces, nBFaces,
+ cellIDMap, vertexIDMap,
+ minFaceSize, maxFaceSize, cellInfo);
+ }
+
+ // Read the internal faces.
+ // Get the InternalFaces entity.
+ CCMIOGetEntity(&ccmErr, topology, kCCMIOInternalFaces, 0, &faceID);
+ // Get the InternalFaces size (num faces).
+ CCMIOEntitySize(&ccmErr, faceID, &nIFaces, NULL);
+
+ GetFaces(faceID, kCCMIOInternalFaces, nIFaces, cellIDMap, vertexIDMap,
+ minFaceSize, maxFaceSize, cellInfo);
+
+
+ if (find(varsOnSubmesh.begin(), varsOnSubmesh.end(), activeVisItVar)
+ != varsOnSubmesh.end())
+ {
+ // need to reduce the number of cells we actually process.
+ intVector validCells;
+ CellInfoVector vcv;
+ GetCellMapData(dom, activeVisItVar, validCells);
+
+ for (int i = 0; i < cellInfo.size(); ++i)
+ {
+ if (find(validCells.begin(), validCells.end(), cellInfo[i].id)
+ != validCells.end())
+ {
+ vcv.push_back(cellInfo[i]);
+ }
+ }
+ cellInfo = vcv;
+ }
+
+ debug5 << "minFaceSize = " << minFaceSize
+ << ", maxFaceSize = " << maxFaceSize << endl;
+}
+
+#ifdef ENABLE_SUBDIVISION
+// ****************************************************************************
+// Method: ComputePatchCenter
+//
+// Purpose:
+// Computes the center and bounds of the patch.
+//
+// Arguments:
+// centers : The centers for each cell in the mesh.
+// patch : The list of cellids that make up the patch.
+// center : The calculated center of the patch.
+// bounds : The calculated bounds of the patch.
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Oct 6 16:12:08 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+static void
+ComputePatchCenter(const double *centers, const intVector &patch, double *center, double *bounds)
+{
+ center[0] = 0.;
+ center[1] = 0.;
+ center[2] = 0.;
+ int nMatches = 0;
+ for(int i = 0; i < patch.size(); ++i)
+ {
+ const double *c = centers + patch[i] * 3;
+
+ // Compute extents of cell centers.
+ if(nMatches == 0)
+ {
+ bounds[0] = bounds[1] = c[0];
+ bounds[2] = bounds[3] = c[1];
+ bounds[4] = bounds[5] = c[2];
+ }
+ else
+ {
+ if(c[0] < bounds[0])
+ bounds[0] = c[0];
+ if(c[0] > bounds[1])
+ bounds[1] = c[0];
+
+ if(c[1] < bounds[2])
+ bounds[2] = c[1];
+ if(c[1] > bounds[3])
+ bounds[3] = c[1];
+
+ if(c[2] < bounds[4])
+ bounds[4] = c[2];
+ if(c[2] > bounds[5])
+ bounds[5] = c[2];
+ }
+
+ center[0] += c[0];
+ center[1] += c[1];
+ center[2] += c[2];
+ nMatches++;
+ }
+ if(nMatches > 0)
+ {
+ center[0] /= double(nMatches);
+ center[1] /= double(nMatches);
+ center[2] /= double(nMatches);
+ }
+}
+
+// ****************************************************************************
+// Method: DivideLargestPatch
+//
+// Purpose:
+// Divides the largest patch in the patch vector.
+//
+// Arguments:
+// centers : The centers of all cells in the mesh.
+// patches : The list of all patches.
+//
+// Returns:
+//
+// Note: This routine modifies the patches vector by splitting 1 of the
+// patches, replacing the split patch with piece0. piece1 is appended
+// to the patch vector.
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Oct 6 16:13:33 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+static void
+DivideLargestPatch(const double *centers, std::vector<intVector> &patches)
+{
+ // Find the index of the largest patch
+ int maxIndex = 0;
+ for(int i = 1; i < patches.size(); ++i)
+ if(patches[i].size() > patches[maxIndex].size())
+ maxIndex = i;
+
+ // Compute the center at which we will bisect.
+ double center[3], bounds[6];
+ ComputePatchCenter(centers, patches[maxIndex], center, bounds);
+
+ // Figure out the longest dimension since that's the dimension we'll bisect.
+ double dX = bounds[1] - bounds[0];
+ double dY = bounds[3] - bounds[2];
+ double dZ = bounds[5] - bounds[4];
+ int longestDimension = 2;
+ if(dX > dY)
+ {
+ if(dX > dZ)
+ longestDimension = 0;
+ }
+ else
+ {
+ if(dY > dZ)
+ longestDimension = 1;
+ }
+
+ const intVector &patch = patches[maxIndex];
+ intVector piece0, piece1;
+ for(int j = 0; j < patch.size(); ++j)
+ {
+ const double *c = centers + patch[j] * 3;
+ if(c[longestDimension] > center[longestDimension])
+ piece0.push_back(patch[j]);
+ else
+ piece1.push_back(patch[j]);
+ }
+ patches[maxIndex] = piece0;
+ patches.push_back(piece1);
+}
+#endif
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::SelectCellsForThisProcessor
+//
+// Purpose:
+// This routine divides the cells spatially into PAR_Size() different bins
+// and sets all of the ids for the cells in cellInfo to -1 (invalid) unless
+// their bin matches PAR_Rank(). This means that we are marking a subset of
+// the cells in cellInfo as being valid so we can return just a part of the
+// dataset.
+//
+// Arguments:
+// cellInfo : The vector of cell data.
+// points : The points used by the cells.
+//
+// Returns:
+//
+// Note: This method is just used in parallel when we have a single
+// domain dataset that we want to automatically chunk up under the
+// covers.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Oct 1 13:25:17 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::SelectCellsForThisProcessor(CellInfoVector &cellInfo, vtkPoints *points)
+{
+#ifdef ENABLE_SUBDIVISION
+ if(subdividingSingleMesh)
+ {
+ // Compute cell centers
+ double *centers = new double[cellInfo.size() * 3];
+ for(size_t i = 0; i < cellInfo.size(); ++i)
+ cellInfo[i].CellCenter(centers + i * 3, points);
+
+ // Start out with all cells in 1 patch
+ std::vector<intVector> patches;
+ intVector allCells;
+ for(size_t i = 0; i < cellInfo.size(); ++i)
+ allCells.push_back(i);
+ patches.push_back(allCells);
+
+ // Divide the largest patch until we have enough patches.
+ while(patches.size() < PAR_Size())
+ DivideLargestPatch(centers, patches);
+
+ // Set cellid to -1 unless we're on the patch whose id == PAR_Rank.
+ for(size_t p = 0; p < patches.size(); ++p)
+ {
+ if(p == PAR_Rank())
+ continue;
+
+ const intVector &patch = patches[p];
+ for(size_t i = 0; i < patch.size(); ++i)
+ cellInfo[patch[i]].id = -1;
+ }
+
+ delete [] centers;
+ }
+#endif
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::RegisterVariableList
+//
+// Purpose:
+// Records the active variable name so only cells valid for the var
+// will be retrieved during GetMesh calls.
+//
+// Programmer: Kathleen Bonnell
+// Creation: February 28, 2008
+//
+// ***************************************************************************
+
+void
+avtCCMFileFormat::RegisterVariableList(const char *primaryVar,
+ const std::vector<CharStrRef> &)
+{
+ activeVisItVar = primaryVar;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetCellMapData
+//
+// Purpose:
+// Retrieves the map data for a given var, which specifies the cell
+// ids the var is defined upon.
+//
+// Arguments:
+// domain The required domain.
+// var The requested variable name.
+// mapData A place to store the retrieved map data.
+//
+// Programmer: Kathleen Bonnell
+// Creation: February 28, 2008
+//
+// ***************************************************************************
+
+void
+avtCCMFileFormat::GetCellMapData(const int domain, const string &var,
+ intVector &mapData)
+{
+ VarFieldMap::const_iterator pos = varsToFields.find(var.c_str());
+ if (pos == varsToFields.end())
+ EXCEPTION1(InvalidVariableException, var);
+
+ mapData.clear();
+
+ CCMIOSize n;
+ CCMIOIndex fmax;
+ CCMIOError ferr = kCCMIONoErr;
+ int j = 0, cnt = 0;
+ CCMIOID fieldData, mapID;
+ CCMIODimensionality cdims;
+ CCMIOID field = pos->second;
+
+ CCMIOReadField(&ferr, field, NULL, NULL, &cdims, NULL);
+ if (cdims == kCCMIOVector)
+ {
+ CCMIOID scalar;
+ ferr = kCCMIONoErr;
+ CCMIOReadMultiDimensionalFieldData(&ferr, field,
+ kCCMIOVectorX, &scalar);
+ if (ferr == kCCMIONoErr)
+ {
+ // componentized vector, use the X component scalar
+ // in order to retrieve more info about the vector
+ field = scalar;
+ }
+ }
+ ferr = kCCMIONoErr;
+ while(CCMIONextEntity(NULL, field, kCCMIOFieldData, &j, &fieldData)
+ == kCCMIONoErr)
+ {
+ CCMIOEntitySize(&ferr, fieldData, &n, &fmax);
+ CCMIOReadFieldDataf(&ferr, fieldData, &mapID, NULL, NULL,
+ kCCMIOStart, kCCMIOEnd);
+ mapData.resize(cnt+n);
+ CCMIOReadMap(&ferr, mapID, &mapData[cnt], kCCMIOStart, kCCMIOEnd);
+ cnt += n;
+ }
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetVar
+//
+// Purpose:
+// Gets a scalar variable associated with this file. Although VTK has
+// support for many different types, the best bet is vtkFloatArray, since
+// that is supported everywhere through VisIt.
+//
+// Arguments:
+// domain The index of the domain. If there are NDomains, this
+// value is guaranteed to be between 0 and NDomains-1,
+// regardless of block origin.
+// varname The name of the variable requested.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Aug 2 15:01:17 PST 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// avtOriginalCellNumbers array now contains CCM cellid, so ensure that
+// the data array is indexed-into correctly, by mapping the original cell
+//
+// Brad Whitlock, Thu Oct 1 13:08:35 PDT 2009
+// Set domain to 0 if we're subdividing a single mesh.
+//
+// ****************************************************************************
+
+vtkDataArray *
+avtCCMFileFormat::GetVar(int domain, const char *varname)
+{
+ const char *mName = "avtCCMFileFormat::GetVar: ";
+ // Override domain if we're automatically dividing the data
+ domain = subdividingSingleMesh ? 0 : domain;
+
+ VarFieldMap::const_iterator pos = varsToFields.find(varname);
+ if (pos == varsToFields.end())
+ EXCEPTION1(InvalidVariableException, varname);
+
+ intVector mapData;
+ floatVector data;
+ ReadScalar(pos->second, mapData, data);
+
+ if (ccmErr != kCCMIONoErr)
+ EXCEPTION1(InvalidVariableException, varname);
+
+ unsigned int nvalues = data.size();
+
+ vtkFloatArray *rv = vtkFloatArray::New();
+ if (originalCells[domain] == NULL)
+ {
+ rv->SetNumberOfValues(nvalues);
+ for (unsigned int i = 0 ; i < nvalues ; ++i)
+ {
+ rv->SetValue(i, data[i]);
+ }
+ }
+ else
+ {
+ // We've tesselated, and have more cells than nvalues. Need to
+ // duplicate values for all cells that share the same original
+ // cell number!
+ vtkUnsignedIntArray *ocarray =
+ vtkUnsignedIntArray::SafeDownCast(originalCells[domain]);
+ int numCells = ocarray->GetNumberOfTuples();
+
+ debug4 << mName << "numCells = " << numCells << endl;
+ unsigned int *oc = ocarray->GetPointer(0);
+ rv->SetNumberOfValues(numCells);
+ std::map<int, int> cellIdMap;
+ for (unsigned int i = 0; i < mapData.size(); i++)
+ {
+ cellIdMap[mapData[i]] = i;
+ }
+
+ for (unsigned int i = 0 ; i < numCells ; i++)
+ {
+ rv->SetValue(i, data[cellIdMap[oc[i*2+1]]]);
+ }
+ }
+
+ return rv;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::GetVectorVar
+//
+// Purpose:
+// Gets a vector variable associated with this file. Although VTK has
+// support for many different types, the best bet is vtkFloatArray, since
+// that is supported everywhere through VisIt.
+//
+// Arguments:
+// domain The index of the domain. If there are NDomains, this
+// value is guaranteed to be between 0 and NDomains-1,
+// regardless of block origin.
+// varname The name of the variable requested.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Aug 2 15:01:17 PST 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// avtOriginalCellNumbers array now contains CCM cellid, so enusre that
+// the data array is indexed-into correctly, by mapping the original cell
+// id through an cellIdMap of valid cell ids (obtained from MapData).
+//
+// Dave Bremer, Fri Apr 11 16:49:45 PDT 2008
+// Initialize the err variable.
+//
+// Brad Whitlock, Thu Oct 1 13:08:35 PDT 2009
+// Set domain to 0 if we're subdividing a single mesh.
+//
+// ****************************************************************************
+
+vtkDataArray *
+avtCCMFileFormat::GetVectorVar(int domain, const char *varname)
+{
+ // Override domain if we're automatically dividing the data
+ domain = subdividingSingleMesh ? 0 : domain;
+
+ VarFieldMap::const_iterator pos = varsToFields.find(varname);
+ if (pos == varsToFields.end())
+ EXCEPTION1(InvalidVariableException, varname);
+
+ intVector mapData;
+ floatVector u, v, w, data;
+ CCMIOID scalar;
+ CCMIOID field = pos->second;
+ CCMIOError err = kCCMIONoErr;
+
+ CCMIOReadMultiDimensionalFieldData(&err, field, kCCMIOVectorX, &scalar);
+
+ if (err == kCCMIOVersionErr)
+ {
+ // If we are reading an older version of the file,
+ // where vectors are stored as vectors, not components,
+ // we need to call CCMIOReadFieldData*(), which is
+ // all that ReadScalar() does.
+ err = kCCMIONoErr;
+ ReadScalar(field, mapData, data, true);
+ }
+ else
+ {
+ ReadScalar(scalar, mapData, u);
+ CCMIOReadMultiDimensionalFieldData(&err, field, kCCMIOVectorY, &scalar);
+ ReadScalar(scalar, mapData, v);
+ CCMIOReadMultiDimensionalFieldData(&err, field, kCCMIOVectorZ, &scalar);
+ ReadScalar(scalar, mapData, w);
+ data.resize(3 * u.size());
+ for (unsigned int k = 0; k < u.size(); ++k)
+ {
+ data[3 * k ] = u[k];
+ data[3 * k + 1] = v[k];
+ data[3 * k + 2] = w[k];
+ }
+ }
+
+ vtkFloatArray *rv = vtkFloatArray::New();
+ rv->SetNumberOfComponents(3);
+ if (originalCells[domain] == NULL)
+ {
+ unsigned int nvalues = data.size();
+ rv->SetNumberOfTuples(nvalues/3);
+ float *v = rv->WritePointer(0, nvalues);
+ for (unsigned int i = 0 ; i < nvalues ; i++)
+ {
+ v[i] = data[i];
+ }
+ }
+ else
+ {
+ // We've tesselated, and have more cells than nvalues. Need to
+ // duplicate values for all cells that share the same original
+ // cell number!
+ vtkUnsignedIntArray *ocarray =
+ vtkUnsignedIntArray::SafeDownCast(originalCells[domain]);
+ int numCells = ocarray->GetNumberOfTuples();
+
+ unsigned int *oc = ocarray->GetPointer(0);
+ rv->SetNumberOfTuples(numCells);
+ float *v = rv->WritePointer(0, numCells*3);
+
+ std::map<int, int> cellIdMap;
+ for (unsigned int i = 0; i < mapData.size(); ++i)
+ {
+ cellIdMap[mapData[i]] = i;
+ }
+
+ for (unsigned int i = 0 ; i < numCells; ++i)
+ {
+ unsigned int id = cellIdMap[oc[i*2+1]];
+ v[i*3+0] = data[id*3+0];
+ v[i*3+1] = data[id*3+1];
+ v[i*3+2] = data[id*3+2];
+ }
+ }
+ return rv;
+}
+
+
+// ***************************************************************************
+// Method: avtCCMFileFormat::ReadScalar
+//
+// Purpose:
+// Reads scalar data from the file.
+//
+// Arguments:
+// field The ID of the field to read.
+// mapData A place to store data specifying which cells the scalar
+// is defined upon.
+// data A place to store the scalar data.
+// readingVector Indicates if an old-style vector is being read.
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// Loop through field entities, as a var may be delineated by cell types,
+// thus being represented by multiple 'fields'. Correctly resize mapData
+// and data each time through the loop. Combined code for Cell and Node
+// data as they were exactly the same.
+//
+// ***************************************************************************
+
+void
+avtCCMFileFormat::ReadScalar(CCMIOID field, intVector &mapData,
+ floatVector &data, bool readingVector)
+
+{
+ const char *mName = "avtCCMFileFormat::ReadScalar: ";
+ CCMIOSize n;
+ CCMIOIndex fmax;
+ int j = 0, cnt = 0;
+ CCMIOID fieldData, mapID;
+ CCMIODataLocation type;
+
+ mapData.clear();
+ data.clear();
+ // Read each piece of field data
+ // The fields may be delineated by cell types, so there may be
+ // multiple fields for the same variable.
+ while (CCMIONextEntity(NULL, field, kCCMIOFieldData, &j, &fieldData)
+ == kCCMIONoErr)
+ {
+ // Figure out how big this data is so we can read it. If we were
+ // storing this information permanently we might use a sparse
+ // array, in which case we would need to find the maximum ID and
+ // make the array that size.
+ CCMIOEntitySize(&ccmErr, fieldData, &n, &fmax);
+ CCMIOReadFieldDataf(&ccmErr, fieldData, &mapID, &type, NULL,
+ kCCMIOStart, kCCMIOEnd);
+
+ if (type == kCCMIOCell || type == kCCMIOVertex)
+ {
+ mapData.resize(cnt +n);
+ CCMIOReadMap(&ccmErr, mapID, &mapData[cnt], kCCMIOStart, kCCMIOEnd);
+ if (readingVector)
+ data.resize(cnt +(3 * n));
+ else
+ data.resize(cnt + n);
+ if (type == kCCMIOCell)
+ debug4 << mName << "Reading cell data n= " << n << endl;
+ else
+ debug4 << mName << "Reading node data n= " << n << endl;
+ // If we want double precision we should use
+ // CCMIOReadFieldDatad().
+ CCMIOReadFieldDataf(&ccmErr, fieldData, &mapID, NULL,
+ &data[cnt], kCCMIOStart, kCCMIOEnd);
+ cnt += n;
+ }
+#if 0
+ else if (type == kCCMIOFace)
+ {
+ debug3 << "\tReadScalar found type kCCMIOFace" << endl;
+ }
+#endif
+
+ if (ccmErr != kCCMIONoErr)
+ debug1 << " Error reading scalar data " << ccmErr << endl;
+ }
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::TesselateCell
+//
+// Purpose:
+//
+//
+//
+//
+// Arguments:
+// civ Contains cell information.
+// points The points comprising the dataset.
+// ugrid The unstructured grid we are building.
+//
+//
+// Programmer: Kathleen Bonnell
+// Creation: October 1, 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// Use cellid stored in CellInfo in avtOriginalCellNumbers array.
+//
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Change fbounds to doubleVector to get around compiler problem on Windows.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::TesselateCell(const int domain, const CellInfoVector &civ,
+ vtkPoints *points, vtkUnstructuredGrid *ugrid)
+{
+#ifndef MDSERVER
+ int dom = subdividingSingleMesh ? 0 : domain;
+
+ const char *mName = "avtCCMFileFormat::TesselateCell: ";
+ unsigned int i, j, k;
+ unsigned int tetCount = 0;
+ vtkPoints *pts = vtkPoints::New();
+ pts->Allocate(points->GetNumberOfPoints());
+ VertexManager uniqueVerts(pts);
+ ccmPolygonToTriangles tess(&uniqueVerts);
+ unsigned int oc[2] = {dom, 0};
+
+ int useCount = 0;
+ for(i = 0; i < civ.size(); ++i)
+ useCount += (civ[i].id != -1) ? 1 : 0;
+
+ originalCells[dom] = vtkUnsignedIntArray::New();
+ originalCells[dom]->SetName("avtOriginalCellNumbers");
+ originalCells[dom]->SetNumberOfComponents(2);
+ originalCells[dom]->Allocate(useCount * 3);
+
+ for (i = 0; i < civ.size(); ++i)
+ {
+ const CellInfo &ci = civ[i];
+ if(ci.id == -1)
+ continue;
+
+ oc[1] = ci.id;
+ int nFaces = ci.faces.size();
+ int nPts = 0;
+ doubleVector fbounds;
+ for (j = 0; j < nFaces; ++j)
+ {
+ nPts += ci.faces[j].nodes.size();
+ fbounds.push_back(VTK_LARGE_FLOAT);
+ fbounds.push_back(-VTK_LARGE_FLOAT);
+ fbounds.push_back(VTK_LARGE_FLOAT);
+ fbounds.push_back(-VTK_LARGE_FLOAT);
+ fbounds.push_back(VTK_LARGE_FLOAT);
+ fbounds.push_back(-VTK_LARGE_FLOAT);
+ }
+ double *pt;
+ double cbounds[6] = {VTK_LARGE_FLOAT, -VTK_LARGE_FLOAT,
+ VTK_LARGE_FLOAT, -VTK_LARGE_FLOAT,
+ VTK_LARGE_FLOAT, -VTK_LARGE_FLOAT};
+
+ int cnt = 0;
+ for (j = 0; j < nFaces; ++j)
+ {
+ const intVector &nodes = ci.faces[j].nodes;
+
+ for (k = 0; k < nodes.size(); ++k)
+ {
+ cnt++;
+ pt = points->GetPoint(nodes[k]);
+
+ if (pt[0] < cbounds[0])
+ cbounds[0] = pt[0];
+ if (pt[0] > cbounds[1])
+ cbounds[1] = pt[0];
+ if (pt[1] < cbounds[2])
+ cbounds[2] = pt[1];
+ if (pt[1] > cbounds[3])
+ cbounds[3] = pt[1];
+ if (pt[2] < cbounds[4])
+ cbounds[4] = pt[2];
+ if (pt[2] > cbounds[5])
+ cbounds[5] = pt[2];
+
+ if (pt[0] < fbounds[j*6+0])
+ fbounds[j*6+0] = pt[0];
+ if (pt[0] > fbounds[j*6+1])
+ fbounds[j*6+1] = pt[0];
+ if (pt[1] < fbounds[j*6+2])
+ fbounds[j*6+2] = pt[1];
+ if (pt[1] > fbounds[j*6+3])
+ fbounds[j*6+3] = pt[1];
+ if (pt[2] < fbounds[j*6+4])
+ fbounds[j*6+4] = pt[2];
+ if (pt[2] > fbounds[j*6+5])
+ fbounds[j*6+5] = pt[2];
+ } // k nodes
+ } // j faces
+
+ double cc[3] = {0.,0.,0.};
+ double fc[3] = {0.,0.,0.};
+ for (j = 0; j < 3; ++j)
+ cc[j] = (cbounds[2*j+1]+cbounds[2*j])/2.0;
+ int centerId = uniqueVerts.GetVertexId(cc);
+
+ for (j = 0; j < nFaces; ++j)
+ {
+ // Find the face center
+ const intVector &nodes = ci.faces[j].nodes;
+ double fc[3] = {0.,0.,0.};
+ for (k = 0; k < 3; ++k)
+ fc[k] = (fbounds[2*k+1+(6*j)]+fbounds[2*k+(6*j)])/2.0;
+
+ // Tesselate the face
+ double n[3] = {(cc[0] - fc[0]), (cc[1] - fc[1]), (cc[2] - fc[2])};
+ tess.SetNormal(n);
+ tess.BeginPolygon();
+ tess.BeginContour();
+ for (k = 0; k < nodes.size(); ++k)
+ {
+ cnt++;
+ pt = points->GetPoint(nodes[k]);
+ tess.AddVertex(pt);
+ } // k nodes
+ tess.EndContour();
+ tess.EndPolygon();
+
+ // Make a tet for each triangle in the face to the cell center.
+ vtkIdType verts[4];
+ verts[3] = centerId;
+ if (tess.GetNumTriangles() > 0)
+ {
+ for (k = 0; k < tess.GetNumTriangles(); ++k)
+ {
+ int a, b, c;
+ tess.GetTriangle(k, a, b, c);
+ verts[0] = a;
+ verts[1] = b;
+ verts[2] = c;
+ ugrid->InsertNextCell(VTK_TETRA, 4, verts);
+ ((vtkUnsignedIntArray*)originalCells[dom])->
+ InsertNextTupleValue(oc);
+ }
+ tetCount += tess.GetNumTriangles();
+ }
+ // prepare for next cell
+ tess.ClearTriangles();
+ } // end face
+ }
+ pts->Squeeze();
+
+ ugrid->SetPoints(pts);
+ pts->Delete();
+ ugrid->GetCellData()->AddArray(originalCells[dom]);
+ ugrid->GetCellData()->CopyFieldOn("avtOriginalCellNumbers");
+
+ debug4 << mName << "Input number of polyhedral cells: " << civ.size()
+ << endl;
+ debug4 << mName << "Output tetrahedral cells: " << tetCount << endl;
+#endif
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::TesselateCells2D
+//
+// Purpose:
+// Adds the 2D cells to the unstructured grid, tessellating them as needed.
+//
+// Arguments:
+// dom : The domain number
+// civ : The cell information vector.
+// points : The points that we'll allocate during tessellation.
+// ugrid : The return unstructured grid object.
+//
+// Returns:
+//
+// Note:
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Dec 18 12:54:56 PST 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Feb 28 15:00:24 PST 2008
+// Use cellid stored in CellInfo in avtOriginalCellNumbers array.
+//
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Remove unused variables.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+typedef std::pair<int,int> edge_pair;
+
+void
+avtCCMFileFormat::TesselateCells2D(const int domain, const CellInfoVector &civ,
+ vtkPoints *points, vtkUnstructuredGrid *ugrid)
+{
+#ifndef MDSERVER
+ int dom = subdividingSingleMesh ? 0 : domain;
+
+ unsigned int i, k;
+ vtkPoints *pts = vtkPoints::New();
+ pts->Allocate(points->GetNumberOfPoints());
+ VertexManager uniqueVerts(pts);
+ ccmPolygonToTriangles tess(&uniqueVerts);
+ unsigned int oc[2] = {dom, 0};
+
+ int useCount = 0;
+ for(i = 0; i < civ.size(); ++i)
+ useCount += (civ[i].id != -1) ? 1 : 0;
+
+ originalCells[dom] = vtkUnsignedIntArray::New();
+ originalCells[dom]->SetName("avtOriginalCellNumbers");
+ originalCells[dom]->SetNumberOfComponents(2);
+ originalCells[dom]->Allocate(useCount*3);
+
+ const double n[3] = {0., 0., 1.};
+
+ for (i = 0; i < civ.size(); ++i)
+ {
+ const CellInfo &ci = civ[i];
+ if(ci.id == -1)
+ continue;
+
+ oc[1] = ci.id;
+ tess.SetNormal(n);
+ tess.BeginPolygon();
+
+ // We need to sort the edges around the polygon so we can tessellate
+ // by adding the first vertex of each edge to define the contour. We
+ // could do a different approach making n-1 triangles for the number
+ // of edges in the shape but this way supports concave polygons whereas
+ // the other approach does not.
+
+ // Put all of the line segments in a pool of free edges.
+ std::set<edge_pair> freeEdges;
+ for (int f = 0; f < ci.faces.size(); ++f)
+ {
+ edge_pair e01(ci.faces[f].nodes[0], ci.faces[f].nodes[1]);
+ freeEdges.insert(e01);
+ }
+
+ while(freeEdges.size() > 0)
+ {
+ std::deque<int> shape;
+
+ // Get seed edge and remove it from the pool
+ edge_pair currentEdge;
+ if(freeEdges.begin() != freeEdges.end())
+ {
+ currentEdge = *freeEdges.begin();
+ freeEdges.erase(freeEdges.begin());
+ }
+ shape.push_back(currentEdge.first);
+ shape.push_back(currentEdge.second);
+
+ // Now, look for edges that contain either of the points in
+ // the current edge.
+ bool found;
+ do
+ {
+ found = false;
+ for(std::set<edge_pair>::iterator pos = freeEdges.begin();
+ pos != freeEdges.end() && !found; ++pos)
+ {
+ if(currentEdge.first == pos->first)
+ {
+ currentEdge.first = pos->second;
+ shape.push_front(pos->second);
+ freeEdges.erase(pos);
+ found = true;
+ }
+ else if(currentEdge.first == pos->second)
+ {
+ currentEdge.first = pos->first;
+ shape.push_front(pos->first);
+ freeEdges.erase(pos);
+ found = true;
+ }
+ else if(currentEdge.second == pos->first)
+ {
+ currentEdge.second = pos->second;
+ shape.push_back(pos->second);
+ freeEdges.erase(pos);
+ found = true;
+ }
+ else if(currentEdge.second == pos->second)
+ {
+ currentEdge.second = pos->first;
+ shape.push_back(pos->first);
+ freeEdges.erase(pos);
+ found = true;
+ }
+ }
+ } while(found);
+
+ if(shape.size() > 2)
+ {
+ tess.BeginContour();
+ for(int v = 0; v < shape.size(); ++v)
+ tess.AddVertex(points->GetPoint(shape[v]));
+ tess.EndContour();
+ }
+ }
+ tess.EndPolygon();
+
+ vtkIdType verts[4];
+ if (tess.GetNumTriangles() > 0)
+ {
+ for (k = 0; k < tess.GetNumTriangles(); ++k)
+ {
+ int a, b, c;
+ tess.GetTriangle(k, a, b, c);
+ verts[0] = a;
+ verts[1] = b;
+ verts[2] = c;
+ ugrid->InsertNextCell(VTK_TRIANGLE, 3, verts);
+
+ ((vtkUnsignedIntArray*)originalCells[dom])->
+ InsertNextTupleValue(oc);
+ }
+ }
+
+ // prepare for next cell
+ tess.ClearTriangles();
+ }
+
+ pts->Squeeze();
+ ugrid->SetPoints(pts);
+ pts->Delete();
+ ugrid->GetCellData()->AddArray(originalCells[dom]);
+ ugrid->GetCellData()->CopyFieldOn("avtOriginalCellNumbers");
+#endif
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::BuildHex
+//
+// Purpose:
+// Creates a VTK hex cell from CCM faces.
+//
+// Programmer: Kathleen Bonnell
+// Creation: October 1, 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Remove unused variables.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::BuildHex(const CellInfo &ci, vtkCellArray *cellArray,
+ intVector &cellTypes)
+{
+ unsigned int i, j;
+ const FaceInfoVector &faces = ci.faces;
+ intVector uniqueNodes;
+ bool useface;
+ bool usedBF = false;
+ int nnodes = 0;
+ vtkIdList *cellNodes = vtkIdList::New();
+
+ for (i = 0; i < faces.size(); ++i)
+ {
+ useface = true;
+ int nnodes = faces[i].nodes.size();
+ if (nnodes != 4)
+ {
+ return;
+ }
+ for (j = 0; useface && j <nnodes; ++j)
+ {
+ useface = (count(uniqueNodes.begin(), uniqueNodes.end(),
+ faces[i].nodes[j]) == 0);
+ }
+ if (useface)
+ {
+ if ((ci.faceTypes[i] == 0 && !usedBF) || (ci.faceTypes[i] == 2))
+ {
+ usedBF = (ci.faceTypes[i] == 0);
+ // reorder this face
+ for (j = 0; j < nnodes; ++j)
+ {
+ uniqueNodes.push_back(faces[i].nodes[j]);
+ }
+ cellNodes->InsertNextId(faces[i].nodes[0]);
+ cellNodes->InsertNextId(faces[i].nodes[3]);
+ cellNodes->InsertNextId(faces[i].nodes[2]);
+ cellNodes->InsertNextId(faces[i].nodes[1]);
+ }
+ else
+ {
+ for (j = 0; j < nnodes; ++j)
+ {
+ uniqueNodes.push_back(faces[i].nodes[j]);
+ cellNodes->InsertNextId(faces[i].nodes[j]);
+ }
+ }
+ }
+ }
+ cellArray->InsertNextCell(cellNodes);
+ cellTypes.push_back(VTK_HEXAHEDRON);
+ cellNodes->Delete();
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::BuildTet
+//
+// Purpose:
+// Creates a VTK tet cell from CCM faces.
+//
+// Programmer: Kathleen Bonnell
+// Creation: October 1, 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Fix '==' '=' mix-up.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::BuildTet(const CellInfo &ci, vtkCellArray *cellArray,
+ intVector &cellTypes)
+{
+ unsigned int i, j;
+ const FaceInfoVector &faces = ci.faces;
+ intVector uniqueNodes;
+ bool lastNodeFound = false;
+ vtkIdList *cellNodes = vtkIdList::New();
+
+ for (i = 0; i < faces.size(); ++i)
+ {
+ if (faces[i].nodes.size() != 3)
+ {
+ return;
+ }
+ }
+
+ if (ci.faceTypes[0] != 1)
+ {
+ uniqueNodes.push_back(faces[0].nodes[0]);
+ uniqueNodes.push_back(faces[0].nodes[2]);
+ uniqueNodes.push_back(faces[0].nodes[1]);
+ cellNodes->InsertNextId(faces[0].nodes[0]);
+ cellNodes->InsertNextId(faces[0].nodes[2]);
+ cellNodes->InsertNextId(faces[0].nodes[1]);
+ }
+ else
+ {
+ uniqueNodes.push_back(faces[0].nodes[0]);
+ uniqueNodes.push_back(faces[0].nodes[1]);
+ uniqueNodes.push_back(faces[0].nodes[2]);
+ cellNodes->InsertNextId(faces[0].nodes[0]);
+ cellNodes->InsertNextId(faces[0].nodes[1]);
+ cellNodes->InsertNextId(faces[0].nodes[2]);
+ }
+
+ for (i = 1; !lastNodeFound && i < faces.size(); ++i)
+ {
+ for (j = 0; !lastNodeFound && j <faces[i].nodes.size(); ++j)
+ {
+ if ((count(uniqueNodes.begin(), uniqueNodes.end(),
+ faces[i].nodes[j]) == 0))
+ {
+ lastNodeFound = true;
+ uniqueNodes.push_back(faces[i].nodes[j]);
+ cellNodes->InsertNextId(faces[i].nodes[j]);
+ }
+ }
+
+ }
+ cellArray->InsertNextCell(cellNodes);
+ cellTypes.push_back(VTK_TETRA);
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::BuildPyramid
+//
+// Purpose:
+// Creates a VTK pyramid cell from CCM faces.
+//
+// Programmer: Kathleen Bonnell
+// Creation: October 1, 2007
+//
+// Modifications:
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::BuildPyramid(const CellInfo &ci, vtkCellArray *cellArray,
+ intVector &cellTypes)
+{
+ unsigned int i, j;
+ const FaceInfoVector &faces = ci.faces;
+ intVector uniqueNodes;
+ int baseFace = 0;
+ vtkIdList *cellNodes = vtkIdList::New();
+
+ for (i = 0; i < faces.size(); ++i)
+ {
+ if (faces[i].nodes.size() == 4)
+ {
+ baseFace = i;
+ for (j = 0; j < 4; ++j)
+ {
+ uniqueNodes.push_back(faces[i].nodes[j]);
+ cellNodes->InsertNextId(faces[i].nodes[j]);
+ }
+ }
+ }
+ int triFace = (baseFace+1) % 4;
+ for (i = 0; i < 3; ++i) // only 3 nodes in the triangle-face
+ {
+ if ((count(uniqueNodes.begin(), uniqueNodes.end(),
+ faces[triFace].nodes[i])) == 0)
+ {
+ cellNodes->InsertNextId(faces[triFace].nodes[i]);
+ break;
+ }
+ }
+
+ cellArray->InsertNextCell(cellNodes);
+ cellTypes.push_back(VTK_PYRAMID);
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::BuildWedge
+//
+// Purpose:
+// Creates a VTK wedge cell from CCM faces.
+//
+// Programmer: Kathleen Bonnell
+// Creation: October 1, 2007
+//
+// Modifications:
+// Kathleen Bonnell, Thu Mar 6 09:21:02 PST 2008
+// Remove unused variables.
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Fixed a bug in which cell and vertex ids were mistaken for 1-based
+// indices.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::BuildWedge(const CellInfo &ci, vtkCellArray *cellArray,
+ intVector &cellTypes)
+{
+ unsigned int i;
+ const FaceInfoVector &faces = ci.faces;
+ vtkIdList *cellNodes = vtkIdList::New();
+ for (i = 0; i < faces.size(); ++i)
+ {
+ if (faces[i].nodes.size() == 3)
+ {
+ cellNodes->InsertNextId(faces[i].nodes[0]);
+ cellNodes->InsertNextId(faces[i].nodes[1]);
+ cellNodes->InsertNextId(faces[i].nodes[2]);
+ }
+ }
+ cellArray->InsertNextCell(cellNodes);
+ cellTypes.push_back(VTK_WEDGE);
+}
+
+
+
+
+
+
+//
+// avtCCMFileFormat::FaceInfo class
+//
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::FaceInfo::FaceInfo
+//
+// Purpose:
+// Constructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Removed unused face id.
+//
+// ****************************************************************************
+
+avtCCMFileFormat::FaceInfo::FaceInfo()
+{
+ //id = -1;
+ cells[0] = -1;
+ cells[1] = -1;
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::FaceInfo::FaceInfo
+//
+// Purpose:
+// Copy Constructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Removed unused face id.
+//
+// ****************************************************************************
+
+avtCCMFileFormat::FaceInfo::FaceInfo(const avtCCMFileFormat::FaceInfo &obj)
+{
+ //id = obj.id;
+ cells[0] = obj.cells[0];
+ cells[1] = obj.cells[1];
+ nodes = obj.nodes;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::FaceInfo::~FaceInfo
+//
+// Purpose:
+// Destructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+avtCCMFileFormat::FaceInfo::~FaceInfo()
+{
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::FaceInfo::operator =
+//
+// Purpose:
+// Assignment operator
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// Dave Bremer, Fri Apr 4 16:29:49 PDT 2008
+// Removed unused face id.
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::FaceInfo::operator =(const avtCCMFileFormat::FaceInfo &obj)
+{
+ //id = obj.id;
+ cells[0] = obj.cells[0];
+ cells[1] = obj.cells[1];
+ nodes = obj.nodes;
+}
+
+//
+// avtCCMFileFormat::CellInfo class
+//
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::CellInfo
+//
+// Purpose:
+// Constructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+avtCCMFileFormat::CellInfo::CellInfo()
+{
+ id = -1;
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::CellInfo
+//
+// Purpose:
+// Copy Constructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+avtCCMFileFormat::CellInfo::CellInfo(const avtCCMFileFormat::CellInfo &obj)
+{
+ id = obj.id;
+ faceTypes = obj.faceTypes;
+ faces = obj.faces;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::~CellInfo
+//
+// Purpose:
+// Destructor
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+avtCCMFileFormat::CellInfo::~CellInfo()
+{
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::operator =
+//
+// Purpose:
+// Assignment operator
+//
+// Programmer: Kathleen Bonnell
+// Creation: September 6, 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::CellInfo::operator =(const avtCCMFileFormat::CellInfo &obj)
+{
+ id = obj.id;
+ faceTypes = obj.faceTypes;
+ faces = obj.faces;
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::CellCenter
+//
+// Purpose:
+// Determines the cell center from its nodes.
+//
+// Arguments:
+// center : The returned center.
+// pts : The global points array that stores all of the nodes.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Oct 1 14:02:37 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::CellInfo::CellCenter(double *center, vtkPoints *pts) const
+{
+ int npts = 0;
+ double c[3] = {0.,0.,0.};
+ for(int i = 0; i < faces.size(); ++i)
+ {
+ for(int j = 0; j < faces[i].nodes.size(); ++j, ++npts)
+ {
+ double *pt = pts->GetPoint(faces[i].nodes[j]);
+ c[0] += pt[0];
+ c[1] += pt[1];
+ c[2] += pt[2];
+ }
+ }
+ center[0] = c[0] / double(npts);
+ center[1] = c[1] / double(npts);
+ center[2] = c[2] / double(npts);
+}
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::CellInfo::UseNodes
+//
+// Purpose:
+// Sets a true value into a domain-node-sized bool array for each node that
+// the cell uses.
+//
+// Arguments:
+// pts : The array containing the values for whether a node is used.
+//
+// Programmer: Brad Whitlock
+// Creation: Thu Oct 1 14:03:22 PDT 2009
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::CellInfo::UseNodes(bool *pts) const
+{
+ for(int i = 0; i < faces.size(); ++i)
+ for(int j = 0; j < faces[i].nodes.size(); ++j)
+ pts[faces[i].nodes[j]] = true;
+}
+
+// ****************************************************************************
+// Method: operator << (ostream &os, CCMIOEntity e)
+//
+// Purpose:
+// Prints a CCMIOEntity
+//
+// Arguments:
+// os : The stream to which we'll print.
+// e : The value to print.
+//
+// Returns: The input stream.
+//
+// Note:
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Dec 18 15:16:45 PST 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+ostream &
+operator << (ostream &os, CCMIOEntity e)
+{
+ switch(e)
+ {
+ case kCCMIOBoundaryRegion:
+ os << "kCCMIOBoundaryRegion"; break;
+ case kCCMIOFieldData:
+ os << "kCCMIOFieldData"; break;
+ case kCCMIOFieldPhase:
+ os << "kCCMIOFieldPhase"; break;
+ case kCCMIOInternalFaces:
+ os << "kCCMIOInternalFaces"; break;
+ case kCCMIOMaxEntity:
+ os << "kCCMIOMaxEntity"; break;
+ case kCCMIOProblemDescription:
+ os << "kCCMIOProblemDescription"; break;
+ case kCCMIOReferenceData:
+ os << "kCCMIOReferenceData"; break;
+ case kCCMIOMap:
+ os << "kCCMIOMap"; break;
+ case kCCMIONull:
+ os << "kCCMIONull"; break;
+ case kCCMIOTopology:
+ os << "kCCMIOTopology"; break;
+ case kCCMIOVertices:
+ os << "kCCMIOVertices"; break;
+ case kCCMIOBoundaryFaces :
+ os << "kCCMIOBoundaryFaces "; break;
+ case kCCMIOCellType:
+ os << "kCCMIOCellType"; break;
+ case kCCMIOCells:
+ os << "kCCMIOCells"; break;
+ case kCCMIOField:
+ os << "kCCMIOField"; break;
+ case kCCMIOFieldSet :
+ os << "kCCMIOFieldSet "; break;
+ case kCCMIOInterfaces:
+ os << "kCCMIOInterfaces"; break;
+ case kCCMIOLagrangianData :
+ os << "kCCMIOLagrangianData "; break;
+ case kCCMIOModelConstants :
+ os << "kCCMIOModelConstants "; break;
+ case kCCMIOProcessor :
+ os << "kCCMIOProcessor "; break;
+ case kCCMIOProstarSet:
+ os << "kCCMIOProstarSet"; break;
+ case kCCMIORestart :
+ os << "kCCMIORestart "; break;
+ case kCCMIORestartData:
+ os << "kCCMIORestartData"; break;
+ case kCCMIOState :
+ os << "kCCMIOState "; break;
+ }
+
+ return os;
+}
+
+// ****************************************************************************
+// Method: operator << (ostream &os, const CCMIONode &node)
+//
+// Purpose:
+// Prints a CCMIONode
+//
+// Arguments:
+// os : The stream to which we'll print.
+// e : The value to print.
+//
+// Returns: The input stream.
+//
+// Note:
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Dec 18 15:16:45 PST 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+ostream &
+operator << (ostream &os, const CCMIONode &node)
+{
+ _CCMIONode *n = (_CCMIONode *)&node;
+ os << "(node=" << n->node << ", parent=" << n->parent << ")";
+ return os;
+}
+
+// ****************************************************************************
+// Method: operator << (ostream &os, const CCMIOID &e)
+//
+// Purpose:
+// Prints a CCMIOID
+//
+// Arguments:
+// os : The stream to which we'll print.
+// e : The value to print.
+//
+// Returns: The input stream.
+//
+// Note:
+//
+// Programmer: Brad Whitlock
+// Creation: Tue Dec 18 15:16:45 PST 2007
+//
+// Modifications:
+//
+// ****************************************************************************
+
+ostream &
+operator << (ostream &os, const CCMIOID &id)
+{
+ os << "(root=" << id.root
+ << ", node=" << id.node
+ << ", id=" << id.id
+ << ", type=" << id.type
+ << ", version=" << id.version << ")";
+ return os;
+}
+
+
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::IDMap::IDMap
+//
+// Purpose:
+// Constructor
+//
+// Programmer: Dave Bremer
+// Creation: Fri Apr 4 16:29:49 PDT 2008
+//
+// Modifications:
+//
+// ****************************************************************************
+
+avtCCMFileFormat::IDMap::IDMap()
+{
+ bSequential = false;
+ bReverseMap = false;
+ iFirstElem = 0;
+ numIDs = 0;
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::IDMap::SetIDs
+//
+// Purpose:
+// Analyze a set of IDs and build a structure for mapping them back
+// to the index corresponding to their position in v.
+//
+// Programmer: Dave Bremer
+// Creation: Fri Apr 4 16:29:49 PDT 2008
+//
+// Modifications:
+//
+// ****************************************************************************
+
+void
+avtCCMFileFormat::IDMap::SetIDs(const intVector &v)
+{
+ numIDs = v.size();
+
+ bSequential = true;
+ bReverseMap = false;
+ iFirstElem = v[0];
+
+ int min = v[0], max = v[0];
+
+ int ii;
+ for (ii = 1; ii < v.size(); ii++)
+ {
+ if (v[ii-1]+1 != v[ii])
+ bSequential = false;
+ if (v[ii] < min)
+ min = v[ii];
+ if (v[ii] > max)
+ max = v[ii];
+ }
+
+ //Don't bother copying data in, in this case.
+ if (bSequential)
+ return;
+
+ if (max-min+1 <= v.size()*2)
+ {
+ bReverseMap = true;
+ iFirstElem = min;
+
+ ids.resize(max-min+1, -1);
+ for (ii = 0; ii < v.size(); ii++)
+ {
+ ids[v[ii]-iFirstElem] = ii;
+ }
+ }
+ else
+ {
+ bReverseMap = false;
+
+ ids.resize(v.size()*2);
+
+ for (ii = 0; ii < v.size(); ii++)
+ {
+ ids[ii*2] = v[ii];
+ ids[ii*2+1] = ii;
+ }
+ qsort( &ids[0], v.size(), sizeof(int)*2, avtCCMFileFormat::IDMap::compare);
+ }
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::IDMap::compare
+//
+// Purpose:
+// compare function for qsort
+//
+// Programmer: Dave Bremer
+// Creation: Fri Apr 4 16:29:49 PDT 2008
+//
+// Modifications:
+//
+// ****************************************************************************
+
+int
+avtCCMFileFormat::IDMap::compare(const void *a, const void *b)
+{
+ return ( *((int *)a) - *((int *)b) );
+}
+
+
+// ****************************************************************************
+// Method: avtCCMFileFormat::IDMap::IDtoIndex
+//
+// Purpose:
+// Map an id to an array index. Return -1 if the id is not in the map.
+//
+// Programmer: Dave Bremer
+// Creation: Fri Apr 4 16:29:49 PDT 2008
+//
+// Modifications:
+//
+// ****************************************************************************
+
+int
+avtCCMFileFormat::IDMap::IDtoIndex(int id) const
+{
+ if (bSequential)
+ {
+ int r = id - iFirstElem;
+ if (r >= 0 && r < numIDs)
+ return r;
+ }
+ else if (bReverseMap)
+ {
+ if (id-iFirstElem >= 0 && id-iFirstElem < ids.size())
+ return ids[id-iFirstElem];
+ }
+ else
+ {
+ int min = 0, max = ids.size()/2 - 1;
+ int mid = (min+max)/2;
+
+ while (min <= max)
+ {
+ if (id < ids[mid*2])
+ {
+ max = mid-1;
+ mid = (min+max)/2;
+ }
+ else if (id > ids[mid*2])
+ {
+ min = mid+1;
+ mid = (min+max)/2;
+ }
+ else
+ {
+ return ids[mid*2+1];
+ }
+ }
+ }
+ return -1;
+}
+
+
+*/
Added: MOAB/trunk/ReadCCMIO.hpp
===================================================================
--- MOAB/trunk/ReadCCMIO.hpp (rev 0)
+++ MOAB/trunk/ReadCCMIO.hpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -0,0 +1,123 @@
+#ifndef READCCMIO_HPP
+#define READCCMIO_HPP
+
+#ifndef IS_BUILDING_MB
+#error "ReadCCMIO.hpp isn't supposed to be included into an application"
+#endif
+
+#include <vector>
+#include <set>
+#include <map>
+#include <string>
+
+#include "MBForward.hpp"
+#include "MBRange.hpp"
+#include "ExoIIInterface.hpp"
+#include "MBReaderIface.hpp"
+#include "TupleList.hpp"
+#include "ccmio.h"
+
+#define TupleList std::map<int, std::vector<MBEntityHandle> >
+#define SenseList std::map<int, std::vector<int> >
+#undef TUPLE_LIST
+
+class MBReadUtilIface;
+
+class MB_DLL_EXPORT ReadCCMIO : public MBReaderIface
+{
+
+public:
+
+ //! Constructor
+ ReadCCMIO(MBInterface *impl);
+
+ //! Destructor
+ virtual ~ReadCCMIO();
+
+ static MBReaderIface* factory( MBInterface* );
+
+ MBErrorCode load_file(const char *file_name,
+ const MBEntityHandle* file_set,
+ const FileOptions& opts,
+ const MBReaderIface::IDTag* subset_list,
+ int subset_list_length,
+ const MBTag* file_id_tag);
+
+private:
+
+ MBErrorCode read_processor(CCMIOID rootID, CCMIOID processorID, CCMIOSize_t proc);
+
+
+ MBErrorCode read_cells(CCMIOSize_t proc, CCMIOID processorID,
+ CCMIOID verticesID, CCMIOID topologyID,
+ CCMIOID solutionID, bool has_solution,
+ TupleList &vert_map);
+
+
+ MBErrorCode construct_cells(TupleList &face_map,
+#ifndef TUPLE_LIST
+ SenseList &sense_map,
+#endif
+ TupleList &vert_map, MBRange &new_cells);
+
+
+ MBErrorCode create_cell_from_faces(std::vector<MBEntityHandle> &facehs,
+ std::vector<int> &senses,
+ MBEntityHandle &cell);
+
+
+ MBErrorCode read_all_faces(CCMIOID topologyID, TupleList &vert_map,
+ TupleList &face_map
+#ifndef TUPLE_LIST
+ ,SenseList &sense_map
+#endif
+ );
+
+
+ MBErrorCode read_faces(CCMIOID faceID, CCMIOEntity bdy_or_int,
+ TupleList &vert_map,
+ TupleList &face_map
+#ifndef TUPLE_LIST
+ ,SenseList &sense_map
+#endif
+ );
+
+ MBErrorCode make_faces(int *farray,
+ TupleList &vert_map,
+ std::vector<MBEntityHandle> &new_faces);
+
+ MBErrorCode read_vertices(CCMIOSize_t proc, CCMIOID processorID, CCMIOID verticesID,
+ CCMIOID topologyID, CCMIOID solutionID, bool has_solution,
+ MBRange &verts, TupleList &vert_map);
+
+
+ MBErrorCode get_processors(CCMIOID stateID, CCMIOID &processorID,
+ std::set<CCMIOSize_t> &procs);
+
+
+ MBErrorCode get_state(CCMIOID rootID, CCMIOID &problemID, CCMIOID &stateID);
+
+
+ virtual MBErrorCode read_tag_values( const char* file_name,
+ const char* tag_name,
+ const FileOptions& opts,
+ std::vector<int>& tag_values_out,
+ const IDTag* subset_list = 0,
+ int subset_list_length = 0 );
+
+ //! Cached tags for reading. Note that all these tags are defined when the
+ //! core is initialized.
+ MBTag mMaterialSetTag;
+ MBTag mDirichletSetTag;
+ MBTag mNeumannSetTag;
+ MBTag mHasMidNodesTag;
+ MBTag mGlobalIdTag;
+
+ MBInterface *mbImpl;
+
+ MBReadUtilIface* readMeshIface;
+
+ bool hasSolution;
+};
+
+#endif
Modified: MOAB/trunk/ReadSmf.cpp
===================================================================
--- MOAB/trunk/ReadSmf.cpp 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/ReadSmf.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -132,7 +132,7 @@
mPartitionTagName = partition_tag_name;
std::ifstream smfFile;
- smfFile.open( filename , std::ifstream::in);
+ //smfFile.open( filename , std::ifstream::in);
if (!smfFile.is_open())
{
return MB_FILE_DOES_NOT_EXIST;
Modified: MOAB/trunk/SMF_State.cpp
===================================================================
--- MOAB/trunk/SMF_State.cpp 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/SMF_State.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -1,8 +1,10 @@
#include "SMF_State.hpp"
+#include <cstring>
+#include <cstdlib>
-inline int streq(const char *a,const char *b) { return std::strcmp(a,b)==0; }
+inline int streq(const char *a,const char *b) { return strcmp(a,b)==0; }
SMF_State::SMF_State(const SMF_ivars& ivar, SMF_State *link)
Added: MOAB/trunk/WriteCCMIO.cpp
===================================================================
--- MOAB/trunk/WriteCCMIO.cpp (rev 0)
+++ MOAB/trunk/WriteCCMIO.cpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -0,0 +1,1035 @@
+/**
+ * MOAB, a Mesh-Oriented datABase, is a software component for creating,
+ * storing and accessing finite element mesh data.
+ *
+ * Copyright 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government
+ * retains certain rights in this software.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ */
+
+
+#ifdef WIN32
+#ifdef _DEBUG
+// turn off warnings that say they debugging identifier has been truncated
+// this warning comes up when using some STL containers
+#pragma warning(disable : 4786)
+#endif
+#endif
+
+
+#include "WriteCCMIO.hpp"
+#include "ccmio.h"
+#include "ccmioutility.h"
+#include "ccmiocore.h"
+#include <utility>
+#include <algorithm>
+#include <time.h>
+#include <string>
+#include <vector>
+#include <stdio.h>
+#include <iostream>
+#include <algorithm>
+
+#include "MBInterface.hpp"
+#include "MBRange.hpp"
+#include "MBCN.hpp"
+#include "assert.h"
+#include "MBInternals.hpp"
+#include "ExoIIUtil.hpp"
+#include "MBTagConventions.hpp"
+#include "MBWriteUtilIface.hpp"
+
+static char const kStateName[] = "default";
+
+static const int ccm_types[] = {
+ 1, // MBVERTEX
+ 2, // MBEDGE
+ -1, // MBTRI
+ -1, // MBQUAD
+ -1, // MBPOLYGON
+ 13, // MBTET
+ 14, // MBPYRAMID
+ 12, // MBPRISM
+ -1, // MBKNIFE
+ 11, // MBHEX
+ 255 // MBPOLYHEDRON
+};
+
+#define INS_ID(stringvar, prefix, id) \
+ sprintf(stringvar, prefix, id)
+
+MBWriterIface* WriteCCMIO::factory( MBInterface* iface )
+{ return new WriteCCMIO( iface ); }
+
+WriteCCMIO::WriteCCMIO(MBInterface *impl)
+ : mbImpl(impl), mCurrentMeshHandle(0)
+{
+ assert(impl != NULL);
+
+ void* ptr = 0;
+ impl->query_interface( "MBWriteUtilIface", &ptr );
+ mWriteIface = reinterpret_cast<MBWriteUtilIface*>(ptr);
+
+ // initialize in case tag_get_handle fails below
+ //! get and cache predefined tag handles
+ int dum_val = 0;
+ MBErrorCode result = impl->tag_get_handle(MATERIAL_SET_TAG_NAME, mMaterialSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(MATERIAL_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, mMaterialSetTag,
+ &dum_val);
+
+ result = impl->tag_get_handle(DIRICHLET_SET_TAG_NAME, mDirichletSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(DIRICHLET_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, mDirichletSetTag,
+ &dum_val);
+
+ result = impl->tag_get_handle(NEUMANN_SET_TAG_NAME, mNeumannSetTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(NEUMANN_SET_TAG_NAME, sizeof(int), MB_TAG_SPARSE, mNeumannSetTag,
+ &dum_val);
+
+#ifdef USE_MPI
+ result = impl->tag_get_handle(PARALLEL_PARTITION_TAG_NAME, mPartitionSetTag);
+ // no need to check result, if it's not there, we don't create one
+#endif
+
+ result = impl->tag_get_handle(HAS_MID_NODES_TAG_NAME, mHasMidNodesTag);
+ if (MB_TAG_NOT_FOUND == result) {
+ int dum_val_array[] = {0, 0, 0, 0};
+ result = impl->tag_create(HAS_MID_NODES_TAG_NAME, 4*sizeof(int), MB_TAG_SPARSE, mHasMidNodesTag,
+ dum_val_array);
+ }
+
+ result = impl->tag_get_handle(GLOBAL_ID_TAG_NAME, mGlobalIdTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create(GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_SPARSE, mGlobalIdTag,
+ &dum_val);
+
+ dum_val = -1;
+ result = impl->tag_get_handle("__matSetIdTag", mMatSetIdTag);
+ if (MB_TAG_NOT_FOUND == result)
+ result = impl->tag_create("__matSetIdTag", sizeof(int), MB_TAG_DENSE, mMatSetIdTag,
+ &dum_val);
+
+
+ impl->tag_create("WriteCCMIO element mark", 1, MB_TAG_BIT, mEntityMark, NULL);
+
+}
+
+WriteCCMIO::~WriteCCMIO()
+{
+ std::string iface_name = "MBWriteUtilIface";
+ mbImpl->release_interface(iface_name, mWriteIface);
+
+ mbImpl->tag_delete(mEntityMark);
+
+}
+
+void WriteCCMIO::reset_matset(std::vector<WriteCCMIO::MaterialSetData> &matset_info)
+{
+ std::vector<WriteCCMIO::MaterialSetData>::iterator iter;
+
+ for (iter = matset_info.begin(); iter != matset_info.end(); iter++)
+ {
+ delete (*iter).elements;
+ }
+}
+
+MBErrorCode WriteCCMIO::write_file(const char *file_name,
+ const bool /* overwrite (commented out to remove warning) */,
+ const FileOptions& opts,
+ const MBEntityHandle *ent_handles,
+ const int num_sets,
+ const std::vector<std::string>&,
+ const MBTag* ,
+ int ,
+ int )
+{
+ assert(0 != mMaterialSetTag &&
+ 0 != mNeumannSetTag &&
+ 0 != mDirichletSetTag);
+
+ // check the file name
+ if (NULL == strstr(file_name, ".ccmio"))
+ return MB_FAILURE;
+
+ std::vector<MBEntityHandle> matsets, dirsets, neusets, partsets, entities;
+
+ fileName = file_name;
+
+ // separate into material sets, dirichlet sets, neumann sets
+
+ if (num_sets == 0) {
+ // default to all defined sets
+ MBRange this_range;
+ mbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range);
+ std::copy(this_range.begin(), this_range.end(), std::back_inserter(matsets));
+ this_range.clear();
+ mbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range);
+ std::copy(this_range.begin(), this_range.end(), std::back_inserter(dirsets));
+ this_range.clear();
+ mbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range);
+ std::copy(this_range.begin(), this_range.end(), std::back_inserter(neusets));
+ if (mPartitionSetTag) {
+ this_range.clear();
+ mbImpl->get_entities_by_type_and_tag(0, MBENTITYSET, &mPartitionSetTag, NULL, 1, this_range);
+ std::copy(this_range.begin(), this_range.end(), std::back_inserter(partsets));
+ }
+ }
+
+ else {
+ int dummy;
+ for (const MBEntityHandle *iter = ent_handles; iter < ent_handles+num_sets; iter++)
+ {
+ if (MB_SUCCESS == mbImpl->tag_get_data(mMaterialSetTag, &(*iter), 1, &dummy))
+ matsets.push_back(*iter);
+ else if (MB_SUCCESS == mbImpl->tag_get_data(mDirichletSetTag, &(*iter), 1, &dummy))
+ dirsets.push_back(*iter);
+ else if (MB_SUCCESS == mbImpl->tag_get_data(mNeumannSetTag, &(*iter), 1, &dummy))
+ neusets.push_back(*iter);
+ else if (mPartitionSetTag &&
+ MB_SUCCESS == mbImpl->tag_get_data(mPartitionSetTag, &(*iter), 1, &dummy))
+ partsets.push_back(*iter);
+ }
+ }
+
+ MBErrorCode result = mbImpl->tag_get_handle(HAS_MID_NODES_TAG_NAME, mHasMidNodesTag);
+ if (MB_TAG_NOT_FOUND == result) {
+ int dum_val_array[] = {0, 0, 0, 0};
+ result = mbImpl->tag_create(HAS_MID_NODES_TAG_NAME, 4*sizeof(int), MB_TAG_SPARSE, mHasMidNodesTag,
+ dum_val_array);
+ }
+
+ // if there is nothing to write just return.
+ if (matsets.empty() && dirsets.empty() && neusets.empty() && partsets.empty())
+ return MB_FILE_WRITE_ERROR;
+
+ std::vector<WriteCCMIO::MaterialSetData> matset_info;
+ std::vector<WriteCCMIO::DirichletSetData> dirset_info;
+ std::vector<WriteCCMIO::NeumannSetData> neuset_info;
+
+ MeshInfo mesh_info;
+
+ matset_info.clear();
+ if(gather_mesh_information(mesh_info, matset_info, neuset_info, dirset_info,
+ matsets, neusets, dirsets) != MB_SUCCESS)
+ {
+ reset_matset(matset_info);
+ return MB_FAILURE;
+ }
+
+
+ // try to open the file after gather mesh info succeeds
+ CCMIOSize_t i = CCMIOSIZEC(0);
+ CCMIOID stateID, processorID, rootID;
+ CCMIOError error = kCCMIONoErr;
+
+ CCMIOOpenFile(&error, file_name, kCCMIOWrite, &rootID);
+ if(kCCMIONoErr != error)
+ {
+ mWriteIface->report_error("Cannot open %s", file_name);
+ return MB_FAILURE;
+ }
+
+ // Create a new state (or re-use an existing one).
+ CCMIONewState(&error, rootID, kStateName, NULL, NULL, &stateID);
+ if (kCCMIONoErr != error) {
+ reset_matset(matset_info);
+ return MB_FAILURE;
+ }
+
+// for (; i < CCMIOSIZEC(partsets.size()); i++) {
+ if (CCMIONextEntity(NULL, stateID, kCCMIOProcessor, &i, &processorID) != kCCMIONoErr)
+ CCMIONewEntity(&error, stateID, kCCMIOProcessor, NULL, &processorID);
+
+ // Get rid of any data that may be in this processor (if the state was
+ // not new).
+ else
+ CCMIOClearProcessor(&error, stateID, processorID, TRUE, TRUE, TRUE, TRUE,
+ TRUE);
+// }
+
+ int *vgids;
+ if( write_nodes(rootID, mesh_info.nodes, mesh_info.num_dim, vgids)
+ != MB_SUCCESS )
+ {
+ reset_matset(matset_info);
+ return MB_FAILURE;
+ }
+
+ if( write_matsets(mesh_info, matset_info, neuset_info, mesh_info.nodes, vgids) )
+ {
+ reset_matset(matset_info);
+ return MB_FAILURE;
+ }
+
+ if (write_problem_description(rootID, stateID)) {
+ return MB_FAILURE;
+ }
+
+ CCMIOCloseFile(&error, rootID);
+
+ if (error != kCCMIONoErr)
+ {
+ return MB_FAILURE;
+ }
+
+ // The CCMIO library uses ADF to store the actual data. Unfortunately,
+ // ADF leaks disk space; deleting a node does not recover all the disk
+ // space. Now that everything is successfully written it might be useful
+ // to call CCMIOCompress() here to ensure that the file is as small as
+ // possible. Please see the Core API documentation for caveats on its
+ // usage.
+ if (CCMIOCompress(NULL, const_cast<char*>(file_name)) != kCCMIONoErr)
+ {
+ std::cout << "Error compressing file. Check that you have "
+ << "adequate disk space " << std::endl << "and that you have write "
+ << "permission to the current directory." << std::endl;
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::write_problem_description(CCMIOID rootID, CCMIOID stateID)
+{
+ // Write out a dummy problem description. If we happen to know that
+ // there already is a problem description previously recorded that
+ // is valid we could skip this step.
+ CCMIOID problem, constants, id;
+ CCMIOError error = kCCMIONoErr;
+
+ CCMIONewEntity(&error, rootID, kCCMIOProblemDescription, "Dummy description",
+ &problem);
+ CCMIONewIndexedEntity(&error, problem, kCCMIOCellType, 1, "Dummy celltypes", &id);
+ CCMIOWriteOptstr(&error, id, "MaterialType", "solid");
+ CCMIONewIndexedEntity(&error, problem, kCCMIOCellType, 2, "Dummy celltypes", &id);
+ CCMIOWriteOptstr(&error, id, "MaterialType", "solid");
+
+ CCMIONewEntity(&error, problem, kCCMIOModelConstants, "Constant values",
+ &constants);
+ CCMIOWriteOptf(&error, constants, "Gravity", 9.82);
+ CCMIOWriteOptf(&error, constants, "B.P. of water", 373);
+
+ // We have problem description recorded but our state does not know
+ // about it. So tell the state that it has a problem description.
+ CCMIOWriteState(&error, stateID, problem, "Example state");
+
+ return MB_SUCCESS;
+}
+
+
+MBErrorCode WriteCCMIO::gather_mesh_information(MeshInfo &mesh_info,
+ std::vector<WriteCCMIO::MaterialSetData> &matset_info,
+ std::vector<WriteCCMIO::NeumannSetData> &neuset_info,
+ std::vector<WriteCCMIO::DirichletSetData> &dirset_info,
+ std::vector<MBEntityHandle> &matsets,
+ std::vector<MBEntityHandle> &neusets,
+ std::vector<MBEntityHandle> &dirsets)
+{
+
+ std::vector<MBEntityHandle>::iterator vector_iter, end_vector_iter;
+
+ mesh_info.num_nodes = 0;
+ mesh_info.num_elements = 0;
+ mesh_info.num_matsets = 0;
+
+ int id = 0;
+
+ vector_iter= matsets.begin();
+ end_vector_iter = matsets.end();
+
+ mesh_info.num_matsets = matsets.size();
+
+ std::vector<MBEntityHandle> parent_meshsets;
+
+ // clean out the bits for the element mark
+ mbImpl->tag_delete(mEntityMark);
+ mbImpl->tag_create("WriteCCMIO element mark", 1, MB_TAG_BIT, mEntityMark, NULL);
+
+ int highest_dimension_of_element_matsets = 0;
+
+ for(vector_iter = matsets.begin(); vector_iter != matsets.end(); vector_iter++)
+ {
+
+ WriteCCMIO::MaterialSetData matset_data;
+ matset_data.elements = new MBRange;
+
+ //for the purpose of qa records, get the parents of these matsets
+ if( mbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS )
+ return MB_FAILURE;
+
+ // get all Entity Handles in the mesh set
+ MBRange dummy_range;
+ mbImpl->get_entities_by_handle(*vector_iter, dummy_range, true );
+
+ // find the dimension of the last entity in this range
+ MBRange::iterator entity_iter = dummy_range.end();
+ entity_iter = dummy_range.end();
+ entity_iter--;
+ int this_dim = MBCN::Dimension(TYPE_FROM_HANDLE(*entity_iter));
+ entity_iter = dummy_range.begin();
+ while (entity_iter != dummy_range.end() &&
+ MBCN::Dimension(TYPE_FROM_HANDLE(*entity_iter)) != this_dim)
+ entity_iter++;
+
+ if (entity_iter != dummy_range.end())
+ std::copy(entity_iter, dummy_range.end(), mb_range_inserter(*(matset_data.elements)));
+
+ assert(matset_data.elements->begin() == matset_data.elements->end() ||
+ MBCN::Dimension(TYPE_FROM_HANDLE(*(matset_data.elements->begin()))) == this_dim);
+
+ // get the matset's id
+ if(mbImpl->tag_get_data(mMaterialSetTag, &(*vector_iter), 1, &id) != MB_SUCCESS ) {
+ mWriteIface->report_error("Couldn't get matset id from a tag for an element matset.");
+ return MB_FAILURE;
+ }
+
+ matset_data.id = id;
+ matset_data.number_attributes = 0;
+
+ // iterate through all the elements in the meshset
+ MBRange::iterator elem_range_iter, end_elem_range_iter;
+ elem_range_iter = matset_data.elements->begin();
+ end_elem_range_iter = matset_data.elements->end();
+
+ // get the entity type for this matset, verifying that it's the same for all elements
+ // THIS ASSUMES HANDLES SORT BY TYPE!!!
+ MBEntityType entity_type = TYPE_FROM_HANDLE(*elem_range_iter);
+ end_elem_range_iter--;
+ if (entity_type != TYPE_FROM_HANDLE(*(end_elem_range_iter++))) {
+ mWriteIface->report_error("Entities in matset %i not of common type", id);
+ return MB_FAILURE;
+ }
+
+ int dimension = MBCN::Dimension(entity_type);
+
+ if( dimension > highest_dimension_of_element_matsets )
+ highest_dimension_of_element_matsets = dimension;
+
+ matset_data.moab_type = mbImpl->type_from_handle(*(matset_data.elements->begin()));
+ if (MBMAXTYPE == matset_data.moab_type) return MB_FAILURE;
+
+ std::vector<MBEntityHandle> tmp_conn;
+ mbImpl->get_connectivity(&(*(matset_data.elements->begin())), 1, tmp_conn);
+ matset_data.element_type =
+ ExoIIUtil::get_element_type_from_num_verts(tmp_conn.size(), entity_type, dimension);
+
+ if (matset_data.element_type == EXOII_MAX_ELEM_TYPE) {
+ mWriteIface->report_error("Element type in matset %i didn't get set correctly", id);
+ return MB_FAILURE;
+ }
+
+ matset_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[matset_data.element_type];
+
+ // number of nodes for this matset
+ matset_data.number_elements = matset_data.elements->size();
+
+ // total number of elements
+ mesh_info.num_elements += matset_data.number_elements;
+
+ // get the nodes for the elements
+ mWriteIface->gather_nodes_from_elements(*matset_data.elements, mEntityMark, mesh_info.nodes);
+
+ if(!neusets.empty())
+ {
+ // if there are neusets, keep track of which elements are being written out
+ for(MBRange::iterator iter = matset_data.elements->begin();
+ iter != matset_data.elements->end(); ++iter)
+ {
+ unsigned char bit = 0x1;
+ mbImpl->tag_set_data(mEntityMark, &(*iter), 1, &bit);
+ }
+ }
+
+ matset_info.push_back( matset_data );
+
+ }
+
+
+ //if user hasn't entered dimension, we figure it out
+ if( mesh_info.num_dim == 0 )
+ {
+ //never want 1 or zero dimensions
+ if( highest_dimension_of_element_matsets < 2 )
+ mesh_info.num_dim = 3;
+ else
+ mesh_info.num_dim = highest_dimension_of_element_matsets;
+ }
+
+ MBRange::iterator range_iter, end_range_iter;
+ range_iter = mesh_info.nodes.begin();
+ end_range_iter = mesh_info.nodes.end();
+
+ mesh_info.num_nodes = mesh_info.nodes.size();
+
+ //------dirsets--------
+
+ vector_iter= dirsets.begin();
+ end_vector_iter = dirsets.end();
+
+ for(; vector_iter != end_vector_iter; vector_iter++)
+ {
+
+ WriteCCMIO::DirichletSetData dirset_data;
+ dirset_data.id = 0;
+ dirset_data.number_nodes = 0;
+
+ // get the dirset's id
+ if(mbImpl->tag_get_data(mDirichletSetTag,&(*vector_iter), 1,&id) != MB_SUCCESS) {
+ mWriteIface->report_error("Couldn't get id tag for dirset %i", id);
+ return MB_FAILURE;
+ }
+
+ dirset_data.id = id;
+
+ std::vector<MBEntityHandle> node_vector;
+ //get the nodes of the dirset that are in mesh_info.nodes
+ if( mbImpl->get_entities_by_handle(*vector_iter, node_vector, true) != MB_SUCCESS ) {
+ mWriteIface->report_error("Couldn't get nodes in dirset %i", id);
+ return MB_FAILURE;
+ }
+
+ std::vector<MBEntityHandle>::iterator iter, end_iter;
+ iter = node_vector.begin();
+ end_iter= node_vector.end();
+
+ int j=0;
+ unsigned char node_marked = 0;
+ MBErrorCode result;
+ for(; iter != end_iter; iter++)
+ {
+ if (TYPE_FROM_HANDLE(*iter) != MBVERTEX) continue;
+ result = mbImpl->tag_get_data(mEntityMark, &(*iter), 1, &node_marked);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get mark data.");
+ return result;
+ }
+
+ if(node_marked == 0x1) dirset_data.nodes.push_back( *iter );
+ j++;
+ }
+
+ dirset_data.number_nodes = dirset_data.nodes.size();
+ dirset_info.push_back( dirset_data );
+ }
+
+ //------neusets--------
+ vector_iter= neusets.begin();
+ end_vector_iter = neusets.end();
+
+ for(; vector_iter != end_vector_iter; vector_iter++)
+ {
+ WriteCCMIO::NeumannSetData neuset_data;
+
+ // get the neuset's id
+ if(mbImpl->tag_get_data(mNeumannSetTag,&(*vector_iter), 1,&id) != MB_SUCCESS)
+ return MB_FAILURE;
+
+ neuset_data.id = id;
+ neuset_data.mesh_set_handle = *vector_iter;
+
+ //get the sides in two lists, one forward the other reverse; starts with forward sense
+ // by convention
+ MBRange forward_elems, reverse_elems;
+ if(get_neuset_elems(*vector_iter, 0, forward_elems, reverse_elems) == MB_FAILURE)
+ return MB_FAILURE;
+
+ MBErrorCode result = get_valid_sides(forward_elems, 1, neuset_data);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get valid sides data.");
+ return result;
+ }
+ result = get_valid_sides(reverse_elems, -1, neuset_data);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get valid sides data.");
+ return result;
+ }
+
+ neuset_data.number_elements = neuset_data.elements.size();
+ neuset_info.push_back( neuset_data );
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::get_valid_sides(MBRange &elems, const int sense,
+ WriteCCMIO::NeumannSetData &neuset_data)
+{
+ // this is where we see if underlying element of side set element is included in output
+
+ unsigned char element_marked = 0;
+ MBErrorCode result;
+ for(MBRange::iterator iter = elems.begin(); iter != elems.end(); iter++)
+ {
+ // should insert here if "side" is a quad/tri on a quad/tri mesh
+ result = mbImpl->tag_get_data(mEntityMark, &(*iter), 1, &element_marked);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get mark data.");
+ return result;
+ }
+
+ if(element_marked == 0x1)
+ {
+ neuset_data.elements.push_back( *iter );
+
+ // TJT TODO: the sense should really be # edges + 1or2
+ neuset_data.side_numbers.push_back((sense == 1 ? 1 : 2));
+ }
+ else //then "side" is probably a quad/tri on a hex/tet mesh
+ {
+ std::vector<MBEntityHandle> parents;
+ int dimension = MBCN::Dimension( TYPE_FROM_HANDLE(*iter));
+
+ //get the adjacent parent element of "side"
+ if( mbImpl->get_adjacencies( &(*iter), 1, dimension+1, false, parents) != MB_SUCCESS ) {
+ mWriteIface->report_error("Couldn't get adjacencies for neuset.");
+ return MB_FAILURE;
+ }
+
+ if(!parents.empty())
+ {
+ //make sure the adjacent parent element will be output
+ for(unsigned int k=0; k<parents.size(); k++)
+ {
+ result = mbImpl->tag_get_data(mEntityMark, &(parents[k]), 1, &element_marked);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get mark data.");
+ return result;
+ }
+
+ int side_no, this_sense, this_offset;
+ if(element_marked == 0x1 &&
+ mbImpl->side_number(parents[k], *iter, side_no,
+ this_sense, this_offset) == MB_SUCCESS &&
+ this_sense == sense) {
+ neuset_data.elements.push_back(parents[k]);
+ neuset_data.side_numbers.push_back(side_no+1);
+ break;
+ }
+ }
+ }
+ else
+ {
+ mWriteIface->report_error("No parent element exists for element in neuset %i", neuset_data.id);
+ return MB_FAILURE;
+ }
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::get_gids(const MBRange &ents, int *&gids,
+ int &minid, int &maxid)
+{
+ int num_ents = ents.size();
+ gids = new int[num_ents];
+ MBErrorCode result = mbImpl->tag_get_data(mGlobalIdTag, ents, &gids[0]);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get global id data.");
+ return result;
+ }
+ minid = *std::min_element(gids, gids+num_ents);
+ maxid = *std::max_element(gids, gids+num_ents);
+ if (0 == minid) {
+ // gids need to be assigned
+ for (int i = 1; i <= num_ents; i++) gids[i] = i;
+ result = mbImpl->tag_set_data(mGlobalIdTag, ents, &gids[0]);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't set global id data.");
+ return result;
+ }
+ maxid = num_ents;
+ }
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::write_nodes(CCMIOID rootID, const MBRange& nodes,
+ const int dimension, int *&vgids)
+{
+ // get/write map (global ids) first
+ const int num_nodes = nodes.size();
+ int minid, maxid;
+ MBErrorCode result = get_gids(nodes, vgids, minid, maxid);
+ if (MB_SUCCESS != result) return result;
+
+ CCMIOID mapID;
+ CCMIOError error;
+ CCMIONewEntity(&error, rootID, kCCMIOMap, "Vertex map", &mapID);
+ CCMIOWriteMap(&error, mapID, CCMIOSIZEC(nodes.size()),
+ CCMIOSIZEC(maxid), vgids,
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+
+ // get the vertex locations
+ double *coords = new double[dimension*num_nodes];
+ std::vector<double*> coord_arrays(3);
+ coord_arrays[0] = coords;
+ coord_arrays[1] = coords+num_nodes;
+ coord_arrays[2] = (dimension == 3 ? coords+2*num_nodes : NULL);
+ result = mWriteIface->get_node_arrays(dimension, num_nodes, nodes,
+ mGlobalIdTag, 0, coord_arrays);
+ if(result != MB_SUCCESS)
+ {
+ delete [] coords;
+ return result;
+ }
+
+ // transform coordinates, if necessary
+ result = transform_coords(dimension, num_nodes, coords);
+ if(result != MB_SUCCESS)
+ {
+ delete [] coords;
+ return result;
+ }
+
+
+ // write the vertices
+ CCMIOID verticesID;
+ CCMIONewEntity(&error, rootID, kCCMIOVertices, "Vertices", &verticesID);
+ CCMIOWriteVerticesd(&error, verticesID,
+ CCMIOSIZEC(dimension*num_nodes), 1.0, mapID, coords,
+ CCMIOINDEXC(1), CCMIOINDEXC(dimension));
+ if (error != kCCMIONoErr) {
+ mWriteIface->report_error("CCMIOWriteVertices failed.");
+ return result;
+ }
+
+ // clean up
+ delete [] coords;
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::transform_coords(const int dimension, const int num_nodes, double *coords)
+{
+ MBTag trans_tag;
+ MBErrorCode result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, trans_tag);
+ if( result == MB_TAG_NOT_FOUND ) return MB_SUCCESS;
+ double trans_matrix[16];
+ result = mbImpl->tag_get_data( trans_tag, NULL, 0, trans_matrix );
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get transform data.");
+ return result;
+ }
+
+ double *tmp_coords = coords;
+ for( int i=0; i<num_nodes; i++, tmp_coords += 1) {
+ double vec1[3] = {0.0, 0.0, 0.0};
+ for( int row=0; row<3; row++ ) {
+ vec1[row] += ( trans_matrix[ (row*4)+0 ] * coords[0]);
+ vec1[row] += ( trans_matrix[ (row*4)+1 ] * coords[num_nodes]);
+ if (3 == dimension) vec1[row] += ( trans_matrix[ (row*4)+2 ] * coords[2*num_nodes]);
+ }
+
+ coords[0] = vec1[0];
+ coords[num_nodes] = vec1[1];
+ coords[2*num_nodes] = vec1[2];
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::write_matsets(MeshInfo & /* mesh_info (commented out to remove warning) */,
+ std::vector<WriteCCMIO::MaterialSetData> &matset_data,
+ std::vector<WriteCCMIO::NeumannSetData> &/* neuset_data */,
+ // (commented out to remove warning)
+ MBRange &verts,
+ const int *vgids)
+{
+ std::vector<int> connect;
+ MBErrorCode result;
+ CCMIOID rootID, cellMapID, topologyID, id;
+
+ // don't usually have anywhere near 31 nodes per element
+ connect.reserve(31);
+ MBRange::const_iterator rit;
+
+ MBRange all_elems;
+ for (unsigned int i = 0; i < matset_data.size(); i++)
+ all_elems.merge(*(matset_data[i].elements));
+
+
+ const int num_elems = all_elems.size();
+ int *egids, minid, maxid;
+ result = get_gids(all_elems, egids, minid, maxid);
+ if (MB_SUCCESS != result) return result;
+
+ // Write the cells
+ CCMIOError error;
+ CCMIONewEntity(&error, rootID, kCCMIOMap, "Cell map", &cellMapID);
+ CCMIOWriteMap(&error, cellMapID, CCMIOSIZEC(all_elems.size()),
+ CCMIOSIZEC(maxid), egids,
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ CCMIONewEntity(&error, rootID, kCCMIOTopology, "Topology", &topologyID);
+ CCMIONewEntity(&error, topologyID, kCCMIOCells, "Cells", &id);
+
+ // get cell types
+ int *ctypes = new int[all_elems.size()];
+ int i = 0;
+ rit = all_elems.begin();
+ for (; i < num_elems; i++, rit++) {
+ ctypes[i] = ccm_types[mbImpl->type_from_handle(*rit)];
+ assert(-1 != ctypes[i]);
+ }
+
+ CCMIOWriteCells(&error, id, cellMapID, ctypes,
+ CCMIOINDEXC(1), CCMIOINDEXC(num_elems));
+ delete [] ctypes;
+
+ // Write the faces
+ // first, allocate a tag of length 6 (max # faces per region, except
+ // for polyhedra)
+ MBTag mark_tag;
+ short int def_val = 0;
+ result = mbImpl->tag_create("__mark", 1, MB_TAG_DENSE, MB_TYPE_OPAQUE,
+ mark_tag, &def_val);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't create mark tag.");
+ return result;
+ }
+
+ // now faces
+ unsigned char markt;
+ std::vector<MBEntityHandle> tmp_face_cells, storage;
+ std::vector<int> iface_connect, iface_cells;
+ std::vector<int> eface_connect, eface_cells;
+ MBEntityHandle tmp_connect[MBCN::MAX_NODES_PER_ELEMENT]; // tmp connect vector
+ const MBEntityHandle *connectc; int num_connectc; // cell connectivity
+ const MBEntityHandle *connectf; int num_connectf; // face connectivity
+ i = 0;
+ rit = all_elems.begin();
+ for (; i < num_elems; i++, rit++) {
+ // if not polyh, get mark
+ MBEntityType etype = TYPE_FROM_HANDLE(*rit);
+ if (MBPOLYHEDRON != etype && MBPOLYGON != etype) {
+ result = mbImpl->tag_get_data(mark_tag, &(*rit), 1, &markt);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get mark data.");
+ return result;
+ }
+ }
+
+ // get connect
+ result = mbImpl->get_connectivity(*rit, connectc, num_connectc, false, &storage);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get entity connectivity.");
+ return result;
+ }
+
+ // if polyh, write faces directly
+ bool is_polyh = (MBPOLYHEDRON == etype);
+
+ int num_faces = MBCN::NumSubEntities(etype, 2);
+ // for each face (from MBCN)
+ for (int f = 0; f < num_faces; f++) {
+ // if this face marked, skip
+ if (!is_polyh && ((markt >> f) & 0x1)) continue;
+
+ // get face connect
+ if (!is_polyh) {
+ // (from MBCN)
+ MBCN::SubEntityConn(connectc, etype, 2, f, tmp_connect, num_connectf);
+ connectf = tmp_connect;
+ }
+ else {
+ // get face connectivity
+ result = mbImpl->get_connectivity(connectc[f], connectf, num_connectf, false);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get polyhedron connectivity.");
+ return result;
+ }
+ }
+
+ // get adj cells from face connect
+ tmp_face_cells.clear();
+ result = mbImpl->get_adjacencies(connectf, num_connectf, 3, false, tmp_face_cells);
+ if (MB_SUCCESS != result || tmp_face_cells.empty()) {
+ mWriteIface->report_error("Error getting adj hexes.");
+ return result;
+ }
+
+ bool is_internal = (tmp_face_cells.size() == 2);
+ if (!is_polyh && is_internal) {
+ // make sure 1st is forward sense
+ int side_num, sense, offset;
+ MBCN::SideNumber(etype, connectc, connectf, num_connectf,
+ 2, side_num, sense, offset);
+ if (sense == 1 && tmp_face_cells[0] != *rit) {
+ assert(2 == tmp_face_cells.size());
+ MBEntityHandle tmph = tmp_face_cells[0];
+ tmp_face_cells[1] = tmp_face_cells[0];
+ tmp_face_cells[0] = tmph;
+ }
+ }
+
+ // get ids of cells in all_elems
+ std::vector<int> *fcells_ptr, *fconnect_ptr;
+ fcells_ptr = (is_internal ? &iface_cells : &eface_cells);
+ fconnect_ptr = (is_internal ? &iface_connect : &eface_connect);
+ fcells_ptr->push_back(egids[all_elems.index(tmp_face_cells[0])]);
+ if (is_internal) fcells_ptr->push_back(egids[all_elems.index(tmp_face_cells[1])]);
+ fconnect_ptr->push_back(num_connectf);
+
+ // get indices of face vertices, add one
+ for (int fv = 0; fv < num_connectf; fv++)
+ fconnect_ptr->push_back(vgids[verts.index(connectf[fv])]);
+
+ if (!is_polyh && is_internal) {
+ // mark other cell for this face, if there is another cell
+ MBEntityHandle other_cell = tmp_face_cells[0];
+ const MBEntityHandle *connecto; int num_connecto;
+ if (other_cell == *rit) other_cell = tmp_face_cells[1];
+ result = mbImpl->get_connectivity(other_cell, connecto, num_connecto,
+ false, &storage);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get other entity connectivity.");
+ return result;
+ }
+ // get side number
+ int side_num, sense, offset;
+ MBCN::SideNumber(TYPE_FROM_HANDLE(other_cell), connecto, connectf, num_connectf,
+ 2, side_num, sense, offset);
+ // set mark for this face
+ short int tmp_mark, tmp_mark2;
+ result = mbImpl->tag_get_data(mark_tag, &other_cell, 1, &tmp_mark);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't get mark data for other cell.");
+ return result;
+ }
+ tmp_mark2 = 0x1 << (unsigned int)side_num;
+ assert("mark for this side on other entity shouldn't be set already" &&
+ !(tmp_mark & tmp_mark2));
+ tmp_mark |= tmp_mark2;
+ result = mbImpl->tag_set_data(mark_tag, &other_cell, 1, &tmp_mark);
+ if (MB_SUCCESS != result) {
+ mWriteIface->report_error("Couldn't set mark data for other cell.");
+ return result;
+ }
+ } // !is_polyh
+ } // loop over faces in elem
+ } // loop over elems
+
+ int num_ifaces = iface_cells.size()/2,
+ num_efaces = eface_cells.size();
+
+ // write internal faces
+ CCMIOID mapID;
+ CCMIONewEntity(&error, rootID, kCCMIOMap, NULL, &mapID);
+ // set gids for internal faces
+ if ((int)all_elems.size() < (num_ifaces + num_efaces)) {
+ delete [] egids;
+ egids = new int[num_ifaces + num_efaces];
+ }
+ for (int i = 1; i <= (num_ifaces+num_efaces); i++) egids[i-1] = i;
+ CCMIOWriteMap(&error, mapID, CCMIOSIZEC(num_ifaces),
+ CCMIOSIZEC(num_ifaces+num_efaces),
+ &egids[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ CCMIONewEntity(&error, topologyID, kCCMIOInternalFaces, "Internal faces", &id);
+ CCMIOWriteFaces(&error, id, kCCMIOInternalFaces, mapID,
+ CCMIOSIZEC(iface_connect.size()), &iface_connect[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ CCMIOWriteFaceCells(&error, id, kCCMIOInternalFaces, mapID, &iface_cells[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+
+ // now external faces
+ CCMIONewEntity(&error, rootID, kCCMIOMap, NULL, &mapID);
+ CCMIOWriteMap(&error, mapID, CCMIOSIZEC(num_efaces),
+ CCMIOSIZEC(num_efaces), &egids[num_ifaces],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ CCMIONewIndexedEntity(&error, topologyID, kCCMIOBoundaryFaces, 0,
+ "Boundary faces", &id);
+ CCMIOWriteFaces(&error, id, kCCMIOBoundaryFaces, mapID,
+ CCMIOSIZEC(eface_connect.size()), &eface_connect[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+ CCMIOWriteFaceCells(&error, id, kCCMIOBoundaryFaces, mapID, &eface_cells[0],
+ CCMIOINDEXC(kCCMIOStart), CCMIOINDEXC(kCCMIOEnd));
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::open_file(const char* filename)
+{
+ // not a valid filname
+ if(strlen((const char*)filename) == 0)
+ {
+ mWriteIface->report_error("Output filename not specified");
+ return MB_FAILURE;
+ }
+
+ /* template - open file & store somewhere */
+
+ // file couldn't be opened
+ return MB_SUCCESS;
+}
+
+MBErrorCode WriteCCMIO::get_neuset_elems(MBEntityHandle neuset, int current_sense,
+ MBRange &forward_elems, MBRange &reverse_elems)
+{
+ MBRange neuset_elems, neuset_meshsets;
+
+ // get the sense tag; don't need to check return, might be an error if the tag
+ // hasn't been created yet
+ MBTag sense_tag = 0;
+ mbImpl->tag_get_handle("SENSE", sense_tag);
+
+ // get the entities in this set
+ MBErrorCode result = mbImpl->get_entities_by_handle(neuset, neuset_elems, true);
+ if (MB_FAILURE == result) return result;
+
+ // now remove the meshsets into the neuset_meshsets; first find the first meshset,
+ MBRange::iterator range_iter = neuset_elems.begin();
+ while (TYPE_FROM_HANDLE(*range_iter) != MBENTITYSET && range_iter != neuset_elems.end())
+ range_iter++;
+
+ // then, if there are some, copy them into neuset_meshsets and erase from neuset_elems
+ if (range_iter != neuset_elems.end()) {
+ std::copy(range_iter, neuset_elems.end(), mb_range_inserter(neuset_meshsets));
+ neuset_elems.erase(range_iter, neuset_elems.end());
+ }
+
+
+ // ok, for the elements, check the sense of this set and copy into the right range
+ // (if the sense is 0, copy into both ranges)
+
+ // need to step forward on list until we reach the right dimension
+ MBRange::iterator dum_it = neuset_elems.end();
+ dum_it--;
+ int target_dim = MBCN::Dimension(TYPE_FROM_HANDLE(*dum_it));
+ dum_it = neuset_elems.begin();
+ while (target_dim != MBCN::Dimension(TYPE_FROM_HANDLE(*dum_it)) &&
+ dum_it != neuset_elems.end())
+ dum_it++;
+
+ if (current_sense == 1 || current_sense == 0)
+ std::copy(dum_it, neuset_elems.end(), mb_range_inserter(forward_elems));
+ if (current_sense == -1 || current_sense == 0)
+ std::copy(dum_it, neuset_elems.end(), mb_range_inserter(reverse_elems));
+
+ // now loop over the contained meshsets, getting the sense of those and calling this
+ // function recursively
+ for (range_iter = neuset_meshsets.begin(); range_iter != neuset_meshsets.end(); range_iter++) {
+
+ // first get the sense; if it's not there, by convention it's forward
+ int this_sense;
+ if (0 == sense_tag ||
+ MB_FAILURE == mbImpl->tag_get_data(sense_tag, &(*range_iter), 1, &this_sense))
+ this_sense = 1;
+
+ // now get all the entities on this meshset, with the proper (possibly reversed) sense
+ get_neuset_elems(*range_iter, this_sense*current_sense,
+ forward_elems, reverse_elems);
+ }
+
+ return result;
+}
+
+
+
Added: MOAB/trunk/WriteCCMIO.hpp
===================================================================
--- MOAB/trunk/WriteCCMIO.hpp (rev 0)
+++ MOAB/trunk/WriteCCMIO.hpp 2010-01-29 05:24:00 UTC (rev 3508)
@@ -0,0 +1,196 @@
+/**
+ * MOAB, a Mesh-Oriented datABase, is a software component for creating,
+ * storing and accessing finite element mesh data.
+ *
+ * Copyright 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government
+ * retains certain rights in this software.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ */
+
+//-------------------------------------------------------------------------
+// Filename : WriteCCMIO.hpp
+//
+// Purpose : ExodusII writer
+//
+// Special Notes : Lots of code taken from verde implementation
+//
+// Creator : Corey Ernst
+//
+// Date : 8/02
+//
+// Owner : Corey Ernst
+//-------------------------------------------------------------------------
+
+#ifndef WRITECCMIO_HPP
+#define WRITECCMIO_HPP
+
+#ifndef IS_BUILDING_MB
+#error "WriteCCMIO.hpp isn't supposed to be included into an application"
+#endif
+
+#include <vector>
+#include <string>
+
+#include "MBForward.hpp"
+#include "MBRange.hpp"
+#include "ExoIIInterface.hpp"
+#include "MBWriterIface.hpp"
+#include "ccmio.h"
+
+class MBWriteUtilIface;
+
+class MB_DLL_EXPORT WriteCCMIO : public MBWriterIface
+{
+
+public:
+
+ //! Constructor
+ WriteCCMIO(MBInterface *impl);
+
+ //! Destructor
+ virtual ~WriteCCMIO();
+
+ static MBWriterIface* factory( MBInterface* );
+
+ //! writes out a file
+ MBErrorCode write_file(const char *file_name,
+ const bool overwrite,
+ const FileOptions& opts,
+ const MBEntityHandle *output_list,
+ const int num_sets,
+ const std::vector<std::string>& qa_list,
+ const MBTag* tag_list,
+ int num_tags,
+ int export_dimension);
+
+//! struct used to hold data for each block to be output; used by
+//! initialize_file to initialize the file header for increased speed
+ struct MaterialSetData
+ {
+ int id;
+ int number_elements;
+ int number_nodes_per_element;
+ int number_attributes;
+ ExoIIElementType element_type;
+ MBEntityType moab_type;
+ MBRange *elements;
+ };
+
+//! struct used to hold data for each nodeset to be output; used by
+//! initialize_file to initialize the file header for increased speed
+ struct DirichletSetData
+ {
+ int id;
+ int number_nodes;
+ std::vector< MBEntityHandle > nodes;
+ std::vector< double > node_dist_factors;
+
+ };
+
+//! struct used to hold data for each sideset to be output; used by
+//! initialize_file to initialize the file header for increased speed
+ struct NeumannSetData
+ {
+ int id;
+ int number_elements;
+ std::vector<MBEntityHandle> elements;
+ std::vector<int> side_numbers;
+ MBEntityHandle mesh_set_handle;
+ };
+
+
+protected:
+
+ //! number of dimensions in this file
+ //int number_dimensions();
+
+ //! open a file for writing
+ MBErrorCode open_file(const char *filename);
+
+ //! contains the general information about a mesh
+ class MeshInfo
+ {
+ public:
+ unsigned int num_dim;
+ unsigned int num_nodes;
+ unsigned int num_elements;
+ unsigned int num_matsets;
+ unsigned int num_dirsets;
+ unsigned int num_neusets;
+ MBRange nodes;
+
+ MeshInfo()
+ : num_dim(0), num_nodes(0), num_elements(0), num_matsets(0),
+ num_dirsets(0), num_neusets(0)
+ {}
+
+ };
+
+private:
+
+ //! interface instance
+ MBInterface *mbImpl;
+ MBWriteUtilIface* mWriteIface;
+
+ //! file name
+ std::string fileName;
+
+ //! Meshset Handle for the mesh that is currently being read
+ MBEntityHandle mCurrentMeshHandle;
+
+ //! Cached tags for reading. Note that all these tags are defined when the
+ //! core is initialized.
+ MBTag mMaterialSetTag;
+ MBTag mDirichletSetTag;
+ MBTag mNeumannSetTag;
+ MBTag mPartitionSetTag;
+ MBTag mHasMidNodesTag;
+ MBTag mGlobalIdTag;
+ MBTag mMatSetIdTag;
+
+ MBTag mEntityMark; //used to say whether an entity will be exported
+
+ MBErrorCode gather_mesh_information(MeshInfo &mesh_info,
+ std::vector<MaterialSetData> &matset_info,
+ std::vector<NeumannSetData> &neuset_info,
+ std::vector<DirichletSetData> &dirset_info,
+ std::vector<MBEntityHandle> &matsets,
+ std::vector<MBEntityHandle> &neusets,
+ std::vector<MBEntityHandle> &dirsets);
+
+ MBErrorCode initialize_file(MeshInfo &mesh_info);
+
+ MBErrorCode write_nodes(CCMIOID rootID, const MBRange& nodes,
+ const int dimension, int *&vgids);
+
+ // get global ids for these entities; allocates gids and passes back,
+ // caller is responsible for deleting
+ MBErrorCode get_gids(const MBRange &ents, int *&gids,
+ int &minid, int &maxid);
+
+ MBErrorCode write_matsets(MeshInfo &mesh_info,
+ std::vector<MaterialSetData> &matset_data,
+ std::vector<NeumannSetData> &neuset_data,
+ MBRange &verts,
+ const int *vgids);
+
+ MBErrorCode get_valid_sides(MBRange &elems, const int sense,
+ WriteCCMIO::NeumannSetData &neuset_data);
+
+ void reset_matset(std::vector<MaterialSetData> &matset_info);
+
+ MBErrorCode get_neuset_elems(MBEntityHandle neuset, int current_sense,
+ MBRange &forward_elems, MBRange &reverse_elems);
+
+ MBErrorCode transform_coords(const int dimension, const int num_nodes, double *coords);
+
+ MBErrorCode write_problem_description(CCMIOID rootID, CCMIOID stateID);
+};
+
+#endif
Modified: MOAB/trunk/configure.ac
===================================================================
--- MOAB/trunk/configure.ac 2010-01-28 21:02:20 UTC (rev 3507)
+++ MOAB/trunk/configure.ac 2010-01-29 05:24:00 UTC (rev 3508)
@@ -375,6 +375,23 @@
fi
################################################################################
+# CCMIO OPTIONS
+################################################################################
+
+old_LDFLAGS="$LDFLAGS"
+LDFLAGS="$LDFLAGS $HDF5_LDFLAGS"
+FATHOM_CHECK_CCMIO
+LDFLAGS="$old_LDFLAGS"
+if test "xno" != "x$HAVE_CCMIO"; then
+ DEFINES="$DEFINES -DCCMIO_FILE"
+fi
+AM_CONDITIONAL(CCMIO_FILE, [test "xno" != "x$HAVE_CCMIO"])
+AM_CPPFLAGS="$CCMIO_CPPFLAGS $AM_CPPFLAGS"
+EXPORT_LDFLAGS="$CCMIO_LDFLAGS $EXPORT_LDFLAGS"
+AC_SUBST(CCMIO_LIBS)
+
+
+################################################################################
# NetCDF OPTIONS
################################################################################
Added: MOAB/trunk/m4/ccmio.m4
===================================================================
--- MOAB/trunk/m4/ccmio.m4 (rev 0)
+++ MOAB/trunk/m4/ccmio.m4 2010-01-29 05:24:00 UTC (rev 3508)
@@ -0,0 +1,73 @@
+#######################################################################################
+# Check for NetCDF library ((C++)
+# Sets HAVE_NETCDF to 'yes' or 'no'
+# If HAVE_NETCDF == yes, then exports:
+# NETCDF_CPPFLAGS
+# NETCDF_LDFLAGS
+# NETCDF_LIBS
+#######################################################################################
+AC_DEFUN([FATHOM_CHECK_CCMIO],[
+
+AC_MSG_CHECKING([if CCMIO support is enabled])
+AC_ARG_WITH(ccmio,
+[AC_HELP_STRING([--with-ccmio=DIR], [Specify CCMIO location])
+AC_HELP_STRING([--without-ccmio], [Disable support for CCMIO file format])],
+[CCMIO_ARG=$withval
+DISTCHECK_CONFIGURE_FLAGS="$DISTCHECK_CONFIGURE_FLAGS --with-ccmio=\"${withval}\""
+]
+, [CCMIO_ARG=])
+if test "xno" = "x$CCMIO_ARG"; then
+ AC_MSG_RESULT([no])
+else
+ AC_MSG_RESULT([yes])
+fi
+
+ # if CCMIO support is not disabled
+HAVE_CCMIO=no
+if test "xno" != "x$CCMIO_ARG"; then
+ HAVE_CCMIO=yes
+
+ # if a path is specified, update LIBS and INCLUDES accordingly
+ if test "xyes" != "x$CCMIO_ARG" && test "x" != "x$CCMIO_ARG"; then
+ if test -d "${CCMIO_ARG}/linux/64/lib"; then
+ CCMIO_LDFLAGS="-L${CCMIO_ARG}/linux/64/lib"
+ elif test -d "${CCMIO_ARG}"; then
+ CCMIO_LDFLAGS="-L${CCMIO_ARG}/lib"
+ elif test -d "${CCMIO_ARG}/lib"; then
+ CCMIO_LDFLAGS="-L${CCMIO_ARG}"
+ else
+ AC_MSG_ERROR("$CCMIO_ARG is not a directory.")
+ fi
+ if test -d "${CCMIO_ARG}/include"; then
+ CCMIO_CPPFLAGS="-I${CCMIO_ARG}/include"
+ else
+ CCMIO_CPPFLAGS="-I${CCMIO_ARG}"
+ fi
+ fi
+
+ old_CPPFLAGS="$CPPFLAGS"
+ CPPFLAGS="$CCMIO_CPPFLAGS $CPPFLAGS"
+ old_LDFLAGS="$LDFLAGS"
+ LDFLAGS="$CCMIO_LDFLAGS $LDFLAGS"
+
+ # Check for C library
+ AC_CHECK_HEADERS( [ccmio.h], [], [AC_MSG_WARN([[Ccmio header not found.]]); HAVE_CCMIO=no] )
+ AC_CHECK_HEADERS( [ccmioutility.h], [], [AC_MSG_WARN([[Ccmio header not found.]]); HAVE_CCMIO=no] )
+ AC_CHECK_HEADERS( [ccmiocore.h], [], [AC_MSG_WARN([[Ccmio header not found.]]); HAVE_CCMIO=no] )
+ CPPFLAGS="$old_CPPFLAGS"
+ LDFLAGS="$old_LDFLAGS"
+
+ if test "x$HAVE_CCMIO" = "xno"; then
+ if test "x$CCMIO_ARG" != "x"; then
+ AC_MSG_ERROR("Ccmio not found or not working")
+ else
+ AC_MSG_WARN("Ccmio support disabled")
+ fi
+ CCMIO_CPPFLAGS=
+ CCMIO_LDFLAGS=
+ else
+ CCMIO_LIBS="-lccmio -ladf"
+ fi
+fi
+
+]) # FATHOM_HAVE_CCMIO
More information about the moab-dev
mailing list