[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sat Jul 12 17:49:39 CDT 2014
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/f34c493671a4/
Changeset: f34c493671a4
Branch: None
User: judajake
Date: 2014-04-05 01:10:51
Summary: update vtk moab
I updated the vtk moab to match https://github.com/robertmaynard/Sandbox.
I did one small change to fix a bug with reading materials
Affected #: 19 files
diff --git a/tools/vtkMOABReaderNew/CMakeLists.txt.in b/tools/vtkMOABReaderNew/CMakeLists.txt.in
index 8b8f888..e005289 100644
--- a/tools/vtkMOABReaderNew/CMakeLists.txt.in
+++ b/tools/vtkMOABReaderNew/CMakeLists.txt.in
@@ -13,13 +13,22 @@ find_package(ParaView REQUIRED)
include(${PARAVIEW_USE_FILE})
include_directories(${PARAVIEW_INCLUDE_DIRS})
include_directories(@srcdir@)
+include_directories(@srcdir@/detail)
set(headers
- @srcdir@/SimpleMoab.h
- @srcdir@/CellTypeToType.h
+ @srcdir@/CellSets.h
@srcdir@/DataSetConverter.h
- @srcdir@/MixedCellConnectivity.h
- @srcdir@/vtkMoabReader.h
+ @srcdir@/detail/CellTypeToType.h
+ @srcdir@/detail/ContinousCellInfo.h
+ @srcdir@/detail/LoadGeometry.h
+ @srcdir@/detail/MixedCellConnectivity.h
+ @srcdir@/detail/ReadSparseTag.h
+ @srcdir@/detail/ReduceSpectralMesh.h
+ @srcdir@/detail/ReduceSpectralMesh.h
+ @srcdir@/detail/UsageTable.h
+ @srcdir@/ExtractShell.h
+ @srcdir@/FaceSets.h
+ @srcdir@/SimpleMoab.h
)
add_paraview_plugin(vtkMoabReaderPlugin "5.0"
diff --git a/tools/vtkMOABReaderNew/CellSets.h b/tools/vtkMOABReaderNew/CellSets.h
new file mode 100644
index 0000000..6171aec
--- /dev/null
+++ b/tools/vtkMOABReaderNew/CellSets.h
@@ -0,0 +1,75 @@
+
+#ifndef __smoab_CellSets_h
+#define __smoab_CellSets_h
+
+#include "SimpleMoab.h"
+
+namespace smoab
+{
+//----------------------------------------------------------------------------
+class CellSet
+{
+public:
+ CellSet(smoab::EntityHandle p,const smoab::Range& cells):
+ Entity(p),
+ Cells(cells)
+ {}
+
+ const smoab::Range& cells() const { return this->Cells; }
+ EntityHandle entity() const { return this->Entity; }
+
+ bool contains(smoab::EntityHandle c) const
+ {
+ return this->Cells.find(c) != this->Cells.end();
+ }
+
+ void erase(smoab::Range cells)
+ {
+ //seems that erase() has a bug, so use subtract
+ this->Cells = smoab::subtract(this->Cells,cells);
+ }
+
+private:
+ smoab::EntityHandle Entity;
+ smoab::Range Cells;
+};
+
+//----------------------------------------------------------------------------
+//CellSets are just a vector of CellSets
+typedef std::vector<CellSet> CellSets;
+
+//----------------------------------------------------------------------------
+//templated so it works with FaceCellSets and CellSets
+template<typename T>
+smoab::Range getParents(const T& set)
+{
+ typedef typename T::const_iterator iterator;
+ smoab::Range result;
+
+ for(iterator i=set.begin(); i != set.end(); ++i)
+ {
+ result.insert(i->entity());
+ }
+ return result;
+}
+
+//----------------------------------------------------------------------------
+//templated so it works with FaceCellSets and CellSets
+template<typename T>
+smoab::Range getAllCells(const T& set)
+{
+ typedef typename T::const_iterator iterator;
+ smoab::Range result;
+
+ for(iterator i=set.begin(); i != set.end(); ++i)
+ {
+ smoab::Range c = i->cells();
+ result.insert(c.begin(),c.end());
+ }
+ return result;
+}
+
+
+}
+
+#endif
diff --git a/tools/vtkMOABReaderNew/CellTypeToType.h b/tools/vtkMOABReaderNew/CellTypeToType.h
deleted file mode 100644
index 6ff91af..0000000
--- a/tools/vtkMOABReaderNew/CellTypeToType.h
+++ /dev/null
@@ -1,59 +0,0 @@
-#ifndef __smoab_CellTypeToType_h
-#define __smoab_CellTypeToType_h
-
-#include "SimpleMoab.h"
-#include "vtkCellType.h"
-
-namespace smoab
-{
-template<int> struct CellTypeToType;
-
-int vtkCellType(moab::EntityType t, int &num_connect)
- {
- int ctype = -1;
- switch (t)
- {
- case moab::MBEDGE:
- if (num_connect == 2) ctype = VTK_LINE;
- else if (num_connect == 3) ctype = VTK_QUADRATIC_EDGE;
- break;
- case moab::MBTRI:
- if (num_connect == 3) ctype = VTK_TRIANGLE;
- else if (num_connect == 6) ctype = VTK_QUADRATIC_TRIANGLE;
- else if (num_connect == 7) ctype = VTK_BIQUADRATIC_TRIANGLE;
- break;
- case moab::MBQUAD:
- if (num_connect == 4) ctype = VTK_QUAD;
- else if (num_connect == 8) ctype = VTK_QUADRATIC_QUAD;
- else if (num_connect == 9) ctype = VTK_BIQUADRATIC_QUAD;
- break;
- case moab::MBPOLYGON:
- if (num_connect == 4) ctype = VTK_POLYGON;
- break;
- case moab::MBTET:
- if (num_connect == 4) ctype = VTK_TETRA;
- else if (num_connect == 10) ctype = VTK_QUADRATIC_TETRA;
- break;
- case moab::MBPYRAMID:
- if (num_connect == 5) ctype = VTK_PYRAMID;
- else if (num_connect == 13) ctype = VTK_QUADRATIC_PYRAMID;
- break;
- case moab::MBPRISM:
- if (num_connect == 6) ctype = VTK_WEDGE;
- else if (num_connect == 15) ctype = VTK_QUADRATIC_WEDGE;
- break;
- case moab::MBHEX:
- if (num_connect == 8) ctype = VTK_HEXAHEDRON;
- else if (num_connect == 20) ctype = VTK_QUADRATIC_HEXAHEDRON;
- else if (num_connect == 21) ctype = VTK_QUADRATIC_HEXAHEDRON, num_connect = 20;
- else if (num_connect == 27) ctype = VTK_TRIQUADRATIC_HEXAHEDRON;
- break;
- default:
- ctype = -1;
- break;
- }
- return ctype;
- }
-}
-
-#endif // CELLTYPETOTYPE_H
diff --git a/tools/vtkMOABReaderNew/DataSetConverter.h b/tools/vtkMOABReaderNew/DataSetConverter.h
index bfae086..3edadf0 100644
--- a/tools/vtkMOABReaderNew/DataSetConverter.h
+++ b/tools/vtkMOABReaderNew/DataSetConverter.h
@@ -2,20 +2,16 @@
#define __smoab_DataSetConverter_h
#include "SimpleMoab.h"
-#include "CellTypeToType.h"
-#include "MixedCellConnectivity.h"
+#include "CellSets.h"
+#include "detail/LoadGeometry.h"
+#include "detail/ReadSparseTag.h"
-#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkDoubleArray.h>
#include <vtkFieldData.h>
#include <vtkIntArray.h>
#include <vtkIdTypeArray.h>
#include <vtkNew.h>
-#include <vtkPointData.h>
-#include <vtkPoints.h>
-#include <vtkUnsignedCharArray.h>
-#include <vtkUnstructuredGrid.h>
#include <algorithm>
@@ -35,17 +31,13 @@ public:
Moab(interface.Moab),
Tag(tag),
ReadMaterialIds(false),
- ReadProperties(false),
- MaterialName("Material")
+ ReadProperties(false)
{
}
void readMaterialIds(bool add) { this->ReadMaterialIds = add; }
bool readMaterialIds() const { return this->ReadMaterialIds; }
- void materialIdName(const std::string& name) { this->MaterialName = name; }
- const std::string& materialIdName() const { return this->MaterialName; }
-
void readProperties(bool readProps) { this->ReadProperties = readProps; }
bool readProperties() const { return this->ReadProperties; }
@@ -54,109 +46,65 @@ public:
//grid. Currently doesn't support reading properties.
//Will read in material ids, if no material id is assigned to an entity,
//its cells will be given an unique id
+ template<typename VTKGridType>
bool fill(const smoab::Range& entities,
- vtkUnstructuredGrid* grid) const
+ VTKGridType* grid) const
{
//create a helper datastructure which can determines all the unique point ids
//and converts moab connecitvity info to vtk connectivity
- moab::Range cells;
- //append all the entities cells together into a single range
+
+ //get all the cells for each parent entity and create
+ // an entity set of those items
+ int dim = this->Tag->value();
typedef smoab::Range::const_iterator iterator;
+ smoab::CellSets entitySets;
for(iterator i=entities.begin(); i!= entities.end(); ++i)
{
+ smoab::Range entitiesCells;
if(this->Tag->isComparable())
{
//if we are comparable only find the cells that match our tags dimension
- smoab::Range entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());
- cells.insert(entitiesCells.begin(),entitiesCells.end());
+ entitiesCells = this->Interface.findEntitiesWithDimension(*i,dim,true);
}
else
{
- //this is a bad representation of all other tags, but we are presuming that
- //neuman and dirichlet are on entitysets with no children
- this->Moab->get_entities_by_handle(*i,cells);
+ entitiesCells = this->Interface.findHighestDimensionEntities(*i,true);
}
+ smoab::CellSet set(*i,entitiesCells);
+ entitySets.push_back(set);
}
- smoab::Range points;
- this->loadCellsAndPoints(cells,points,grid);
-
- if(this->readMaterialIds())
- {
- typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
- typedef std::vector<int>::const_iterator IdConstIterator;
- typedef std::vector<int>::iterator IdIterator;
-
- std::vector<smoab::EntityHandle> searchableCells;
- searchableCells.reserve(grid->GetNumberOfCells());
- std::copy(cells.begin(),cells.end(),std::back_inserter(searchableCells));
- cells.clear(); //release memory we don't need
+ moab::Range cells = smoab::getAllCells(entitySets);
- std::vector<int> materialIds(entities.size());
- //first off iterate the entities and determine which ones
- //have moab material ids
+ //convert the datastructure from a list of cells to a vtk data set
+ detail::LoadGeometry loadGeom(cells,dim,this->Interface);
+ loadGeom.fill(grid);
- //wrap this area with scope, to remove local variables
+ if(this->readMaterialIds())
{
- smoab::MaterialTag tag;
- IdIterator materialIndex = materialIds.begin();
- for(iterator i=entities.begin();
- i != entities.end();
- ++i, ++materialIndex)
- {
- moab::Tag mtag = this->Interface.getMoabTag(tag);
-
- int value=-1;
- this->Moab->tag_get_data(mtag,&(*i),1,&value);
- *materialIndex=static_cast<int>(value);
- }
-
- //now determine ids for all entities that don't have materials
- IdConstIterator maxPos = std::max_element(materialIds.begin(),
- materialIds.end());
- int maxMaterial = *maxPos;
- for(IdIterator i=materialIds.begin(); i!= materialIds.end(); ++i)
- {
- if(*i==-1)
- {
- *i = ++maxMaterial;
- }
- }
- }
- //now we create the material field, and set all the values
- vtkNew<vtkIntArray> materialSet;
- materialSet->SetName(this->materialIdName().c_str());
- materialSet->SetNumberOfValues(grid->GetNumberOfCells());
+ detail::ReadSparseTag materialTagReading(entitySets,
+ cells,
+ this->Interface);
- IdConstIterator materialValue = materialIds.begin();
- for(iterator i=entities.begin(); i!= entities.end(); ++i, ++materialValue)
- {
- //this is a time vs memory trade off, I don't want to store
- //the all the cell ids twice over, lets use more time
- smoab::Range entitiesCells;
- if(this->Tag->isComparable())
- {entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());}
- else
- {this->Moab->get_entities_by_handle(*i,entitiesCells);}
-
- EntityHandleIterator s_begin = searchableCells.begin();
- EntityHandleIterator s_end = searchableCells.end();
- for(iterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
- {
- EntityHandleIterator result = std::lower_bound(s_begin,
- s_end,
- *j);
- std::size_t newId = std::distance(s_begin,result);
- materialSet->SetValue(static_cast<int>(newId), *materialValue);
- }
- }
+ smoab::MaterialTag mtag;
+ vtkNew<vtkIntArray> materials;
+ materials->SetName(mtag.name());
+ materialTagReading.fill(materials.GetPointer(),&mtag);
+ grid->GetCellData()->AddArray(materials.GetPointer());
+ }
- grid->GetCellData()->AddArray(materialSet.GetPointer());
+ //by default we always try to load the default tag
+ detail::ReadSparseTag sTagReading(entitySets,
+ cells,
+ this->Interface);
- }
+ vtkNew<vtkIntArray> sparseTagData;
+ sparseTagData->SetName(this->Tag->name());
+ sTagReading.fill(sparseTagData.GetPointer(),this->Tag);
+ grid->GetCellData()->AddArray(sparseTagData.GetPointer());
return true;
}
@@ -165,28 +113,32 @@ public:
//given a single entity handle create a unstructured grid from it.
//optional third parameter is the material id to use if readMaterialIds
//is on, and no material sparse tag is found for this entity
+ template<typename VTKGridType>
bool fill(const smoab::EntityHandle& entity,
- vtkUnstructuredGrid* grid,
+ VTKGridType* grid,
const int materialId=0) const
{
//create a helper datastructure which can determines all the unique point ids
//and converts moab connecitvity info to vtk connectivity
- moab::Range cells;
+
+ smoab::Range cells;
+ int dim = this->Tag->value();
if(this->Tag->isComparable())
{
//if we are comparable only find the cells that match our tags dimension
- cells = this->Interface.findEntitiesWithDimension(entity,Tag->value());
+ cells = this->Interface.findEntitiesWithDimension(entity,dim,true);
}
else
{
- //this is a bad representation of all other tags, but we are presuming that
- //neuman and dirichlet are on entitysets with no children
- this->Moab->get_entities_by_handle(entity,cells);
+ //load subentities
+ cells = this->Interface.findHighestDimensionEntities(entity,true);
}
- smoab::Range points;
- this->loadCellsAndPoints(cells,points,grid);
+ //convert the datastructure from a list of cells to a vtk data set
+ detail::LoadGeometry loadGeom(cells,dim,this->Interface);
+ loadGeom.fill(grid);
+ const smoab::Range& points = loadGeom.moabPoints();
if(this->readProperties())
{
@@ -194,53 +146,34 @@ public:
this->readProperties(points,grid->GetPointData());
}
+ smoab::CellSets cellSets;
+ smoab::CellSet set(entity,cells);
+ cellSets.push_back(set);
if(this->readMaterialIds())
{
- this->readSparseTag(smoab::MaterialTag(),entity,
- grid->GetNumberOfCells(),
- grid->GetCellData(),
- materialId);
- }
- return true;
- }
+ smoab::MaterialTag mtag;
+ detail::ReadSparseTag materialTagReading(cellSets,
+ cells,
+ this->Interface);
- //----------------------------------------------------------------------------
- void loadCellsAndPoints(const smoab::Range& cells,
- smoab::Range& points,
- vtkUnstructuredGrid* grid) const
- {
+ vtkNew<vtkIntArray> materials;
+ materials->SetName(mtag.name());
+ materialTagReading.fill(materials.GetPointer(),&mtag);
+ grid->GetCellData()->AddArray(materials.GetPointer());
- smoab::MixedCellConnectivity mixConn(cells,this->Moab);
-
- //now that mixConn has all the cells properly stored, lets fixup
- //the ids so that they start at zero and keep the same logical ordering
- //as before.
- vtkIdType numCells, connLen;
- mixConn.compactIds(numCells,connLen);
- this->setGridsTopology(mixConn,grid,numCells,connLen);
+ }
- mixConn.moabPoints(points);
+ //by default we always try to load the default tag
+ detail::ReadSparseTag sTagReading(cellSets,
+ cells,
+ this->Interface);
- vtkNew<vtkPoints> newPoints;
- this->addCoordinates(points,newPoints.GetPointer());
- grid->SetPoints(newPoints.GetPointer());
+ vtkNew<vtkIntArray> sparseTagData;
+ sparseTagData->SetName(this->Tag->name());
+ sTagReading.fill(sparseTagData.GetPointer(),this->Tag);
+ grid->GetCellData()->AddArray(sparseTagData.GetPointer());
- }
-
- //----------------------------------------------------------------------------
- void addCoordinates(smoab::Range pointEntities, vtkPoints* pointContainer) const
- {
- //since the smoab::range are always unique and sorted
- //we can use the more efficient coords_iterate
- //call in moab, which returns moab internal allocated memory
- pointContainer->SetDataTypeToDouble();
- pointContainer->SetNumberOfPoints(pointEntities.size());
-
- //need a pointer to the allocated vtkPoints memory so that we
- //don't need to use an extra copy and we can bypass all vtk's check
- //on out of bounds
- double *rawPoints = static_cast<double*>(pointContainer->GetVoidPointer(0));
- this->Moab->get_coords(pointEntities,rawPoints);
+ return true;
}
private:
@@ -260,36 +193,6 @@ private:
}
//----------------------------------------------------------------------------
- bool readSparseTag(smoab::Tag tag,
- smoab::EntityHandle const& entity,
- vtkIdType length,
- vtkFieldData* field,
- vtkIdType defaultValue) const
- {
-
- typedef std::vector<moab::Tag>::const_iterator iterator;
- moab::Tag mtag = this->Interface.getMoabTag(tag);
-
- int value=0;
- moab::ErrorCode rval = this->Moab->tag_get_data(mtag,&entity,1,&value);
- if(rval!=moab::MB_SUCCESS)
- {
- value = defaultValue;
- }
-
- vtkNew<vtkIntArray> materialSet;
- materialSet->SetNumberOfValues(length);
- materialSet->SetName(this->materialIdName().c_str());
-
- int *raw = static_cast<int*>(materialSet->GetVoidPointer(0));
- std::fill(raw,raw+length,value);
-
- field->AddArray(materialSet.GetPointer());
-
- return true;
- }
-
- //----------------------------------------------------------------------------
void readDenseTags(std::vector<moab::Tag> &tags,
smoab::Range const& entities,
vtkFieldData* field) const
@@ -359,36 +262,7 @@ private:
}
}
- //----------------------------------------------------------------------------
- void setGridsTopology(smoab::MixedCellConnectivity const& mixedCells,
- vtkUnstructuredGrid* grid,
- vtkIdType numCells,
- vtkIdType numConnectivity) const
- {
- //correct the connectivity size to account for the vtk padding
- const vtkIdType vtkConnectivity = numCells + numConnectivity;
- vtkNew<vtkIdTypeArray> cellArray;
- vtkNew<vtkIdTypeArray> cellLocations;
- vtkNew<vtkUnsignedCharArray> cellTypes;
-
- cellArray->SetNumberOfValues(vtkConnectivity);
- cellLocations->SetNumberOfValues(numCells);
- cellTypes->SetNumberOfValues(numCells);
-
- vtkIdType* rawArray = static_cast<vtkIdType*>(cellArray->GetVoidPointer(0));
- vtkIdType* rawLocations = static_cast<vtkIdType*>(cellLocations->GetVoidPointer(0));
- unsigned char* rawTypes = static_cast<unsigned char*>(cellTypes->GetVoidPointer(0));
-
- mixedCells.copyToVtkCellInfo(rawArray,rawLocations,rawTypes);
-
- vtkNew<vtkCellArray> cells;
- cells->SetCells(numCells,cellArray.GetPointer());
- grid->SetCells(cellTypes.GetPointer(),
- cellLocations.GetPointer(),
- cells.GetPointer(),
- NULL,NULL);
- }
};
}
diff --git a/tools/vtkMOABReaderNew/ExtractShell.h b/tools/vtkMOABReaderNew/ExtractShell.h
new file mode 100644
index 0000000..359fbce
--- /dev/null
+++ b/tools/vtkMOABReaderNew/ExtractShell.h
@@ -0,0 +1,80 @@
+#ifndef __smoab_ExtractShell_h
+#define __smoab_ExtractShell_h
+
+#include "SimpleMoab.h"
+#include "detail/UsageTable.h"
+
+#include <algorithm>
+
+namespace smoab{
+
+class ExtractShell
+{
+ const smoab::Interface& Interface;
+ smoab::CellSets VCells;
+
+public:
+ ExtractShell(const smoab::CellSets volCells,
+ const smoab::Interface& interface):
+ Interface(interface),
+ VCells(volCells)
+ {
+ }
+
+ bool findSkins(smoab::CellSets &surfaceCellSets);
+};
+
+
+//----------------------------------------------------------------------------
+bool ExtractShell::findSkins(smoab::CellSets &surfaceCellSets)
+{
+ typedef smoab::Range::const_iterator Iterator;
+
+ typedef smoab::CellSets::const_iterator SetIterator;
+
+
+ smoab::Range cellsToRemove;
+ for(SetIterator set = this->VCells.begin();
+ set != this->VCells.end();
+ ++set)
+ {
+ const smoab::Range &cells = set->cells();
+ this->Interface.createAdjacencies(set->cells(),2);
+
+
+ //we create the usage table for each iteration so that we only
+ //get the shell of each cell set. If we used the table between
+ //sets we would get the shell of the combined sets
+ smoab::detail::UsageTable table;
+ for(Iterator i = cells.begin(); i != cells.end(); ++i)
+ {
+ std::vector<smoab::EntityHandle> faceCells =
+ this->Interface.sideElements(*i,2);
+
+ //the usage id allows you to label cells when going into the table
+ //so that you can extract multiple shells where each is based on
+ //a single region id.
+ std::vector<int> regionId(1,faceCells.size());
+ table.incrementUsage(faceCells,regionId);
+ }
+ smoab::Range surfaceCells = table.singleUsage();
+
+ //create a new cell set that
+ smoab::CellSet surfaceSet(set->entity(),surfaceCells);
+ surfaceCellSets.push_back(surfaceSet);
+
+ smoab::Range subsetToRemove = table.multipleUsage();
+ cellsToRemove.insert(subsetToRemove.begin(),subsetToRemove.end());
+ }
+
+ //we will remove all cells that have multiple usages from the moab database
+ //I really don't care if they already existed or not.
+ this->Interface.remove(cellsToRemove);
+ return true;
+}
+
+
+
+}
+
+#endif // __smoab_ExtractShell_h
diff --git a/tools/vtkMOABReaderNew/FaceSets.h b/tools/vtkMOABReaderNew/FaceSets.h
new file mode 100644
index 0000000..536f3f9
--- /dev/null
+++ b/tools/vtkMOABReaderNew/FaceSets.h
@@ -0,0 +1,256 @@
+#ifndef __smoab_FaceSets_h
+#define __smoab_FaceSets_h
+
+#include "CellSets.h"
+#include <set>
+
+namespace smoab
+{
+//----------------------------------------------------------------------------
+class FaceCellSet : public CellSet
+{
+public:
+ FaceCellSet(int id, smoab::EntityHandle p,const smoab::Range& cells):
+ CellSet(p,cells),
+ ID(id)
+ {}
+
+ int faceId() const { return ID; }
+ void overrideFaceId(int i) { ID = i; } //USE AT YOUR OWN RISK
+
+private:
+ int ID;
+};
+
+
+//----------------------------------------------------------------------------
+typedef std::vector<FaceCellSet> FaceCellSets;
+
+
+//----------------------------------------------------------------------------
+//class that store the regions that a faces is adjacent too
+struct FacesAdjRegions
+{
+ FacesAdjRegions(int f, smoab::EntityHandle r0, smoab::EntityHandle r1):
+ FaceId(f),
+ Region0(r0),
+ Region1(r1)
+ {
+ if (r0 > r1)
+ {
+ std::swap(this->Region0,this->Region1);
+ }
+ }
+
+ FacesAdjRegions(int f):
+ FaceId(f),
+ Region0(-3),
+ Region1(-2)
+ {}
+
+ bool operator<(const FacesAdjRegions& other) const
+ {
+ return (this->FaceId < other.FaceId);
+ }
+
+ smoab::EntityHandle otherId(smoab::EntityHandle other) const
+ {
+ if(other == Region0)
+ {
+ return Region1;
+ }
+ return Region0;
+ }
+
+ int FaceId;
+ smoab::EntityHandle Region0;
+ smoab::EntityHandle Region1;
+};
+//----------------------------------------------------------------------------
+smoab::FaceCellSets findFaceSets(smoab::CellSets shells,
+ smoab::CellSets boundaries,
+ std::set<smoab::FacesAdjRegions>& faceMaps)
+{
+ typedef smoab::CellSets::iterator iterator;
+ typedef smoab::FaceCellSets::iterator faceIterator;
+ typedef std::set<smoab::FacesAdjRegions>::const_iterator FaceAdjIterator;
+
+ //we need to properly label each unique face in shells
+ //we do this by intersecting each shell with each other shell
+ //to find shell on shell contact, and than we intersect each
+ //resulting shell with the boundary conditions
+ //the end result of these intersections will be the new modelfaces
+ int faceId = 1;
+ smoab::FaceCellSets shellFaces;
+
+ //first intersect each shell with each other shell
+ std::set<smoab::FacesAdjRegions> shellFaceContacts;
+ for(iterator i=shells.begin();i!= shells.end();++i)
+ {
+ //copy the cells so we can add a face that represents
+ //all the cells of the region that aren't shared with another region
+ int numCells = i->cells().size(); //size() on range is slow, so cache it
+ for(iterator j = i+1;
+ j != shells.end() && numCells > 0;
+ ++j)
+ {
+ //intersect i and j to make a new face
+ smoab::Range intersection = smoab::intersect(i->cells(),j->cells());
+ if(!intersection.empty())
+ {
+ //don't want to increment faceId when the intersection is empty
+ smoab::FaceCellSet face(faceId++,i->entity(),intersection);
+ shellFaces.push_back(face);
+ i->erase(intersection);
+ j->erase(intersection);
+ numCells -= intersection.size();
+
+ //add this to the face map
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),j->entity());
+ shellFaceContacts.insert(faceInfo);
+ }
+ }
+ //if all the cells for shell i are used, don't add a new
+ //empty face
+ if(numCells > 0)
+ {
+ smoab::FaceCellSet face(faceId++,i->entity(),i->cells());
+ shellFaces.push_back(face);
+
+ //add this to the face map
+ smoab::FacesAdjRegions faceInfo(faceId-1,-1,i->entity());
+ shellFaceContacts.insert(faceInfo);
+ }
+ }
+
+ //now we have all the faces that match shell on shell contact
+ //we know process all the new faces to see if they intersect
+ //with any boundary sets. A boundary set can span multiple
+ //shells so we want to process it as a second loop
+
+ //store the end before we start adding boundary faces, which
+ //we don't need to check agianst other boundaries
+ faceId = 1; //reset the faced id
+
+ //store in a new face set, expanding the current one causes incorrect results
+ smoab::FaceCellSets faces;
+ for(faceIterator i=shellFaces.begin();i != shellFaces.end();++i)
+ {
+ //determine from the shell faces if the new face we are creating
+ //is bounded by two regions or just one
+ smoab::FacesAdjRegions idToSearchFor(i->faceId());
+ FaceAdjIterator adjRegions = shellFaceContacts.find(idToSearchFor);
+ smoab::EntityHandle otherRegionId = adjRegions->otherId(i->entity());
+
+ int numCells = i->cells().size(); //size() on range is slow, so cache it
+ for(iterator j=boundaries.begin();j != boundaries.end(); ++j)
+ {
+ smoab::Range intersect = smoab::intersect(i->cells(),j->cells());
+ if(!intersect.empty())
+ {
+ //don't want to increment faceId when the intersection is empty
+ smoab::FaceCellSet face(faceId++,j->entity(),intersect);
+ faces.push_back(face);
+ i->erase(intersect);
+ numCells -= intersect.size();
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),otherRegionId);
+ faceMaps.insert(faceInfo);
+ }
+ }
+ if(numCells > 0)
+ {
+ smoab::FaceCellSet face(faceId++,i->entity(),i->cells());
+ faces.push_back(face);
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),otherRegionId);
+ faceMaps.insert(faceInfo);
+ }
+ }
+ return faces;
+}
+
+//----------------------------------------------------------------------------
+template<typename T>
+std::vector<T> faceIdsPerCell(const smoab::FaceCellSets& faces)
+{
+ typedef smoab::FaceCellSets::const_iterator iterator;
+ typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
+ typedef smoab::Range::const_iterator RangeIterator;
+
+ //find all the cells that are in the faceCellSet, and than map
+ //the proper face id to that relative position, here comes lower_bounds!
+ std::vector<smoab::EntityHandle> searchableCells;
+ smoab::Range faceRange = smoab::getAllCells(faces);
+ smoab::RangeToVector(faceRange,searchableCells);
+
+ //faceIds will be the resulting array
+ std::vector<T> faceIds(searchableCells.size());
+
+ //construct the start and end iterators for the lower bounds call
+
+ EntityHandleIterator s_begin = searchableCells.begin();
+ EntityHandleIterator s_end = searchableCells.end();
+
+ //search the face cell sets
+ for(iterator i=faces.begin(); i!=faces.end(); ++i)
+ {
+ T value = static_cast<T>(i->faceId());
+ const smoab::Range& entitiesCells = i->cells();
+ for(RangeIterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
+ {
+ EntityHandleIterator result = std::lower_bound(s_begin,
+ s_end,
+ *j);
+ std::size_t newId = std::distance(s_begin,result);
+ faceIds[newId] = value;
+ }
+ }
+ return faceIds;
+}
+
+//----------------------------------------------------------------------------
+//given a face adjacency, determine the regions spare tag values
+template<typename T>
+std::pair<T,T> FaceAdjRegionValues(const smoab::FacesAdjRegions& faceAdj,
+ smoab::Tag* t,
+ const smoab::Interface& interface)
+ {
+ std::pair<T,T> returnValue;
+ const int defaultValue = interface.getDefaultTagVaue<int>(*t);
+
+ /*
+ * IF A REGION IS SET TO -1 WE NEED TO PUSH THAT VALUE DOWN
+ * AS THE MATERIAL, SINCE THE MOAB DEFAULT TAG VALUE WILL
+ * BE CONSIDIERED A REGION, AND WE WANT TO SAY IT BOUNDS THE
+ * VOID REGION
+ */
+ int tagValue = defaultValue; //use tagValue to pass in default value
+ if(faceAdj.Region0 != -1)
+ {
+ tagValue = interface.getTagData(*t,faceAdj.Region0,tagValue);
+ }
+ else
+ {
+ tagValue = -1;
+ }
+
+ //set the first region tag value into the pair we are returing
+ returnValue.first = static_cast<T>(tagValue);
+
+ tagValue = defaultValue; //use tagValue to pass in default value
+ tagValue = interface.getTagData(*t,faceAdj.Region1,tagValue);
+ if(faceAdj.Region1 != -1)
+ {
+ tagValue = interface.getTagData(*t,faceAdj.Region1,tagValue);
+ }
+ else
+ {
+ tagValue = -1;
+ }
+ returnValue.second = static_cast<T>(tagValue);
+
+ return returnValue;
+ }
+
+ } //smoab
+
+#endif
diff --git a/tools/vtkMOABReaderNew/MixedCellConnectivity.h b/tools/vtkMOABReaderNew/MixedCellConnectivity.h
deleted file mode 100644
index 1b26d0f..0000000
--- a/tools/vtkMOABReaderNew/MixedCellConnectivity.h
+++ /dev/null
@@ -1,271 +0,0 @@
-#ifndef __smoab_MixedCellConnectivity_h
-#define __smoab_MixedCellConnectivity_h
-
-#include "vtkCellType.h"
-#include <algorithm>
-
-namespace
-{
-
-template<int N> struct QuadratricOrdering{};
-
-template<> struct QuadratricOrdering<VTK_QUADRATIC_WEDGE>
-{
- static const int NUM_VERTS = 15;
- void reorder(vtkIdType* connectivity) const
- {
- std::swap_ranges(connectivity+9,connectivity+12,connectivity+12);
- }
-};
-
-template<> struct QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON>
-{
- static const int NUM_VERTS = 27;
- void reorder(vtkIdType* connectivity) const
- {
- std::swap_ranges(connectivity+12,connectivity+16,connectivity+16);
-
- //move 20 to 22
- //move 22 to 23
- //move 23 to 20
-
- //swap 20 with 22
- std::swap(connectivity[20],connectivity[23]);
-
- //swap 22 with 23
- std::swap(connectivity[22],connectivity[23]);
- }
-};
-
-template<typename QuadraticOrdering>
-void FixQuadraticIdOrdering(vtkIdType* connectivity, vtkIdType numCells,
- QuadraticOrdering& ordering)
-{
- //skip the first index that holds the length of the cells
- //if we skip it once here, and than properly increment it makes the code
- //far easier
- connectivity+=1;
- for(vtkIdType i=0; i < numCells; ++i)
- {
- ordering.reorder(connectivity);
- connectivity += ordering.NUM_VERTS + 1;
- }
-}
-}
-
-namespace smoab
-{
-
-class MixedCellConnectivity
-{
-public:
- MixedCellConnectivity(smoab::Range const& cells, moab::Interface* moab):
- Connectivity(),
- UniquePoints(),
- Info()
- {
- int count = 0;
- const std::size_t cellSize=cells.size();
- while(count != cellSize)
- {
- EntityHandle* connectivity;
- int numVerts=0, iterationCount=0;
- //use the highly efficent calls, since we know that are of the same dimension
- moab->connect_iterate(cells.begin()+count,
- cells.end(),
- connectivity,
- numVerts,
- iterationCount);
- //if we didn't read anything, break!
- if(iterationCount == 0)
- {
- break;
- }
-
- //identify the cell type that we currently have,
- //store that along with the connectivity in a temp storage vector
- const moab::EntityType type = moab->type_from_handle(*cells.begin()+count);
-
- //while all these cells are contiously of the same type,
- //quadric hexs in vtk have 20 points, but moab has 21 so we
- //need to store this difference
- int numVTKVerts = numVerts;
- int vtkCellType = smoab::vtkCellType(type,numVTKVerts);
-
- RunLengthInfo info = { vtkCellType, numVerts, (numVerts-numVTKVerts), iterationCount };
- this->Info.push_back(info);
- this->Connectivity.push_back(connectivity);
-
- count += iterationCount;
- }
- }
-
- //----------------------------------------------------------------------------
- void compactIds(vtkIdType& numCells, vtkIdType& connectivityLength)
- {
- //converts all the ids to be ordered starting at zero, and also
- //keeping the orginal logical ordering. Stores the result of this
- //operation in the unstrucutred grid that is passed in
-
- //lets determine the total length of the connectivity
- connectivityLength = 0;
- numCells = 0;
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i)
- {
- connectivityLength += (*i).numCells * (*i).numVerts;
- numCells += (*i).numCells;
- }
-
- this->UniquePoints.reserve(connectivityLength);
-
- this->copyConnectivity(this->UniquePoints);
- std::sort(this->UniquePoints.begin(),this->UniquePoints.end());
-
- typedef std::vector<EntityHandle>::iterator EntityIterator;
- EntityIterator newEnd = std::unique(this->UniquePoints.begin(),
- this->UniquePoints.end());
-
- const std::size_t newSize = std::distance(this->UniquePoints.begin(),newEnd);
- this->UniquePoints.resize(newSize);
- }
-
- //----------------------------------------------------------------------------
- void moabPoints(smoab::Range& range) const
- {
- //from the documentation a reverse iterator is the fastest way
- //to insert into a range.
- std::copy(this->UniquePoints.rbegin(),
- this->UniquePoints.rend(),
- moab::range_inserter(range));
- }
-
- //----------------------------------------------------------------------------
- //copy the connectivity from the moab held arrays to the user input vector
- void copyConnectivity(std::vector<EntityHandle>& output) const
- {
- //walk the info to find the length of each sub connectivity array,
- //and insert them into the vector, ordering is implied by the order
- //the connecitivy sub array are added to this class
- ConnConstIterator c = this->Connectivity.begin();
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i,++c)
- {
- //remember our Connectivity is a vector of pointers whose
- //length is held in the info vector.
- const int numUnusedPoints = (*i).numUnusedVerts;
- if(numUnusedPoints==0)
- {
- const int connLength = (*i).numCells * (*i).numVerts;
- std::copy(*c,*c+connLength,std::back_inserter(output));
- }
- else
- {
- //we have cell connectivity that we need to skip,
- //so we have to manual copy each cells connectivity
- const int size = (*i).numCells;
- const int numPoints = (*i).numVerts;
- for(int i=0; i < size; ++i)
- {
- std::copy(*c,*c+numPoints,std::back_inserter(output));
- }
- c+=numPoints + (*i).numUnusedVerts;
- }
-
- }
- }
-
- //copy the information from this contianer to a vtk cell array, and
- //related lookup information
- void copyToVtkCellInfo(vtkIdType* cellArray,
- vtkIdType* cellLocations,
- unsigned char* cellTypes) const
- {
- vtkIdType currentVtkConnectivityIndex = 0;
- ConnConstIterator c = this->Connectivity.begin();
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i, ++c)
- {
- //for this group of the same cell type we need to fill the cellTypes
- const int numCells = (*i).numCells;
- const int numVerts = (*i).numVerts;
-
- std::fill_n(cellTypes,
- numCells,
- static_cast<unsigned char>((*i).type));
-
- //for each cell in this collection that have the same type
- //grab the raw array now, so we can properly increment for each vert in each cell
- EntityHandle* moabConnectivity = *c;
- for(int j=0;j < numCells; ++j)
- {
- cellLocations[j]= currentVtkConnectivityIndex;
-
- //cell arrays start and end are different, since we
- //have to account for element that states the length of each cell
- cellArray[0]=numVerts;
-
- for(int k=0; k < numVerts; ++k, ++moabConnectivity )
- {
- //this is going to be a root of some failures when we start
- //reading really large datasets under 32bit.
-
-
- //fyi, don't use a range ds for unique points, distance
- //function is horribly slow they need to override it
- EntityConstIterator result = std::lower_bound(
- this->UniquePoints.begin(),
- this->UniquePoints.end(),
- *moabConnectivity);
- std::size_t newId = std::distance(this->UniquePoints.begin(),
- result);
- cellArray[k+1] = static_cast<vtkIdType>(newId);
- }
-
- //skip any extra unused points, which is currnetly only
- //the extra center point in moab quadratic hex
- moabConnectivity+=(*i).numUnusedVerts;
-
- currentVtkConnectivityIndex += numVerts+1;
- cellArray += numVerts+1;
- }
-
- //For Tri-Quadratic-Hex and Quadratric-Wedge Moab and VTK
- //Differ on the order of the edge ids. For wedge we need to swap
- //indices 9,10,11 with 12,13,14 for each cell. For Hex we sawp
- //12,13,14,15 with 16,17,18,19
- int vtkCellType = (*i).type;
- vtkIdType* connectivity = cellArray - (numCells * (numVerts+1));
- if(vtkCellType == VTK_TRIQUADRATIC_HEXAHEDRON)
- {
- ::QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON> newOrdering;
- ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
- }
- else if(vtkCellType == VTK_QUADRATIC_WEDGE)
- {
- ::QuadratricOrdering<VTK_QUADRATIC_WEDGE> newOrdering;
- ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
- }
-
- cellLocations += numCells;
- cellTypes += numCells;
- }
-
- }
-
-private:
- std::vector<EntityHandle*> Connectivity;
- std::vector<EntityHandle> UniquePoints;
-
- struct RunLengthInfo{ int type; int numVerts; int numUnusedVerts; int numCells; };
- std::vector<RunLengthInfo> Info;
-
- typedef std::vector<EntityHandle>::const_iterator EntityConstIterator;
- typedef std::vector<EntityHandle*>::const_iterator ConnConstIterator;
- typedef std::vector<RunLengthInfo>::const_iterator InfoConstIterator;
-};
-}
-#endif // __smoab_MixedCellConnectivity_h
diff --git a/tools/vtkMOABReaderNew/SimpleMoab.h b/tools/vtkMOABReaderNew/SimpleMoab.h
index 8fc54e6..4fb734c 100644
--- a/tools/vtkMOABReaderNew/SimpleMoab.h
+++ b/tools/vtkMOABReaderNew/SimpleMoab.h
@@ -5,6 +5,8 @@
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
+#include "moab/CN.hpp"
+
#include "MBTagConventions.hpp"
#include <iostream>
@@ -26,12 +28,6 @@ using moab::intersect;
using moab::subtract;
using moab::unite;
-//forward declare this->Moab for Tag
-struct Interface;
-
-//forward declar the DataSetConverter so it can be a friend of Interface
-class DataSetConverter;
-
class Tag
{
const std::string Name_;
@@ -66,10 +62,24 @@ public:
GeomTag(int d):Tag("GEOM_DIMENSION"),dim(d){}
GeomTag():Tag("GEOM_DIMENSION"), dim(0){}
+ virtual ~GeomTag(){}
+
bool isComparable() const { return dim > 0; }
int value() const { return dim; }
};
+
+//forward declare this->Moab for Tag
+struct Interface;
+
+//forward declare the DataSetConverter so it can be a friend of Interface
+class DataSetConverter;
+
+//forward declare the LoadGeometry so it can be a friend of Interface
+namespace detail{ class LoadGeometry; }
+namespace detail{ class LoadPoly; }
+
+
//light weight wrapper on a moab this->Moab that exposes only the reduced class
//that we need
class Interface
@@ -102,6 +112,37 @@ public:
}
//----------------------------------------------------------------------------
+ template<typename T>
+ T getDefaultTagVaue(moab::Tag tag) const
+ {
+ T defaultValue;
+ this->Moab->tag_get_default_value(tag,&defaultValue);
+ return defaultValue;
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getDefaultTagVaue(smoab::Tag tag) const
+ {
+ return this->getDefaultTagVaue<T>(getMoabTag(tag));
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getTagData(moab::Tag tag, const smoab::EntityHandle& entity, T value) const
+ {
+ this->Moab->tag_get_data(tag,&entity,1,&value);
+ return value;
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getTagData(smoab::Tag tag, const smoab::EntityHandle& entity, T value = T()) const
+ {
+ return this->getTagData(getMoabTag(tag),entity,value);
+ }
+
+ //----------------------------------------------------------------------------
//returns the moab name for the given entity handle if it has a sparse Name tag
std::string name(const smoab::EntityHandle& entity) const
{
@@ -119,6 +160,19 @@ public:
return std::string(name);
}
+ //----------------------------------------------------------------------------
+ //returns the geometeric dimension of an entity.
+ int dimension(const smoab::EntityHandle& entity) const
+ {
+ return this->Moab->dimension_from_handle(entity);
+ }
+
+ //----------------------------------------------------------------------------
+ //returns the geometeric dimension of an entity.
+ smoab::EntityType entityType(const smoab::EntityHandle& entity) const
+ {
+ return this->Moab->type_from_handle(entity);
+ }
//----------------------------------------------------------------------------
smoab::EntityHandle getRoot() const { return this->Moab->get_root_set(); }
@@ -133,6 +187,17 @@ public:
}
//----------------------------------------------------------------------------
+ //given a single entity handle find all items in that mesh set that aren't
+ //them selves entitysets. If recurse is true we also recurse sub entitysets
+ smoab::Range findAllMeshEntities(smoab::EntityHandle const& entity,
+ bool recurse=false) const
+ {
+ smoab::Range result;
+ this->Moab->get_entities_by_handle(entity,result,recurse);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
//Find all entities with a given tag. We don't use geom as a tag as that
//isn't a fast operation. Yes finding the intersection of geom entities and
//a material / boundary tag will be more work, but it is rarely done currently
@@ -179,60 +244,55 @@ public:
//----------------------------------------------------------------------------
//Find all entities from a given root of a given dimensionality
smoab::Range findEntitiesWithDimension(const smoab::EntityHandle root,
- int dimension) const
+ const int dimension,
+ bool recurse=false) const
{
typedef smoab::Range::const_iterator iterator;
smoab::Range result;
- this->Moab->get_entities_by_dimension(root,dimension,result);
+ this->Moab->get_entities_by_dimension(root,dimension,result,recurse);
-
- smoab::Range children;
- this->Moab->get_child_meshsets(root,children,0);
- for(iterator i=children.begin(); i !=children.end();++i)
+ if(recurse)
{
- this->Moab->get_entities_by_dimension(*i,dimension,result);
+ smoab::Range children;
+ this->Moab->get_child_meshsets(root,children,0);
+ for(iterator i=children.begin(); i !=children.end();++i)
+ {
+ this->Moab->get_entities_by_dimension(*i,dimension,result);
+ }
}
return result;
}
- //----------------------------------------------------------------------------
- smoab::Range findAdjacentEntities(const smoab::EntityHandle& entity,
- int dimension) const
- {
- const int adjType = static_cast<int>(smoab::INTERSECT);
- smoab::Range result;
- const bool create_if_missing = false;
- this->Moab->get_adjacencies(&entity,
- 1,
- dimension,
- create_if_missing,
- result,
- adjType);
- return result;
- }
- //----------------------------------------------------------------------------
- smoab::Range findAdjacentEntities(const smoab::Range& range,
- int dimension,
- const smoab::adjacency_type type = smoab::UNION) const
+ //----------------------------------------------------------------------------
+ smoab::Range findHighestDimensionEntities(const smoab::EntityHandle& entity,
+ bool recurse=false) const
{
- //the smoab and moab adjacent intersection enums are in the same order
- const int adjType = static_cast<int>(type);
- smoab::Range result;
- const bool create_if_missing = false;
- this->Moab->get_adjacencies(range,dimension,
- create_if_missing,
- result,
- adjType);
+ //the goal is to load all entities that are not entity sets of this
+ //node, while also subsetting by the highest dimension
- return result;
+ //lets find the entities of only the highest dimension
+ int num_ents=0;
+ int dim=3;
+ while(num_ents<=0&&dim>0)
+ {
+ this->Moab->get_number_entities_by_dimension(entity,dim,num_ents,recurse);
+ --dim;
+ }
+ ++dim; //reincrement to correct last decrement
+ if(num_ents > 0)
+ {
+ //we have found entities of a given dimension
+ return this->findEntitiesWithDimension(entity,dim,recurse);
+ }
+ return smoab::Range();
}
//----------------------------------------------------------------------------
//Find all elements in the database that have children and zero parents.
//this doesn't find
- smoab::Range findEntityRootParents(smoab::EntityHandle const& root) const
+ smoab::Range findEntityRootParents(const smoab::EntityHandle& root) const
{
smoab::Range parents;
@@ -258,7 +318,7 @@ public:
//----------------------------------------------------------------------------
//finds entities that have zero children and zero parents
- smoab::Range findDetachedEntities(moab::EntityHandle const& root) const
+ smoab::Range findDetachedEntities(const moab::EntityHandle& root) const
{
smoab::Range detached;
@@ -284,7 +344,7 @@ public:
//----------------------------------------------------------------------------
//find all children of the entity passed in that has multiple parents
- smoab::Range findEntitiesWithMultipleParents(smoab::EntityHandle const& root)
+ smoab::Range findEntitiesWithMultipleParents(const smoab::EntityHandle& root) const
{
smoab::Range multipleParents;
typedef moab::Range::const_iterator iterator;
@@ -305,8 +365,114 @@ public:
}
//----------------------------------------------------------------------------
+ //find all entities that are adjacent to a single entity
+ smoab::Range findAdjacencies(const smoab::EntityHandle& entity,
+ int dimension) const
+ {
+ const int adjType = static_cast<int>(smoab::INTERSECT);
+ smoab::Range result;
+ const bool create_if_missing = false;
+ this->Moab->get_adjacencies(&entity,
+ 1,
+ dimension,
+ create_if_missing,
+ result,
+ adjType);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ smoab::Range findAdjacencies(const smoab::Range& range,
+ int dimension,
+ const smoab::adjacency_type type = smoab::UNION) const
+ {
+ //the smoab and moab adjacent intersection enums are in the same order
+ const int adjType = static_cast<int>(type);
+ smoab::Range result;
+ const bool create_if_missing = false;
+ this->Moab->get_adjacencies(range,dimension,
+ create_if_missing,
+ result,
+ adjType);
+
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ //create adjacencies, only works when the dimension requested is lower than
+ //dimension of the range of entities
+ smoab::Range createAdjacencies(const smoab::Range& range,
+ int dimension,
+ const smoab::adjacency_type type = smoab::UNION) const
+ {
+ //the smoab and moab adjacent intersection enums are in the same order
+ const int adjType = static_cast<int>(type);
+ smoab::Range result;
+ const bool create_if_missing = true;
+ this->Moab->get_adjacencies(range,dimension,
+ create_if_missing,
+ result,
+ adjType);
+
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ int numChildMeshSets(const smoab::EntityHandle& root) const
+ {
+ int numChildren;
+ this->Moab->num_child_meshsets(root,&numChildren);
+ return numChildren;
+ }
+
+ //----------------------------------------------------------------------------
+ smoab::Range getChildSets(const smoab::EntityHandle& root) const
+ {
+ smoab::Range children;
+ this->Moab->get_child_meshsets(root,children,0);
+ return children;
+ }
+
+ //----------------------------------------------------------------------------
+ //remove a collection of entities from the database
+ void remove(smoab::Range const& toDelete) const
+ {
+ this->Moab->delete_entities(toDelete);
+ }
+
+ //----------------------------------------------------------------------------
+ //a entityHandle with value zero means no side element was found
+ smoab::EntityHandle sideElement(smoab::EntityHandle const& cell,
+ int dim, int side) const
+ {
+ smoab::EntityHandle result(0);
+ this->Moab->side_element(cell,dim,side,result);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ //returns all the existing side elements of a cell, elements that
+ //are zero mean that side element doesn't exist
+ std::vector<smoab::EntityHandle> sideElements(
+ smoab::EntityHandle const& cell,
+ int dim) const
+ {
+ const EntityType volumeCellType = this->Moab->type_from_handle(cell);
+ const int numSides = static_cast<int>(moab::CN::NumSubEntities(
+ volumeCellType, dim));
+
+ std::vector<smoab::EntityHandle> result(numSides);
+ for (int side = 0; side < numSides; ++side)
+ {
+ smoab::EntityHandle *sideElem = &result[side]; //get memory of vector
+ this->Moab->side_element(cell,dim,side,*sideElem);
+ }
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
//prints all elements in a range objects
- void printRange(smoab::Range const& range)
+ void printRange(const smoab::Range& range) const
{
typedef Range::const_iterator iterator;
for(iterator i=range.begin(); i!=range.end(); ++i)
@@ -315,12 +481,25 @@ public:
this->Moab->list_entity(*i);
}
}
-
friend class smoab::DataSetConverter;
+ friend class smoab::detail::LoadGeometry;
+ friend class smoab::detail::LoadPoly;
private:
moab::Interface* Moab;
};
+//----------------------------------------------------------------------------
+void RangeToVector(const smoab::Range &range,
+ std::vector<smoab::EntityHandle>& vector )
+{
+ vector.reserve(range.size());
+ std::copy(range.begin(),
+ range.end(),
+ std::back_inserter(vector));
+}
+
+
+
}
#endif
diff --git a/tools/vtkMOABReaderNew/detail/CellTypeToType.h b/tools/vtkMOABReaderNew/detail/CellTypeToType.h
new file mode 100644
index 0000000..b2955e4
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/CellTypeToType.h
@@ -0,0 +1,147 @@
+#ifndef __smoab_detail_CellTypeToType_h
+#define __smoab_detail_CellTypeToType_h
+
+#include "vtkCellType.h"
+#include <algorithm>
+
+namespace smoab{ namespace detail{
+
+template<int N> struct QuadratricOrdering{};
+
+template<> struct QuadratricOrdering<VTK_QUADRATIC_WEDGE>
+{
+ static const int NUM_VERTS = 15;
+ void reorder(vtkIdType* connectivity) const
+ {
+ std::swap_ranges(connectivity+9,connectivity+12,connectivity+12);
+ }
+};
+
+template<> struct QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON>
+{
+ static const int NUM_VERTS = 27;
+ void reorder(vtkIdType* connectivity) const
+ {
+ std::swap_ranges(connectivity+12,connectivity+16,connectivity+16);
+
+ //move 20 to 22
+ //move 22 to 23
+ //move 23 to 20
+
+ //swap 20 with 22
+ std::swap(connectivity[20],connectivity[23]);
+
+ //swap 22 with 23
+ std::swap(connectivity[22],connectivity[23]);
+ }
+};
+
+template<typename QuadraticOrdering>
+void FixQuadraticIdOrdering(vtkIdType* connectivity, vtkIdType numCells,
+ QuadraticOrdering& ordering)
+{
+ //skip the first index that holds the length of the cells
+ //if we skip it once here, and than properly increment it makes the code
+ //far easier
+ connectivity+=1;
+ for(vtkIdType i=0; i < numCells; ++i)
+ {
+ ordering.reorder(connectivity);
+ connectivity += ordering.NUM_VERTS + 1;
+ }
+}
+
+
+int vtkCellType(moab::EntityType t, int &num_connect)
+ {
+ int ctype = -1;
+ switch (t)
+ {
+ case moab::MBEDGE:
+ if (num_connect == 2) ctype = VTK_LINE;
+ else if (num_connect == 3) ctype = VTK_QUADRATIC_EDGE;
+ break;
+ case moab::MBTRI:
+ if (num_connect == 3) ctype = VTK_TRIANGLE;
+ else if (num_connect == 6) ctype = VTK_QUADRATIC_TRIANGLE;
+ else if (num_connect == 7) ctype = VTK_BIQUADRATIC_TRIANGLE;
+ break;
+ case moab::MBQUAD:
+ if (num_connect == 4) ctype = VTK_QUAD;
+ else if (num_connect == 8) ctype = VTK_QUADRATIC_QUAD;
+ else if (num_connect == 9) ctype = VTK_BIQUADRATIC_QUAD;
+ break;
+ case moab::MBPOLYGON:
+ if (num_connect == 4) ctype = VTK_POLYGON;
+ break;
+ case moab::MBTET:
+ if (num_connect == 4) ctype = VTK_TETRA;
+ else if (num_connect == 10) ctype = VTK_QUADRATIC_TETRA;
+ break;
+ case moab::MBPYRAMID:
+ if (num_connect == 5) ctype = VTK_PYRAMID;
+ else if (num_connect == 13) ctype = VTK_QUADRATIC_PYRAMID;
+ break;
+ case moab::MBPRISM:
+ if (num_connect == 6) ctype = VTK_WEDGE;
+ else if (num_connect == 15) ctype = VTK_QUADRATIC_WEDGE;
+ break;
+ case moab::MBHEX:
+ if (num_connect == 8) ctype = VTK_HEXAHEDRON;
+ else if (num_connect == 20) ctype = VTK_QUADRATIC_HEXAHEDRON;
+ else if (num_connect == 21) ctype = VTK_QUADRATIC_HEXAHEDRON, num_connect = 20;
+ else if (num_connect == 27) ctype = VTK_TRIQUADRATIC_HEXAHEDRON;
+ break;
+ default:
+ ctype = -1;
+ break;
+ }
+ return ctype;
+ }
+
+int vtkLinearCellType(moab::EntityType t, int &num_connect)
+ {
+ int ctype = -1;
+ switch (t)
+ {
+ case moab::MBEDGE:
+ ctype = VTK_LINE;
+ num_connect = 2;
+ break;
+ case moab::MBTRI:
+ ctype = VTK_TRIANGLE;
+ num_connect = 3;
+ break;
+ case moab::MBQUAD:
+ ctype = VTK_QUAD;
+ num_connect = 4;
+ break;
+ case moab::MBPOLYGON:
+ ctype = VTK_POLYGON;
+ num_connect = 4;
+ break;
+ case moab::MBTET:
+ ctype = VTK_TETRA;
+ num_connect = 4;
+ break;
+ case moab::MBPYRAMID:
+ ctype = VTK_PYRAMID;
+ num_connect = 5;
+ break;
+ case moab::MBPRISM:
+ ctype = VTK_WEDGE;
+ num_connect = 6;
+ break;
+ case moab::MBHEX:
+ ctype = VTK_HEXAHEDRON;
+ num_connect = 8;
+ break;
+ default:
+ break;
+ }
+ return ctype;
+ }
+
+} } //namespace smaob::detail
+
+#endif // CELLTYPETOTYPE_H
diff --git a/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h b/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h
new file mode 100644
index 0000000..28319df
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h
@@ -0,0 +1,18 @@
+#ifndef __smoab_detail_ContinousCellInfo_h
+#define __smoab_detail_ContinousCellInfo_h
+
+
+namespace smoab { namespace detail {
+
+struct ContinousCellInfo
+{
+ int type;
+ int numVerts;
+ int numUnusedVerts;
+ int numCells;
+
+};
+
+} } //namespace smoab::detail
+
+#endif
diff --git a/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h b/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h
new file mode 100644
index 0000000..1b94548
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h
@@ -0,0 +1,219 @@
+#ifndef __smoab_LinearCellConnectivity_h
+#define __smoab_LinearCellConnectivity_h
+
+#include "CellTypeToType.h"
+#include "ContinousCellInfo.h"
+
+#include <algorithm>
+#include <vector>
+
+namespace smoab { namespace detail {
+
+namespace internal
+{
+ //we want a subset of the real connetivity array,
+ //this does that for use with a super easy wrapper
+ struct SubsetArray
+ {
+ SubsetArray(EntityHandle* realConn,
+ int numCells,
+ int currentVertsPerCell,
+ int newVertsPerCell):
+ Array()
+ {
+ const int size = numCells*newVertsPerCell;
+ this->Array.reserve(size);
+ if(currentVertsPerCell == newVertsPerCell)
+ {
+ std::copy(realConn,realConn+size, std::back_inserter(this->Array));
+ }
+ else
+ {
+ //skip copy only the first N points which we want
+ //since moab stores linear points first per cell
+ EntityHandle *pos = realConn;
+ for(int i=0; i < numCells;++i)
+ {
+ std::copy(pos,pos+newVertsPerCell,std::back_inserter(this->Array));
+ pos += currentVertsPerCell;
+ }
+ }
+ }
+ typedef std::vector<EntityHandle>::const_iterator const_iterator;
+ typedef std::vector<EntityHandle>::iterator iterator;
+
+ const_iterator begin() const { return this->Array.begin(); }
+ iterator begin() { return this->Array.begin(); }
+
+ const_iterator end() const { return this->Array.end(); }
+ iterator end(){ return this->Array.end(); }
+
+ private:
+ std::vector<EntityHandle> Array;
+ };
+}
+
+class LinearCellConnectivity
+{
+public:
+
+ LinearCellConnectivity(smoab::Range const& cells, moab::Interface* moab):
+ Connectivity(),
+ UniquePoints(),
+ Info()
+ {
+ int count = 0;
+ const std::size_t cellSize=cells.size();
+ while(count != cellSize)
+ {
+ EntityHandle* connectivity;
+ int numVerts=0, iterationCount=0;
+ //use the highly efficent calls, since we know that are of the same dimension
+ moab->connect_iterate(cells.begin()+count,
+ cells.end(),
+ connectivity,
+ numVerts,
+ iterationCount);
+ //if we didn't read anything, break!
+ if(iterationCount == 0)
+ {
+ break;
+ }
+
+ //identify the cell type that we currently have,
+ //store that along with the connectivity in a temp storage vector
+ const moab::EntityType type = moab->type_from_handle(*cells.begin()+count);
+
+ int vtkNumVerts;
+ int vtkCellType = smoab::detail::vtkLinearCellType(type,vtkNumVerts);
+
+ ContinousCellInfo info = { vtkCellType, vtkNumVerts, 0, iterationCount };
+ this->Info.push_back(info);
+
+
+ //we need to copy only a subset of the connectivity array
+ internal::SubsetArray conn(connectivity,iterationCount,numVerts,vtkNumVerts);
+ this->Connectivity.push_back(conn);
+
+ count += iterationCount;
+ }
+ }
+
+ //----------------------------------------------------------------------------
+ void compactIds(vtkIdType& numCells, vtkIdType& connectivityLength)
+ {
+ //converts all the ids to be ordered starting at zero, and also
+ //keeping the orginal logical ordering. Stores the result of this
+ //operation in the unstrucutred grid that is passed in
+
+ //lets determine the total length of the connectivity
+ connectivityLength = 0;
+ numCells = 0;
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i)
+ {
+ connectivityLength += (*i).numCells * (*i).numVerts;
+ numCells += (*i).numCells;
+ }
+
+ this->UniquePoints.reserve(connectivityLength);
+
+ this->copyConnectivity(this->UniquePoints);
+ std::sort(this->UniquePoints.begin(),this->UniquePoints.end());
+
+ typedef std::vector<EntityHandle>::iterator EntityIterator;
+ EntityIterator newEnd = std::unique(this->UniquePoints.begin(),
+ this->UniquePoints.end());
+
+ const std::size_t newSize = std::distance(this->UniquePoints.begin(),newEnd);
+ this->UniquePoints.resize(newSize);
+ }
+
+ //----------------------------------------------------------------------------
+ void moabPoints(smoab::Range& range) const
+ {
+ //from the documentation a reverse iterator is the fastest way
+ //to insert into a range.
+ std::copy(this->UniquePoints.rbegin(),
+ this->UniquePoints.rend(),
+ moab::range_inserter(range));
+ }
+
+ //----------------------------------------------------------------------------
+ //copy the connectivity from the moab held arrays to the user input vector
+ void copyConnectivity(std::vector<EntityHandle>& output) const
+ {
+ //walk the info to find the length of each sub connectivity array,
+ //and insert them into the vector, ordering is implied by the order
+ //the connecitivy sub array are added to this class
+ ConnConstIterator c = this->Connectivity.begin();
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i,++c)
+ {
+ //remember our Connectivity is a vector of pointers whose
+ //length is held in the info vector.
+ const int numUnusedPoints = (*i).numUnusedVerts;
+ const int connLength = (*i).numCells * (*i).numVerts;
+ std::copy(c->begin(),c->end(),std::back_inserter(output));
+ }
+ }
+
+ //copy the information from this contianer to a vtk cell array, and
+ //related lookup information
+ void copyToVtkCellInfo(vtkIdType* cellArray) const
+ {
+ ConnConstIterator c = this->Connectivity.begin();
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i, ++c)
+ {
+ //for this group of the same cell type we need to fill the cellTypes
+ const int numCells = (*i).numCells;
+ const int numVerts = (*i).numVerts;
+
+ //for each cell in this collection that have the same type
+ //grab the raw array now, so we can properly increment for each vert in each cell
+ internal::SubsetArray::const_iterator moabConnectivity = c->begin();
+ for(int j=0;j < numCells; ++j)
+ {
+ //cell arrays start and end are different, since we
+ //have to account for element that states the length of each cell
+ cellArray[0]=numVerts;
+
+
+ for(int k=0; k < numVerts; ++k, ++moabConnectivity )
+ {
+ //this is going to be a root of some failures when we start
+ //reading really large datasets under 32bit.
+
+
+ //fyi, don't use a range ds for unique points, distance
+ //function is horribly slow they need to override it
+ EntityConstIterator result = std::lower_bound(
+ this->UniquePoints.begin(),
+ this->UniquePoints.end(),
+ *moabConnectivity);
+ std::size_t newId = std::distance(this->UniquePoints.begin(),
+ result);
+ cellArray[k+1] = static_cast<vtkIdType>(newId);
+ }
+ cellArray += numVerts+1;
+ }
+ }
+ }
+
+private:
+ std::vector<internal::SubsetArray> Connectivity;
+ std::vector<EntityHandle> UniquePoints;
+
+ std::vector<detail::ContinousCellInfo> Info;
+
+ typedef std::vector<EntityHandle>::const_iterator EntityConstIterator;
+ typedef std::vector<internal::SubsetArray>::const_iterator ConnConstIterator;
+ typedef std::vector<detail::ContinousCellInfo>::const_iterator InfoConstIterator;
+};
+} } //namespace smoab::detail
+
+#endif // __smoab_LinearCellConnectivity_h
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/dcd789d7817e/
Changeset: dcd789d7817e
Branch: master
User: vijaysm
Date: 2014-07-13 00:49:31
Summary: Merged in judajake/moab/update_vtk_plugin (pull request #23)
Update VTKReader Paraview plugin for MOAB
Affected #: 19 files
diff --git a/tools/vtkMOABReaderNew/CMakeLists.txt.in b/tools/vtkMOABReaderNew/CMakeLists.txt.in
index 8b8f888..e005289 100644
--- a/tools/vtkMOABReaderNew/CMakeLists.txt.in
+++ b/tools/vtkMOABReaderNew/CMakeLists.txt.in
@@ -13,13 +13,22 @@ find_package(ParaView REQUIRED)
include(${PARAVIEW_USE_FILE})
include_directories(${PARAVIEW_INCLUDE_DIRS})
include_directories(@srcdir@)
+include_directories(@srcdir@/detail)
set(headers
- @srcdir@/SimpleMoab.h
- @srcdir@/CellTypeToType.h
+ @srcdir@/CellSets.h
@srcdir@/DataSetConverter.h
- @srcdir@/MixedCellConnectivity.h
- @srcdir@/vtkMoabReader.h
+ @srcdir@/detail/CellTypeToType.h
+ @srcdir@/detail/ContinousCellInfo.h
+ @srcdir@/detail/LoadGeometry.h
+ @srcdir@/detail/MixedCellConnectivity.h
+ @srcdir@/detail/ReadSparseTag.h
+ @srcdir@/detail/ReduceSpectralMesh.h
+ @srcdir@/detail/ReduceSpectralMesh.h
+ @srcdir@/detail/UsageTable.h
+ @srcdir@/ExtractShell.h
+ @srcdir@/FaceSets.h
+ @srcdir@/SimpleMoab.h
)
add_paraview_plugin(vtkMoabReaderPlugin "5.0"
diff --git a/tools/vtkMOABReaderNew/CellSets.h b/tools/vtkMOABReaderNew/CellSets.h
new file mode 100644
index 0000000..6171aec
--- /dev/null
+++ b/tools/vtkMOABReaderNew/CellSets.h
@@ -0,0 +1,75 @@
+
+#ifndef __smoab_CellSets_h
+#define __smoab_CellSets_h
+
+#include "SimpleMoab.h"
+
+namespace smoab
+{
+//----------------------------------------------------------------------------
+class CellSet
+{
+public:
+ CellSet(smoab::EntityHandle p,const smoab::Range& cells):
+ Entity(p),
+ Cells(cells)
+ {}
+
+ const smoab::Range& cells() const { return this->Cells; }
+ EntityHandle entity() const { return this->Entity; }
+
+ bool contains(smoab::EntityHandle c) const
+ {
+ return this->Cells.find(c) != this->Cells.end();
+ }
+
+ void erase(smoab::Range cells)
+ {
+ //seems that erase() has a bug, so use subtract
+ this->Cells = smoab::subtract(this->Cells,cells);
+ }
+
+private:
+ smoab::EntityHandle Entity;
+ smoab::Range Cells;
+};
+
+//----------------------------------------------------------------------------
+//CellSets are just a vector of CellSets
+typedef std::vector<CellSet> CellSets;
+
+//----------------------------------------------------------------------------
+//templated so it works with FaceCellSets and CellSets
+template<typename T>
+smoab::Range getParents(const T& set)
+{
+ typedef typename T::const_iterator iterator;
+ smoab::Range result;
+
+ for(iterator i=set.begin(); i != set.end(); ++i)
+ {
+ result.insert(i->entity());
+ }
+ return result;
+}
+
+//----------------------------------------------------------------------------
+//templated so it works with FaceCellSets and CellSets
+template<typename T>
+smoab::Range getAllCells(const T& set)
+{
+ typedef typename T::const_iterator iterator;
+ smoab::Range result;
+
+ for(iterator i=set.begin(); i != set.end(); ++i)
+ {
+ smoab::Range c = i->cells();
+ result.insert(c.begin(),c.end());
+ }
+ return result;
+}
+
+
+}
+
+#endif
diff --git a/tools/vtkMOABReaderNew/CellTypeToType.h b/tools/vtkMOABReaderNew/CellTypeToType.h
deleted file mode 100644
index 6ff91af..0000000
--- a/tools/vtkMOABReaderNew/CellTypeToType.h
+++ /dev/null
@@ -1,59 +0,0 @@
-#ifndef __smoab_CellTypeToType_h
-#define __smoab_CellTypeToType_h
-
-#include "SimpleMoab.h"
-#include "vtkCellType.h"
-
-namespace smoab
-{
-template<int> struct CellTypeToType;
-
-int vtkCellType(moab::EntityType t, int &num_connect)
- {
- int ctype = -1;
- switch (t)
- {
- case moab::MBEDGE:
- if (num_connect == 2) ctype = VTK_LINE;
- else if (num_connect == 3) ctype = VTK_QUADRATIC_EDGE;
- break;
- case moab::MBTRI:
- if (num_connect == 3) ctype = VTK_TRIANGLE;
- else if (num_connect == 6) ctype = VTK_QUADRATIC_TRIANGLE;
- else if (num_connect == 7) ctype = VTK_BIQUADRATIC_TRIANGLE;
- break;
- case moab::MBQUAD:
- if (num_connect == 4) ctype = VTK_QUAD;
- else if (num_connect == 8) ctype = VTK_QUADRATIC_QUAD;
- else if (num_connect == 9) ctype = VTK_BIQUADRATIC_QUAD;
- break;
- case moab::MBPOLYGON:
- if (num_connect == 4) ctype = VTK_POLYGON;
- break;
- case moab::MBTET:
- if (num_connect == 4) ctype = VTK_TETRA;
- else if (num_connect == 10) ctype = VTK_QUADRATIC_TETRA;
- break;
- case moab::MBPYRAMID:
- if (num_connect == 5) ctype = VTK_PYRAMID;
- else if (num_connect == 13) ctype = VTK_QUADRATIC_PYRAMID;
- break;
- case moab::MBPRISM:
- if (num_connect == 6) ctype = VTK_WEDGE;
- else if (num_connect == 15) ctype = VTK_QUADRATIC_WEDGE;
- break;
- case moab::MBHEX:
- if (num_connect == 8) ctype = VTK_HEXAHEDRON;
- else if (num_connect == 20) ctype = VTK_QUADRATIC_HEXAHEDRON;
- else if (num_connect == 21) ctype = VTK_QUADRATIC_HEXAHEDRON, num_connect = 20;
- else if (num_connect == 27) ctype = VTK_TRIQUADRATIC_HEXAHEDRON;
- break;
- default:
- ctype = -1;
- break;
- }
- return ctype;
- }
-}
-
-#endif // CELLTYPETOTYPE_H
diff --git a/tools/vtkMOABReaderNew/DataSetConverter.h b/tools/vtkMOABReaderNew/DataSetConverter.h
index bfae086..3edadf0 100644
--- a/tools/vtkMOABReaderNew/DataSetConverter.h
+++ b/tools/vtkMOABReaderNew/DataSetConverter.h
@@ -2,20 +2,16 @@
#define __smoab_DataSetConverter_h
#include "SimpleMoab.h"
-#include "CellTypeToType.h"
-#include "MixedCellConnectivity.h"
+#include "CellSets.h"
+#include "detail/LoadGeometry.h"
+#include "detail/ReadSparseTag.h"
-#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkDoubleArray.h>
#include <vtkFieldData.h>
#include <vtkIntArray.h>
#include <vtkIdTypeArray.h>
#include <vtkNew.h>
-#include <vtkPointData.h>
-#include <vtkPoints.h>
-#include <vtkUnsignedCharArray.h>
-#include <vtkUnstructuredGrid.h>
#include <algorithm>
@@ -35,17 +31,13 @@ public:
Moab(interface.Moab),
Tag(tag),
ReadMaterialIds(false),
- ReadProperties(false),
- MaterialName("Material")
+ ReadProperties(false)
{
}
void readMaterialIds(bool add) { this->ReadMaterialIds = add; }
bool readMaterialIds() const { return this->ReadMaterialIds; }
- void materialIdName(const std::string& name) { this->MaterialName = name; }
- const std::string& materialIdName() const { return this->MaterialName; }
-
void readProperties(bool readProps) { this->ReadProperties = readProps; }
bool readProperties() const { return this->ReadProperties; }
@@ -54,109 +46,65 @@ public:
//grid. Currently doesn't support reading properties.
//Will read in material ids, if no material id is assigned to an entity,
//its cells will be given an unique id
+ template<typename VTKGridType>
bool fill(const smoab::Range& entities,
- vtkUnstructuredGrid* grid) const
+ VTKGridType* grid) const
{
//create a helper datastructure which can determines all the unique point ids
//and converts moab connecitvity info to vtk connectivity
- moab::Range cells;
- //append all the entities cells together into a single range
+
+ //get all the cells for each parent entity and create
+ // an entity set of those items
+ int dim = this->Tag->value();
typedef smoab::Range::const_iterator iterator;
+ smoab::CellSets entitySets;
for(iterator i=entities.begin(); i!= entities.end(); ++i)
{
+ smoab::Range entitiesCells;
if(this->Tag->isComparable())
{
//if we are comparable only find the cells that match our tags dimension
- smoab::Range entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());
- cells.insert(entitiesCells.begin(),entitiesCells.end());
+ entitiesCells = this->Interface.findEntitiesWithDimension(*i,dim,true);
}
else
{
- //this is a bad representation of all other tags, but we are presuming that
- //neuman and dirichlet are on entitysets with no children
- this->Moab->get_entities_by_handle(*i,cells);
+ entitiesCells = this->Interface.findHighestDimensionEntities(*i,true);
}
+ smoab::CellSet set(*i,entitiesCells);
+ entitySets.push_back(set);
}
- smoab::Range points;
- this->loadCellsAndPoints(cells,points,grid);
-
- if(this->readMaterialIds())
- {
- typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
- typedef std::vector<int>::const_iterator IdConstIterator;
- typedef std::vector<int>::iterator IdIterator;
-
- std::vector<smoab::EntityHandle> searchableCells;
- searchableCells.reserve(grid->GetNumberOfCells());
- std::copy(cells.begin(),cells.end(),std::back_inserter(searchableCells));
- cells.clear(); //release memory we don't need
+ moab::Range cells = smoab::getAllCells(entitySets);
- std::vector<int> materialIds(entities.size());
- //first off iterate the entities and determine which ones
- //have moab material ids
+ //convert the datastructure from a list of cells to a vtk data set
+ detail::LoadGeometry loadGeom(cells,dim,this->Interface);
+ loadGeom.fill(grid);
- //wrap this area with scope, to remove local variables
+ if(this->readMaterialIds())
{
- smoab::MaterialTag tag;
- IdIterator materialIndex = materialIds.begin();
- for(iterator i=entities.begin();
- i != entities.end();
- ++i, ++materialIndex)
- {
- moab::Tag mtag = this->Interface.getMoabTag(tag);
-
- int value=-1;
- this->Moab->tag_get_data(mtag,&(*i),1,&value);
- *materialIndex=static_cast<int>(value);
- }
-
- //now determine ids for all entities that don't have materials
- IdConstIterator maxPos = std::max_element(materialIds.begin(),
- materialIds.end());
- int maxMaterial = *maxPos;
- for(IdIterator i=materialIds.begin(); i!= materialIds.end(); ++i)
- {
- if(*i==-1)
- {
- *i = ++maxMaterial;
- }
- }
- }
- //now we create the material field, and set all the values
- vtkNew<vtkIntArray> materialSet;
- materialSet->SetName(this->materialIdName().c_str());
- materialSet->SetNumberOfValues(grid->GetNumberOfCells());
+ detail::ReadSparseTag materialTagReading(entitySets,
+ cells,
+ this->Interface);
- IdConstIterator materialValue = materialIds.begin();
- for(iterator i=entities.begin(); i!= entities.end(); ++i, ++materialValue)
- {
- //this is a time vs memory trade off, I don't want to store
- //the all the cell ids twice over, lets use more time
- smoab::Range entitiesCells;
- if(this->Tag->isComparable())
- {entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());}
- else
- {this->Moab->get_entities_by_handle(*i,entitiesCells);}
-
- EntityHandleIterator s_begin = searchableCells.begin();
- EntityHandleIterator s_end = searchableCells.end();
- for(iterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
- {
- EntityHandleIterator result = std::lower_bound(s_begin,
- s_end,
- *j);
- std::size_t newId = std::distance(s_begin,result);
- materialSet->SetValue(static_cast<int>(newId), *materialValue);
- }
- }
+ smoab::MaterialTag mtag;
+ vtkNew<vtkIntArray> materials;
+ materials->SetName(mtag.name());
+ materialTagReading.fill(materials.GetPointer(),&mtag);
+ grid->GetCellData()->AddArray(materials.GetPointer());
+ }
- grid->GetCellData()->AddArray(materialSet.GetPointer());
+ //by default we always try to load the default tag
+ detail::ReadSparseTag sTagReading(entitySets,
+ cells,
+ this->Interface);
- }
+ vtkNew<vtkIntArray> sparseTagData;
+ sparseTagData->SetName(this->Tag->name());
+ sTagReading.fill(sparseTagData.GetPointer(),this->Tag);
+ grid->GetCellData()->AddArray(sparseTagData.GetPointer());
return true;
}
@@ -165,28 +113,32 @@ public:
//given a single entity handle create a unstructured grid from it.
//optional third parameter is the material id to use if readMaterialIds
//is on, and no material sparse tag is found for this entity
+ template<typename VTKGridType>
bool fill(const smoab::EntityHandle& entity,
- vtkUnstructuredGrid* grid,
+ VTKGridType* grid,
const int materialId=0) const
{
//create a helper datastructure which can determines all the unique point ids
//and converts moab connecitvity info to vtk connectivity
- moab::Range cells;
+
+ smoab::Range cells;
+ int dim = this->Tag->value();
if(this->Tag->isComparable())
{
//if we are comparable only find the cells that match our tags dimension
- cells = this->Interface.findEntitiesWithDimension(entity,Tag->value());
+ cells = this->Interface.findEntitiesWithDimension(entity,dim,true);
}
else
{
- //this is a bad representation of all other tags, but we are presuming that
- //neuman and dirichlet are on entitysets with no children
- this->Moab->get_entities_by_handle(entity,cells);
+ //load subentities
+ cells = this->Interface.findHighestDimensionEntities(entity,true);
}
- smoab::Range points;
- this->loadCellsAndPoints(cells,points,grid);
+ //convert the datastructure from a list of cells to a vtk data set
+ detail::LoadGeometry loadGeom(cells,dim,this->Interface);
+ loadGeom.fill(grid);
+ const smoab::Range& points = loadGeom.moabPoints();
if(this->readProperties())
{
@@ -194,53 +146,34 @@ public:
this->readProperties(points,grid->GetPointData());
}
+ smoab::CellSets cellSets;
+ smoab::CellSet set(entity,cells);
+ cellSets.push_back(set);
if(this->readMaterialIds())
{
- this->readSparseTag(smoab::MaterialTag(),entity,
- grid->GetNumberOfCells(),
- grid->GetCellData(),
- materialId);
- }
- return true;
- }
+ smoab::MaterialTag mtag;
+ detail::ReadSparseTag materialTagReading(cellSets,
+ cells,
+ this->Interface);
- //----------------------------------------------------------------------------
- void loadCellsAndPoints(const smoab::Range& cells,
- smoab::Range& points,
- vtkUnstructuredGrid* grid) const
- {
+ vtkNew<vtkIntArray> materials;
+ materials->SetName(mtag.name());
+ materialTagReading.fill(materials.GetPointer(),&mtag);
+ grid->GetCellData()->AddArray(materials.GetPointer());
- smoab::MixedCellConnectivity mixConn(cells,this->Moab);
-
- //now that mixConn has all the cells properly stored, lets fixup
- //the ids so that they start at zero and keep the same logical ordering
- //as before.
- vtkIdType numCells, connLen;
- mixConn.compactIds(numCells,connLen);
- this->setGridsTopology(mixConn,grid,numCells,connLen);
+ }
- mixConn.moabPoints(points);
+ //by default we always try to load the default tag
+ detail::ReadSparseTag sTagReading(cellSets,
+ cells,
+ this->Interface);
- vtkNew<vtkPoints> newPoints;
- this->addCoordinates(points,newPoints.GetPointer());
- grid->SetPoints(newPoints.GetPointer());
+ vtkNew<vtkIntArray> sparseTagData;
+ sparseTagData->SetName(this->Tag->name());
+ sTagReading.fill(sparseTagData.GetPointer(),this->Tag);
+ grid->GetCellData()->AddArray(sparseTagData.GetPointer());
- }
-
- //----------------------------------------------------------------------------
- void addCoordinates(smoab::Range pointEntities, vtkPoints* pointContainer) const
- {
- //since the smoab::range are always unique and sorted
- //we can use the more efficient coords_iterate
- //call in moab, which returns moab internal allocated memory
- pointContainer->SetDataTypeToDouble();
- pointContainer->SetNumberOfPoints(pointEntities.size());
-
- //need a pointer to the allocated vtkPoints memory so that we
- //don't need to use an extra copy and we can bypass all vtk's check
- //on out of bounds
- double *rawPoints = static_cast<double*>(pointContainer->GetVoidPointer(0));
- this->Moab->get_coords(pointEntities,rawPoints);
+ return true;
}
private:
@@ -260,36 +193,6 @@ private:
}
//----------------------------------------------------------------------------
- bool readSparseTag(smoab::Tag tag,
- smoab::EntityHandle const& entity,
- vtkIdType length,
- vtkFieldData* field,
- vtkIdType defaultValue) const
- {
-
- typedef std::vector<moab::Tag>::const_iterator iterator;
- moab::Tag mtag = this->Interface.getMoabTag(tag);
-
- int value=0;
- moab::ErrorCode rval = this->Moab->tag_get_data(mtag,&entity,1,&value);
- if(rval!=moab::MB_SUCCESS)
- {
- value = defaultValue;
- }
-
- vtkNew<vtkIntArray> materialSet;
- materialSet->SetNumberOfValues(length);
- materialSet->SetName(this->materialIdName().c_str());
-
- int *raw = static_cast<int*>(materialSet->GetVoidPointer(0));
- std::fill(raw,raw+length,value);
-
- field->AddArray(materialSet.GetPointer());
-
- return true;
- }
-
- //----------------------------------------------------------------------------
void readDenseTags(std::vector<moab::Tag> &tags,
smoab::Range const& entities,
vtkFieldData* field) const
@@ -359,36 +262,7 @@ private:
}
}
- //----------------------------------------------------------------------------
- void setGridsTopology(smoab::MixedCellConnectivity const& mixedCells,
- vtkUnstructuredGrid* grid,
- vtkIdType numCells,
- vtkIdType numConnectivity) const
- {
- //correct the connectivity size to account for the vtk padding
- const vtkIdType vtkConnectivity = numCells + numConnectivity;
- vtkNew<vtkIdTypeArray> cellArray;
- vtkNew<vtkIdTypeArray> cellLocations;
- vtkNew<vtkUnsignedCharArray> cellTypes;
-
- cellArray->SetNumberOfValues(vtkConnectivity);
- cellLocations->SetNumberOfValues(numCells);
- cellTypes->SetNumberOfValues(numCells);
-
- vtkIdType* rawArray = static_cast<vtkIdType*>(cellArray->GetVoidPointer(0));
- vtkIdType* rawLocations = static_cast<vtkIdType*>(cellLocations->GetVoidPointer(0));
- unsigned char* rawTypes = static_cast<unsigned char*>(cellTypes->GetVoidPointer(0));
-
- mixedCells.copyToVtkCellInfo(rawArray,rawLocations,rawTypes);
-
- vtkNew<vtkCellArray> cells;
- cells->SetCells(numCells,cellArray.GetPointer());
- grid->SetCells(cellTypes.GetPointer(),
- cellLocations.GetPointer(),
- cells.GetPointer(),
- NULL,NULL);
- }
};
}
diff --git a/tools/vtkMOABReaderNew/ExtractShell.h b/tools/vtkMOABReaderNew/ExtractShell.h
new file mode 100644
index 0000000..359fbce
--- /dev/null
+++ b/tools/vtkMOABReaderNew/ExtractShell.h
@@ -0,0 +1,80 @@
+#ifndef __smoab_ExtractShell_h
+#define __smoab_ExtractShell_h
+
+#include "SimpleMoab.h"
+#include "detail/UsageTable.h"
+
+#include <algorithm>
+
+namespace smoab{
+
+class ExtractShell
+{
+ const smoab::Interface& Interface;
+ smoab::CellSets VCells;
+
+public:
+ ExtractShell(const smoab::CellSets volCells,
+ const smoab::Interface& interface):
+ Interface(interface),
+ VCells(volCells)
+ {
+ }
+
+ bool findSkins(smoab::CellSets &surfaceCellSets);
+};
+
+
+//----------------------------------------------------------------------------
+bool ExtractShell::findSkins(smoab::CellSets &surfaceCellSets)
+{
+ typedef smoab::Range::const_iterator Iterator;
+
+ typedef smoab::CellSets::const_iterator SetIterator;
+
+
+ smoab::Range cellsToRemove;
+ for(SetIterator set = this->VCells.begin();
+ set != this->VCells.end();
+ ++set)
+ {
+ const smoab::Range &cells = set->cells();
+ this->Interface.createAdjacencies(set->cells(),2);
+
+
+ //we create the usage table for each iteration so that we only
+ //get the shell of each cell set. If we used the table between
+ //sets we would get the shell of the combined sets
+ smoab::detail::UsageTable table;
+ for(Iterator i = cells.begin(); i != cells.end(); ++i)
+ {
+ std::vector<smoab::EntityHandle> faceCells =
+ this->Interface.sideElements(*i,2);
+
+ //the usage id allows you to label cells when going into the table
+ //so that you can extract multiple shells where each is based on
+ //a single region id.
+ std::vector<int> regionId(1,faceCells.size());
+ table.incrementUsage(faceCells,regionId);
+ }
+ smoab::Range surfaceCells = table.singleUsage();
+
+ //create a new cell set that
+ smoab::CellSet surfaceSet(set->entity(),surfaceCells);
+ surfaceCellSets.push_back(surfaceSet);
+
+ smoab::Range subsetToRemove = table.multipleUsage();
+ cellsToRemove.insert(subsetToRemove.begin(),subsetToRemove.end());
+ }
+
+ //we will remove all cells that have multiple usages from the moab database
+ //I really don't care if they already existed or not.
+ this->Interface.remove(cellsToRemove);
+ return true;
+}
+
+
+
+}
+
+#endif // __smoab_ExtractShell_h
diff --git a/tools/vtkMOABReaderNew/FaceSets.h b/tools/vtkMOABReaderNew/FaceSets.h
new file mode 100644
index 0000000..536f3f9
--- /dev/null
+++ b/tools/vtkMOABReaderNew/FaceSets.h
@@ -0,0 +1,256 @@
+#ifndef __smoab_FaceSets_h
+#define __smoab_FaceSets_h
+
+#include "CellSets.h"
+#include <set>
+
+namespace smoab
+{
+//----------------------------------------------------------------------------
+class FaceCellSet : public CellSet
+{
+public:
+ FaceCellSet(int id, smoab::EntityHandle p,const smoab::Range& cells):
+ CellSet(p,cells),
+ ID(id)
+ {}
+
+ int faceId() const { return ID; }
+ void overrideFaceId(int i) { ID = i; } //USE AT YOUR OWN RISK
+
+private:
+ int ID;
+};
+
+
+//----------------------------------------------------------------------------
+typedef std::vector<FaceCellSet> FaceCellSets;
+
+
+//----------------------------------------------------------------------------
+//class that store the regions that a faces is adjacent too
+struct FacesAdjRegions
+{
+ FacesAdjRegions(int f, smoab::EntityHandle r0, smoab::EntityHandle r1):
+ FaceId(f),
+ Region0(r0),
+ Region1(r1)
+ {
+ if (r0 > r1)
+ {
+ std::swap(this->Region0,this->Region1);
+ }
+ }
+
+ FacesAdjRegions(int f):
+ FaceId(f),
+ Region0(-3),
+ Region1(-2)
+ {}
+
+ bool operator<(const FacesAdjRegions& other) const
+ {
+ return (this->FaceId < other.FaceId);
+ }
+
+ smoab::EntityHandle otherId(smoab::EntityHandle other) const
+ {
+ if(other == Region0)
+ {
+ return Region1;
+ }
+ return Region0;
+ }
+
+ int FaceId;
+ smoab::EntityHandle Region0;
+ smoab::EntityHandle Region1;
+};
+//----------------------------------------------------------------------------
+smoab::FaceCellSets findFaceSets(smoab::CellSets shells,
+ smoab::CellSets boundaries,
+ std::set<smoab::FacesAdjRegions>& faceMaps)
+{
+ typedef smoab::CellSets::iterator iterator;
+ typedef smoab::FaceCellSets::iterator faceIterator;
+ typedef std::set<smoab::FacesAdjRegions>::const_iterator FaceAdjIterator;
+
+ //we need to properly label each unique face in shells
+ //we do this by intersecting each shell with each other shell
+ //to find shell on shell contact, and than we intersect each
+ //resulting shell with the boundary conditions
+ //the end result of these intersections will be the new modelfaces
+ int faceId = 1;
+ smoab::FaceCellSets shellFaces;
+
+ //first intersect each shell with each other shell
+ std::set<smoab::FacesAdjRegions> shellFaceContacts;
+ for(iterator i=shells.begin();i!= shells.end();++i)
+ {
+ //copy the cells so we can add a face that represents
+ //all the cells of the region that aren't shared with another region
+ int numCells = i->cells().size(); //size() on range is slow, so cache it
+ for(iterator j = i+1;
+ j != shells.end() && numCells > 0;
+ ++j)
+ {
+ //intersect i and j to make a new face
+ smoab::Range intersection = smoab::intersect(i->cells(),j->cells());
+ if(!intersection.empty())
+ {
+ //don't want to increment faceId when the intersection is empty
+ smoab::FaceCellSet face(faceId++,i->entity(),intersection);
+ shellFaces.push_back(face);
+ i->erase(intersection);
+ j->erase(intersection);
+ numCells -= intersection.size();
+
+ //add this to the face map
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),j->entity());
+ shellFaceContacts.insert(faceInfo);
+ }
+ }
+ //if all the cells for shell i are used, don't add a new
+ //empty face
+ if(numCells > 0)
+ {
+ smoab::FaceCellSet face(faceId++,i->entity(),i->cells());
+ shellFaces.push_back(face);
+
+ //add this to the face map
+ smoab::FacesAdjRegions faceInfo(faceId-1,-1,i->entity());
+ shellFaceContacts.insert(faceInfo);
+ }
+ }
+
+ //now we have all the faces that match shell on shell contact
+ //we know process all the new faces to see if they intersect
+ //with any boundary sets. A boundary set can span multiple
+ //shells so we want to process it as a second loop
+
+ //store the end before we start adding boundary faces, which
+ //we don't need to check agianst other boundaries
+ faceId = 1; //reset the faced id
+
+ //store in a new face set, expanding the current one causes incorrect results
+ smoab::FaceCellSets faces;
+ for(faceIterator i=shellFaces.begin();i != shellFaces.end();++i)
+ {
+ //determine from the shell faces if the new face we are creating
+ //is bounded by two regions or just one
+ smoab::FacesAdjRegions idToSearchFor(i->faceId());
+ FaceAdjIterator adjRegions = shellFaceContacts.find(idToSearchFor);
+ smoab::EntityHandle otherRegionId = adjRegions->otherId(i->entity());
+
+ int numCells = i->cells().size(); //size() on range is slow, so cache it
+ for(iterator j=boundaries.begin();j != boundaries.end(); ++j)
+ {
+ smoab::Range intersect = smoab::intersect(i->cells(),j->cells());
+ if(!intersect.empty())
+ {
+ //don't want to increment faceId when the intersection is empty
+ smoab::FaceCellSet face(faceId++,j->entity(),intersect);
+ faces.push_back(face);
+ i->erase(intersect);
+ numCells -= intersect.size();
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),otherRegionId);
+ faceMaps.insert(faceInfo);
+ }
+ }
+ if(numCells > 0)
+ {
+ smoab::FaceCellSet face(faceId++,i->entity(),i->cells());
+ faces.push_back(face);
+ smoab::FacesAdjRegions faceInfo(faceId-1,i->entity(),otherRegionId);
+ faceMaps.insert(faceInfo);
+ }
+ }
+ return faces;
+}
+
+//----------------------------------------------------------------------------
+template<typename T>
+std::vector<T> faceIdsPerCell(const smoab::FaceCellSets& faces)
+{
+ typedef smoab::FaceCellSets::const_iterator iterator;
+ typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
+ typedef smoab::Range::const_iterator RangeIterator;
+
+ //find all the cells that are in the faceCellSet, and than map
+ //the proper face id to that relative position, here comes lower_bounds!
+ std::vector<smoab::EntityHandle> searchableCells;
+ smoab::Range faceRange = smoab::getAllCells(faces);
+ smoab::RangeToVector(faceRange,searchableCells);
+
+ //faceIds will be the resulting array
+ std::vector<T> faceIds(searchableCells.size());
+
+ //construct the start and end iterators for the lower bounds call
+
+ EntityHandleIterator s_begin = searchableCells.begin();
+ EntityHandleIterator s_end = searchableCells.end();
+
+ //search the face cell sets
+ for(iterator i=faces.begin(); i!=faces.end(); ++i)
+ {
+ T value = static_cast<T>(i->faceId());
+ const smoab::Range& entitiesCells = i->cells();
+ for(RangeIterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
+ {
+ EntityHandleIterator result = std::lower_bound(s_begin,
+ s_end,
+ *j);
+ std::size_t newId = std::distance(s_begin,result);
+ faceIds[newId] = value;
+ }
+ }
+ return faceIds;
+}
+
+//----------------------------------------------------------------------------
+//given a face adjacency, determine the regions spare tag values
+template<typename T>
+std::pair<T,T> FaceAdjRegionValues(const smoab::FacesAdjRegions& faceAdj,
+ smoab::Tag* t,
+ const smoab::Interface& interface)
+ {
+ std::pair<T,T> returnValue;
+ const int defaultValue = interface.getDefaultTagVaue<int>(*t);
+
+ /*
+ * IF A REGION IS SET TO -1 WE NEED TO PUSH THAT VALUE DOWN
+ * AS THE MATERIAL, SINCE THE MOAB DEFAULT TAG VALUE WILL
+ * BE CONSIDIERED A REGION, AND WE WANT TO SAY IT BOUNDS THE
+ * VOID REGION
+ */
+ int tagValue = defaultValue; //use tagValue to pass in default value
+ if(faceAdj.Region0 != -1)
+ {
+ tagValue = interface.getTagData(*t,faceAdj.Region0,tagValue);
+ }
+ else
+ {
+ tagValue = -1;
+ }
+
+ //set the first region tag value into the pair we are returing
+ returnValue.first = static_cast<T>(tagValue);
+
+ tagValue = defaultValue; //use tagValue to pass in default value
+ tagValue = interface.getTagData(*t,faceAdj.Region1,tagValue);
+ if(faceAdj.Region1 != -1)
+ {
+ tagValue = interface.getTagData(*t,faceAdj.Region1,tagValue);
+ }
+ else
+ {
+ tagValue = -1;
+ }
+ returnValue.second = static_cast<T>(tagValue);
+
+ return returnValue;
+ }
+
+ } //smoab
+
+#endif
diff --git a/tools/vtkMOABReaderNew/MixedCellConnectivity.h b/tools/vtkMOABReaderNew/MixedCellConnectivity.h
deleted file mode 100644
index 1b26d0f..0000000
--- a/tools/vtkMOABReaderNew/MixedCellConnectivity.h
+++ /dev/null
@@ -1,271 +0,0 @@
-#ifndef __smoab_MixedCellConnectivity_h
-#define __smoab_MixedCellConnectivity_h
-
-#include "vtkCellType.h"
-#include <algorithm>
-
-namespace
-{
-
-template<int N> struct QuadratricOrdering{};
-
-template<> struct QuadratricOrdering<VTK_QUADRATIC_WEDGE>
-{
- static const int NUM_VERTS = 15;
- void reorder(vtkIdType* connectivity) const
- {
- std::swap_ranges(connectivity+9,connectivity+12,connectivity+12);
- }
-};
-
-template<> struct QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON>
-{
- static const int NUM_VERTS = 27;
- void reorder(vtkIdType* connectivity) const
- {
- std::swap_ranges(connectivity+12,connectivity+16,connectivity+16);
-
- //move 20 to 22
- //move 22 to 23
- //move 23 to 20
-
- //swap 20 with 22
- std::swap(connectivity[20],connectivity[23]);
-
- //swap 22 with 23
- std::swap(connectivity[22],connectivity[23]);
- }
-};
-
-template<typename QuadraticOrdering>
-void FixQuadraticIdOrdering(vtkIdType* connectivity, vtkIdType numCells,
- QuadraticOrdering& ordering)
-{
- //skip the first index that holds the length of the cells
- //if we skip it once here, and than properly increment it makes the code
- //far easier
- connectivity+=1;
- for(vtkIdType i=0; i < numCells; ++i)
- {
- ordering.reorder(connectivity);
- connectivity += ordering.NUM_VERTS + 1;
- }
-}
-}
-
-namespace smoab
-{
-
-class MixedCellConnectivity
-{
-public:
- MixedCellConnectivity(smoab::Range const& cells, moab::Interface* moab):
- Connectivity(),
- UniquePoints(),
- Info()
- {
- int count = 0;
- const std::size_t cellSize=cells.size();
- while(count != cellSize)
- {
- EntityHandle* connectivity;
- int numVerts=0, iterationCount=0;
- //use the highly efficent calls, since we know that are of the same dimension
- moab->connect_iterate(cells.begin()+count,
- cells.end(),
- connectivity,
- numVerts,
- iterationCount);
- //if we didn't read anything, break!
- if(iterationCount == 0)
- {
- break;
- }
-
- //identify the cell type that we currently have,
- //store that along with the connectivity in a temp storage vector
- const moab::EntityType type = moab->type_from_handle(*cells.begin()+count);
-
- //while all these cells are contiously of the same type,
- //quadric hexs in vtk have 20 points, but moab has 21 so we
- //need to store this difference
- int numVTKVerts = numVerts;
- int vtkCellType = smoab::vtkCellType(type,numVTKVerts);
-
- RunLengthInfo info = { vtkCellType, numVerts, (numVerts-numVTKVerts), iterationCount };
- this->Info.push_back(info);
- this->Connectivity.push_back(connectivity);
-
- count += iterationCount;
- }
- }
-
- //----------------------------------------------------------------------------
- void compactIds(vtkIdType& numCells, vtkIdType& connectivityLength)
- {
- //converts all the ids to be ordered starting at zero, and also
- //keeping the orginal logical ordering. Stores the result of this
- //operation in the unstrucutred grid that is passed in
-
- //lets determine the total length of the connectivity
- connectivityLength = 0;
- numCells = 0;
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i)
- {
- connectivityLength += (*i).numCells * (*i).numVerts;
- numCells += (*i).numCells;
- }
-
- this->UniquePoints.reserve(connectivityLength);
-
- this->copyConnectivity(this->UniquePoints);
- std::sort(this->UniquePoints.begin(),this->UniquePoints.end());
-
- typedef std::vector<EntityHandle>::iterator EntityIterator;
- EntityIterator newEnd = std::unique(this->UniquePoints.begin(),
- this->UniquePoints.end());
-
- const std::size_t newSize = std::distance(this->UniquePoints.begin(),newEnd);
- this->UniquePoints.resize(newSize);
- }
-
- //----------------------------------------------------------------------------
- void moabPoints(smoab::Range& range) const
- {
- //from the documentation a reverse iterator is the fastest way
- //to insert into a range.
- std::copy(this->UniquePoints.rbegin(),
- this->UniquePoints.rend(),
- moab::range_inserter(range));
- }
-
- //----------------------------------------------------------------------------
- //copy the connectivity from the moab held arrays to the user input vector
- void copyConnectivity(std::vector<EntityHandle>& output) const
- {
- //walk the info to find the length of each sub connectivity array,
- //and insert them into the vector, ordering is implied by the order
- //the connecitivy sub array are added to this class
- ConnConstIterator c = this->Connectivity.begin();
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i,++c)
- {
- //remember our Connectivity is a vector of pointers whose
- //length is held in the info vector.
- const int numUnusedPoints = (*i).numUnusedVerts;
- if(numUnusedPoints==0)
- {
- const int connLength = (*i).numCells * (*i).numVerts;
- std::copy(*c,*c+connLength,std::back_inserter(output));
- }
- else
- {
- //we have cell connectivity that we need to skip,
- //so we have to manual copy each cells connectivity
- const int size = (*i).numCells;
- const int numPoints = (*i).numVerts;
- for(int i=0; i < size; ++i)
- {
- std::copy(*c,*c+numPoints,std::back_inserter(output));
- }
- c+=numPoints + (*i).numUnusedVerts;
- }
-
- }
- }
-
- //copy the information from this contianer to a vtk cell array, and
- //related lookup information
- void copyToVtkCellInfo(vtkIdType* cellArray,
- vtkIdType* cellLocations,
- unsigned char* cellTypes) const
- {
- vtkIdType currentVtkConnectivityIndex = 0;
- ConnConstIterator c = this->Connectivity.begin();
- for(InfoConstIterator i = this->Info.begin();
- i != this->Info.end();
- ++i, ++c)
- {
- //for this group of the same cell type we need to fill the cellTypes
- const int numCells = (*i).numCells;
- const int numVerts = (*i).numVerts;
-
- std::fill_n(cellTypes,
- numCells,
- static_cast<unsigned char>((*i).type));
-
- //for each cell in this collection that have the same type
- //grab the raw array now, so we can properly increment for each vert in each cell
- EntityHandle* moabConnectivity = *c;
- for(int j=0;j < numCells; ++j)
- {
- cellLocations[j]= currentVtkConnectivityIndex;
-
- //cell arrays start and end are different, since we
- //have to account for element that states the length of each cell
- cellArray[0]=numVerts;
-
- for(int k=0; k < numVerts; ++k, ++moabConnectivity )
- {
- //this is going to be a root of some failures when we start
- //reading really large datasets under 32bit.
-
-
- //fyi, don't use a range ds for unique points, distance
- //function is horribly slow they need to override it
- EntityConstIterator result = std::lower_bound(
- this->UniquePoints.begin(),
- this->UniquePoints.end(),
- *moabConnectivity);
- std::size_t newId = std::distance(this->UniquePoints.begin(),
- result);
- cellArray[k+1] = static_cast<vtkIdType>(newId);
- }
-
- //skip any extra unused points, which is currnetly only
- //the extra center point in moab quadratic hex
- moabConnectivity+=(*i).numUnusedVerts;
-
- currentVtkConnectivityIndex += numVerts+1;
- cellArray += numVerts+1;
- }
-
- //For Tri-Quadratic-Hex and Quadratric-Wedge Moab and VTK
- //Differ on the order of the edge ids. For wedge we need to swap
- //indices 9,10,11 with 12,13,14 for each cell. For Hex we sawp
- //12,13,14,15 with 16,17,18,19
- int vtkCellType = (*i).type;
- vtkIdType* connectivity = cellArray - (numCells * (numVerts+1));
- if(vtkCellType == VTK_TRIQUADRATIC_HEXAHEDRON)
- {
- ::QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON> newOrdering;
- ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
- }
- else if(vtkCellType == VTK_QUADRATIC_WEDGE)
- {
- ::QuadratricOrdering<VTK_QUADRATIC_WEDGE> newOrdering;
- ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
- }
-
- cellLocations += numCells;
- cellTypes += numCells;
- }
-
- }
-
-private:
- std::vector<EntityHandle*> Connectivity;
- std::vector<EntityHandle> UniquePoints;
-
- struct RunLengthInfo{ int type; int numVerts; int numUnusedVerts; int numCells; };
- std::vector<RunLengthInfo> Info;
-
- typedef std::vector<EntityHandle>::const_iterator EntityConstIterator;
- typedef std::vector<EntityHandle*>::const_iterator ConnConstIterator;
- typedef std::vector<RunLengthInfo>::const_iterator InfoConstIterator;
-};
-}
-#endif // __smoab_MixedCellConnectivity_h
diff --git a/tools/vtkMOABReaderNew/SimpleMoab.h b/tools/vtkMOABReaderNew/SimpleMoab.h
index 8fc54e6..4fb734c 100644
--- a/tools/vtkMOABReaderNew/SimpleMoab.h
+++ b/tools/vtkMOABReaderNew/SimpleMoab.h
@@ -5,6 +5,8 @@
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
+#include "moab/CN.hpp"
+
#include "MBTagConventions.hpp"
#include <iostream>
@@ -26,12 +28,6 @@ using moab::intersect;
using moab::subtract;
using moab::unite;
-//forward declare this->Moab for Tag
-struct Interface;
-
-//forward declar the DataSetConverter so it can be a friend of Interface
-class DataSetConverter;
-
class Tag
{
const std::string Name_;
@@ -66,10 +62,24 @@ public:
GeomTag(int d):Tag("GEOM_DIMENSION"),dim(d){}
GeomTag():Tag("GEOM_DIMENSION"), dim(0){}
+ virtual ~GeomTag(){}
+
bool isComparable() const { return dim > 0; }
int value() const { return dim; }
};
+
+//forward declare this->Moab for Tag
+struct Interface;
+
+//forward declare the DataSetConverter so it can be a friend of Interface
+class DataSetConverter;
+
+//forward declare the LoadGeometry so it can be a friend of Interface
+namespace detail{ class LoadGeometry; }
+namespace detail{ class LoadPoly; }
+
+
//light weight wrapper on a moab this->Moab that exposes only the reduced class
//that we need
class Interface
@@ -102,6 +112,37 @@ public:
}
//----------------------------------------------------------------------------
+ template<typename T>
+ T getDefaultTagVaue(moab::Tag tag) const
+ {
+ T defaultValue;
+ this->Moab->tag_get_default_value(tag,&defaultValue);
+ return defaultValue;
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getDefaultTagVaue(smoab::Tag tag) const
+ {
+ return this->getDefaultTagVaue<T>(getMoabTag(tag));
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getTagData(moab::Tag tag, const smoab::EntityHandle& entity, T value) const
+ {
+ this->Moab->tag_get_data(tag,&entity,1,&value);
+ return value;
+ }
+
+ //----------------------------------------------------------------------------
+ template<typename T>
+ T getTagData(smoab::Tag tag, const smoab::EntityHandle& entity, T value = T()) const
+ {
+ return this->getTagData(getMoabTag(tag),entity,value);
+ }
+
+ //----------------------------------------------------------------------------
//returns the moab name for the given entity handle if it has a sparse Name tag
std::string name(const smoab::EntityHandle& entity) const
{
@@ -119,6 +160,19 @@ public:
return std::string(name);
}
+ //----------------------------------------------------------------------------
+ //returns the geometeric dimension of an entity.
+ int dimension(const smoab::EntityHandle& entity) const
+ {
+ return this->Moab->dimension_from_handle(entity);
+ }
+
+ //----------------------------------------------------------------------------
+ //returns the geometeric dimension of an entity.
+ smoab::EntityType entityType(const smoab::EntityHandle& entity) const
+ {
+ return this->Moab->type_from_handle(entity);
+ }
//----------------------------------------------------------------------------
smoab::EntityHandle getRoot() const { return this->Moab->get_root_set(); }
@@ -133,6 +187,17 @@ public:
}
//----------------------------------------------------------------------------
+ //given a single entity handle find all items in that mesh set that aren't
+ //them selves entitysets. If recurse is true we also recurse sub entitysets
+ smoab::Range findAllMeshEntities(smoab::EntityHandle const& entity,
+ bool recurse=false) const
+ {
+ smoab::Range result;
+ this->Moab->get_entities_by_handle(entity,result,recurse);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
//Find all entities with a given tag. We don't use geom as a tag as that
//isn't a fast operation. Yes finding the intersection of geom entities and
//a material / boundary tag will be more work, but it is rarely done currently
@@ -179,60 +244,55 @@ public:
//----------------------------------------------------------------------------
//Find all entities from a given root of a given dimensionality
smoab::Range findEntitiesWithDimension(const smoab::EntityHandle root,
- int dimension) const
+ const int dimension,
+ bool recurse=false) const
{
typedef smoab::Range::const_iterator iterator;
smoab::Range result;
- this->Moab->get_entities_by_dimension(root,dimension,result);
+ this->Moab->get_entities_by_dimension(root,dimension,result,recurse);
-
- smoab::Range children;
- this->Moab->get_child_meshsets(root,children,0);
- for(iterator i=children.begin(); i !=children.end();++i)
+ if(recurse)
{
- this->Moab->get_entities_by_dimension(*i,dimension,result);
+ smoab::Range children;
+ this->Moab->get_child_meshsets(root,children,0);
+ for(iterator i=children.begin(); i !=children.end();++i)
+ {
+ this->Moab->get_entities_by_dimension(*i,dimension,result);
+ }
}
return result;
}
- //----------------------------------------------------------------------------
- smoab::Range findAdjacentEntities(const smoab::EntityHandle& entity,
- int dimension) const
- {
- const int adjType = static_cast<int>(smoab::INTERSECT);
- smoab::Range result;
- const bool create_if_missing = false;
- this->Moab->get_adjacencies(&entity,
- 1,
- dimension,
- create_if_missing,
- result,
- adjType);
- return result;
- }
- //----------------------------------------------------------------------------
- smoab::Range findAdjacentEntities(const smoab::Range& range,
- int dimension,
- const smoab::adjacency_type type = smoab::UNION) const
+ //----------------------------------------------------------------------------
+ smoab::Range findHighestDimensionEntities(const smoab::EntityHandle& entity,
+ bool recurse=false) const
{
- //the smoab and moab adjacent intersection enums are in the same order
- const int adjType = static_cast<int>(type);
- smoab::Range result;
- const bool create_if_missing = false;
- this->Moab->get_adjacencies(range,dimension,
- create_if_missing,
- result,
- adjType);
+ //the goal is to load all entities that are not entity sets of this
+ //node, while also subsetting by the highest dimension
- return result;
+ //lets find the entities of only the highest dimension
+ int num_ents=0;
+ int dim=3;
+ while(num_ents<=0&&dim>0)
+ {
+ this->Moab->get_number_entities_by_dimension(entity,dim,num_ents,recurse);
+ --dim;
+ }
+ ++dim; //reincrement to correct last decrement
+ if(num_ents > 0)
+ {
+ //we have found entities of a given dimension
+ return this->findEntitiesWithDimension(entity,dim,recurse);
+ }
+ return smoab::Range();
}
//----------------------------------------------------------------------------
//Find all elements in the database that have children and zero parents.
//this doesn't find
- smoab::Range findEntityRootParents(smoab::EntityHandle const& root) const
+ smoab::Range findEntityRootParents(const smoab::EntityHandle& root) const
{
smoab::Range parents;
@@ -258,7 +318,7 @@ public:
//----------------------------------------------------------------------------
//finds entities that have zero children and zero parents
- smoab::Range findDetachedEntities(moab::EntityHandle const& root) const
+ smoab::Range findDetachedEntities(const moab::EntityHandle& root) const
{
smoab::Range detached;
@@ -284,7 +344,7 @@ public:
//----------------------------------------------------------------------------
//find all children of the entity passed in that has multiple parents
- smoab::Range findEntitiesWithMultipleParents(smoab::EntityHandle const& root)
+ smoab::Range findEntitiesWithMultipleParents(const smoab::EntityHandle& root) const
{
smoab::Range multipleParents;
typedef moab::Range::const_iterator iterator;
@@ -305,8 +365,114 @@ public:
}
//----------------------------------------------------------------------------
+ //find all entities that are adjacent to a single entity
+ smoab::Range findAdjacencies(const smoab::EntityHandle& entity,
+ int dimension) const
+ {
+ const int adjType = static_cast<int>(smoab::INTERSECT);
+ smoab::Range result;
+ const bool create_if_missing = false;
+ this->Moab->get_adjacencies(&entity,
+ 1,
+ dimension,
+ create_if_missing,
+ result,
+ adjType);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ smoab::Range findAdjacencies(const smoab::Range& range,
+ int dimension,
+ const smoab::adjacency_type type = smoab::UNION) const
+ {
+ //the smoab and moab adjacent intersection enums are in the same order
+ const int adjType = static_cast<int>(type);
+ smoab::Range result;
+ const bool create_if_missing = false;
+ this->Moab->get_adjacencies(range,dimension,
+ create_if_missing,
+ result,
+ adjType);
+
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ //create adjacencies, only works when the dimension requested is lower than
+ //dimension of the range of entities
+ smoab::Range createAdjacencies(const smoab::Range& range,
+ int dimension,
+ const smoab::adjacency_type type = smoab::UNION) const
+ {
+ //the smoab and moab adjacent intersection enums are in the same order
+ const int adjType = static_cast<int>(type);
+ smoab::Range result;
+ const bool create_if_missing = true;
+ this->Moab->get_adjacencies(range,dimension,
+ create_if_missing,
+ result,
+ adjType);
+
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ int numChildMeshSets(const smoab::EntityHandle& root) const
+ {
+ int numChildren;
+ this->Moab->num_child_meshsets(root,&numChildren);
+ return numChildren;
+ }
+
+ //----------------------------------------------------------------------------
+ smoab::Range getChildSets(const smoab::EntityHandle& root) const
+ {
+ smoab::Range children;
+ this->Moab->get_child_meshsets(root,children,0);
+ return children;
+ }
+
+ //----------------------------------------------------------------------------
+ //remove a collection of entities from the database
+ void remove(smoab::Range const& toDelete) const
+ {
+ this->Moab->delete_entities(toDelete);
+ }
+
+ //----------------------------------------------------------------------------
+ //a entityHandle with value zero means no side element was found
+ smoab::EntityHandle sideElement(smoab::EntityHandle const& cell,
+ int dim, int side) const
+ {
+ smoab::EntityHandle result(0);
+ this->Moab->side_element(cell,dim,side,result);
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
+ //returns all the existing side elements of a cell, elements that
+ //are zero mean that side element doesn't exist
+ std::vector<smoab::EntityHandle> sideElements(
+ smoab::EntityHandle const& cell,
+ int dim) const
+ {
+ const EntityType volumeCellType = this->Moab->type_from_handle(cell);
+ const int numSides = static_cast<int>(moab::CN::NumSubEntities(
+ volumeCellType, dim));
+
+ std::vector<smoab::EntityHandle> result(numSides);
+ for (int side = 0; side < numSides; ++side)
+ {
+ smoab::EntityHandle *sideElem = &result[side]; //get memory of vector
+ this->Moab->side_element(cell,dim,side,*sideElem);
+ }
+ return result;
+ }
+
+ //----------------------------------------------------------------------------
//prints all elements in a range objects
- void printRange(smoab::Range const& range)
+ void printRange(const smoab::Range& range) const
{
typedef Range::const_iterator iterator;
for(iterator i=range.begin(); i!=range.end(); ++i)
@@ -315,12 +481,25 @@ public:
this->Moab->list_entity(*i);
}
}
-
friend class smoab::DataSetConverter;
+ friend class smoab::detail::LoadGeometry;
+ friend class smoab::detail::LoadPoly;
private:
moab::Interface* Moab;
};
+//----------------------------------------------------------------------------
+void RangeToVector(const smoab::Range &range,
+ std::vector<smoab::EntityHandle>& vector )
+{
+ vector.reserve(range.size());
+ std::copy(range.begin(),
+ range.end(),
+ std::back_inserter(vector));
+}
+
+
+
}
#endif
diff --git a/tools/vtkMOABReaderNew/detail/CellTypeToType.h b/tools/vtkMOABReaderNew/detail/CellTypeToType.h
new file mode 100644
index 0000000..b2955e4
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/CellTypeToType.h
@@ -0,0 +1,147 @@
+#ifndef __smoab_detail_CellTypeToType_h
+#define __smoab_detail_CellTypeToType_h
+
+#include "vtkCellType.h"
+#include <algorithm>
+
+namespace smoab{ namespace detail{
+
+template<int N> struct QuadratricOrdering{};
+
+template<> struct QuadratricOrdering<VTK_QUADRATIC_WEDGE>
+{
+ static const int NUM_VERTS = 15;
+ void reorder(vtkIdType* connectivity) const
+ {
+ std::swap_ranges(connectivity+9,connectivity+12,connectivity+12);
+ }
+};
+
+template<> struct QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON>
+{
+ static const int NUM_VERTS = 27;
+ void reorder(vtkIdType* connectivity) const
+ {
+ std::swap_ranges(connectivity+12,connectivity+16,connectivity+16);
+
+ //move 20 to 22
+ //move 22 to 23
+ //move 23 to 20
+
+ //swap 20 with 22
+ std::swap(connectivity[20],connectivity[23]);
+
+ //swap 22 with 23
+ std::swap(connectivity[22],connectivity[23]);
+ }
+};
+
+template<typename QuadraticOrdering>
+void FixQuadraticIdOrdering(vtkIdType* connectivity, vtkIdType numCells,
+ QuadraticOrdering& ordering)
+{
+ //skip the first index that holds the length of the cells
+ //if we skip it once here, and than properly increment it makes the code
+ //far easier
+ connectivity+=1;
+ for(vtkIdType i=0; i < numCells; ++i)
+ {
+ ordering.reorder(connectivity);
+ connectivity += ordering.NUM_VERTS + 1;
+ }
+}
+
+
+int vtkCellType(moab::EntityType t, int &num_connect)
+ {
+ int ctype = -1;
+ switch (t)
+ {
+ case moab::MBEDGE:
+ if (num_connect == 2) ctype = VTK_LINE;
+ else if (num_connect == 3) ctype = VTK_QUADRATIC_EDGE;
+ break;
+ case moab::MBTRI:
+ if (num_connect == 3) ctype = VTK_TRIANGLE;
+ else if (num_connect == 6) ctype = VTK_QUADRATIC_TRIANGLE;
+ else if (num_connect == 7) ctype = VTK_BIQUADRATIC_TRIANGLE;
+ break;
+ case moab::MBQUAD:
+ if (num_connect == 4) ctype = VTK_QUAD;
+ else if (num_connect == 8) ctype = VTK_QUADRATIC_QUAD;
+ else if (num_connect == 9) ctype = VTK_BIQUADRATIC_QUAD;
+ break;
+ case moab::MBPOLYGON:
+ if (num_connect == 4) ctype = VTK_POLYGON;
+ break;
+ case moab::MBTET:
+ if (num_connect == 4) ctype = VTK_TETRA;
+ else if (num_connect == 10) ctype = VTK_QUADRATIC_TETRA;
+ break;
+ case moab::MBPYRAMID:
+ if (num_connect == 5) ctype = VTK_PYRAMID;
+ else if (num_connect == 13) ctype = VTK_QUADRATIC_PYRAMID;
+ break;
+ case moab::MBPRISM:
+ if (num_connect == 6) ctype = VTK_WEDGE;
+ else if (num_connect == 15) ctype = VTK_QUADRATIC_WEDGE;
+ break;
+ case moab::MBHEX:
+ if (num_connect == 8) ctype = VTK_HEXAHEDRON;
+ else if (num_connect == 20) ctype = VTK_QUADRATIC_HEXAHEDRON;
+ else if (num_connect == 21) ctype = VTK_QUADRATIC_HEXAHEDRON, num_connect = 20;
+ else if (num_connect == 27) ctype = VTK_TRIQUADRATIC_HEXAHEDRON;
+ break;
+ default:
+ ctype = -1;
+ break;
+ }
+ return ctype;
+ }
+
+int vtkLinearCellType(moab::EntityType t, int &num_connect)
+ {
+ int ctype = -1;
+ switch (t)
+ {
+ case moab::MBEDGE:
+ ctype = VTK_LINE;
+ num_connect = 2;
+ break;
+ case moab::MBTRI:
+ ctype = VTK_TRIANGLE;
+ num_connect = 3;
+ break;
+ case moab::MBQUAD:
+ ctype = VTK_QUAD;
+ num_connect = 4;
+ break;
+ case moab::MBPOLYGON:
+ ctype = VTK_POLYGON;
+ num_connect = 4;
+ break;
+ case moab::MBTET:
+ ctype = VTK_TETRA;
+ num_connect = 4;
+ break;
+ case moab::MBPYRAMID:
+ ctype = VTK_PYRAMID;
+ num_connect = 5;
+ break;
+ case moab::MBPRISM:
+ ctype = VTK_WEDGE;
+ num_connect = 6;
+ break;
+ case moab::MBHEX:
+ ctype = VTK_HEXAHEDRON;
+ num_connect = 8;
+ break;
+ default:
+ break;
+ }
+ return ctype;
+ }
+
+} } //namespace smaob::detail
+
+#endif // CELLTYPETOTYPE_H
diff --git a/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h b/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h
new file mode 100644
index 0000000..28319df
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/ContinousCellInfo.h
@@ -0,0 +1,18 @@
+#ifndef __smoab_detail_ContinousCellInfo_h
+#define __smoab_detail_ContinousCellInfo_h
+
+
+namespace smoab { namespace detail {
+
+struct ContinousCellInfo
+{
+ int type;
+ int numVerts;
+ int numUnusedVerts;
+ int numCells;
+
+};
+
+} } //namespace smoab::detail
+
+#endif
diff --git a/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h b/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h
new file mode 100644
index 0000000..1b94548
--- /dev/null
+++ b/tools/vtkMOABReaderNew/detail/LinearCellConnectivity.h
@@ -0,0 +1,219 @@
+#ifndef __smoab_LinearCellConnectivity_h
+#define __smoab_LinearCellConnectivity_h
+
+#include "CellTypeToType.h"
+#include "ContinousCellInfo.h"
+
+#include <algorithm>
+#include <vector>
+
+namespace smoab { namespace detail {
+
+namespace internal
+{
+ //we want a subset of the real connetivity array,
+ //this does that for use with a super easy wrapper
+ struct SubsetArray
+ {
+ SubsetArray(EntityHandle* realConn,
+ int numCells,
+ int currentVertsPerCell,
+ int newVertsPerCell):
+ Array()
+ {
+ const int size = numCells*newVertsPerCell;
+ this->Array.reserve(size);
+ if(currentVertsPerCell == newVertsPerCell)
+ {
+ std::copy(realConn,realConn+size, std::back_inserter(this->Array));
+ }
+ else
+ {
+ //skip copy only the first N points which we want
+ //since moab stores linear points first per cell
+ EntityHandle *pos = realConn;
+ for(int i=0; i < numCells;++i)
+ {
+ std::copy(pos,pos+newVertsPerCell,std::back_inserter(this->Array));
+ pos += currentVertsPerCell;
+ }
+ }
+ }
+ typedef std::vector<EntityHandle>::const_iterator const_iterator;
+ typedef std::vector<EntityHandle>::iterator iterator;
+
+ const_iterator begin() const { return this->Array.begin(); }
+ iterator begin() { return this->Array.begin(); }
+
+ const_iterator end() const { return this->Array.end(); }
+ iterator end(){ return this->Array.end(); }
+
+ private:
+ std::vector<EntityHandle> Array;
+ };
+}
+
+class LinearCellConnectivity
+{
+public:
+
+ LinearCellConnectivity(smoab::Range const& cells, moab::Interface* moab):
+ Connectivity(),
+ UniquePoints(),
+ Info()
+ {
+ int count = 0;
+ const std::size_t cellSize=cells.size();
+ while(count != cellSize)
+ {
+ EntityHandle* connectivity;
+ int numVerts=0, iterationCount=0;
+ //use the highly efficent calls, since we know that are of the same dimension
+ moab->connect_iterate(cells.begin()+count,
+ cells.end(),
+ connectivity,
+ numVerts,
+ iterationCount);
+ //if we didn't read anything, break!
+ if(iterationCount == 0)
+ {
+ break;
+ }
+
+ //identify the cell type that we currently have,
+ //store that along with the connectivity in a temp storage vector
+ const moab::EntityType type = moab->type_from_handle(*cells.begin()+count);
+
+ int vtkNumVerts;
+ int vtkCellType = smoab::detail::vtkLinearCellType(type,vtkNumVerts);
+
+ ContinousCellInfo info = { vtkCellType, vtkNumVerts, 0, iterationCount };
+ this->Info.push_back(info);
+
+
+ //we need to copy only a subset of the connectivity array
+ internal::SubsetArray conn(connectivity,iterationCount,numVerts,vtkNumVerts);
+ this->Connectivity.push_back(conn);
+
+ count += iterationCount;
+ }
+ }
+
+ //----------------------------------------------------------------------------
+ void compactIds(vtkIdType& numCells, vtkIdType& connectivityLength)
+ {
+ //converts all the ids to be ordered starting at zero, and also
+ //keeping the orginal logical ordering. Stores the result of this
+ //operation in the unstrucutred grid that is passed in
+
+ //lets determine the total length of the connectivity
+ connectivityLength = 0;
+ numCells = 0;
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i)
+ {
+ connectivityLength += (*i).numCells * (*i).numVerts;
+ numCells += (*i).numCells;
+ }
+
+ this->UniquePoints.reserve(connectivityLength);
+
+ this->copyConnectivity(this->UniquePoints);
+ std::sort(this->UniquePoints.begin(),this->UniquePoints.end());
+
+ typedef std::vector<EntityHandle>::iterator EntityIterator;
+ EntityIterator newEnd = std::unique(this->UniquePoints.begin(),
+ this->UniquePoints.end());
+
+ const std::size_t newSize = std::distance(this->UniquePoints.begin(),newEnd);
+ this->UniquePoints.resize(newSize);
+ }
+
+ //----------------------------------------------------------------------------
+ void moabPoints(smoab::Range& range) const
+ {
+ //from the documentation a reverse iterator is the fastest way
+ //to insert into a range.
+ std::copy(this->UniquePoints.rbegin(),
+ this->UniquePoints.rend(),
+ moab::range_inserter(range));
+ }
+
+ //----------------------------------------------------------------------------
+ //copy the connectivity from the moab held arrays to the user input vector
+ void copyConnectivity(std::vector<EntityHandle>& output) const
+ {
+ //walk the info to find the length of each sub connectivity array,
+ //and insert them into the vector, ordering is implied by the order
+ //the connecitivy sub array are added to this class
+ ConnConstIterator c = this->Connectivity.begin();
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i,++c)
+ {
+ //remember our Connectivity is a vector of pointers whose
+ //length is held in the info vector.
+ const int numUnusedPoints = (*i).numUnusedVerts;
+ const int connLength = (*i).numCells * (*i).numVerts;
+ std::copy(c->begin(),c->end(),std::back_inserter(output));
+ }
+ }
+
+ //copy the information from this contianer to a vtk cell array, and
+ //related lookup information
+ void copyToVtkCellInfo(vtkIdType* cellArray) const
+ {
+ ConnConstIterator c = this->Connectivity.begin();
+ for(InfoConstIterator i = this->Info.begin();
+ i != this->Info.end();
+ ++i, ++c)
+ {
+ //for this group of the same cell type we need to fill the cellTypes
+ const int numCells = (*i).numCells;
+ const int numVerts = (*i).numVerts;
+
+ //for each cell in this collection that have the same type
+ //grab the raw array now, so we can properly increment for each vert in each cell
+ internal::SubsetArray::const_iterator moabConnectivity = c->begin();
+ for(int j=0;j < numCells; ++j)
+ {
+ //cell arrays start and end are different, since we
+ //have to account for element that states the length of each cell
+ cellArray[0]=numVerts;
+
+
+ for(int k=0; k < numVerts; ++k, ++moabConnectivity )
+ {
+ //this is going to be a root of some failures when we start
+ //reading really large datasets under 32bit.
+
+
+ //fyi, don't use a range ds for unique points, distance
+ //function is horribly slow they need to override it
+ EntityConstIterator result = std::lower_bound(
+ this->UniquePoints.begin(),
+ this->UniquePoints.end(),
+ *moabConnectivity);
+ std::size_t newId = std::distance(this->UniquePoints.begin(),
+ result);
+ cellArray[k+1] = static_cast<vtkIdType>(newId);
+ }
+ cellArray += numVerts+1;
+ }
+ }
+ }
+
+private:
+ std::vector<internal::SubsetArray> Connectivity;
+ std::vector<EntityHandle> UniquePoints;
+
+ std::vector<detail::ContinousCellInfo> Info;
+
+ typedef std::vector<EntityHandle>::const_iterator EntityConstIterator;
+ typedef std::vector<internal::SubsetArray>::const_iterator ConnConstIterator;
+ typedef std::vector<detail::ContinousCellInfo>::const_iterator InfoConstIterator;
+};
+} } //namespace smoab::detail
+
+#endif // __smoab_LinearCellConnectivity_h
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list