[MOAB-dev] r1797 - MOAB/trunk
bckiedro at mcs.anl.gov
bckiedro at mcs.anl.gov
Tue May 6 15:54:34 CDT 2008
Author: bckiedro
Date: 2008-05-06 15:54:34 -0500 (Tue, 06 May 2008)
New Revision: 1797
Added:
MOAB/trunk/MBElemUtil.cpp
MOAB/trunk/MBElemUtil.hpp
MOAB/trunk/ReadIDEAS.cpp
MOAB/trunk/ReadIDEAS.hpp
Modified:
MOAB/trunk/MBGeomUtil.cpp
MOAB/trunk/MBGeomUtil.hpp
MOAB/trunk/MBReaderWriterSet.cpp
MOAB/trunk/Makefile.am
Log:
05/06/2008 -- BCK
Added the following capabilities:
* The kd-tree can now handle MBHEX elements
* Support for reading IDEAS .unv file format
MBGeomUtil.cpp : in box_elem_overlap, added a case for MBHEX
added box_hex_overlap routine
added point_in_trilinear_hex routine (two versions, one with and without acceleration)
MBElemUtil.cpp : created this new file for code relating to element utilities
currently contains nat_coords_trilinear_hex that gets natural coordinates of a MBHEX
ReadIDEAS.cpp : created this new file as part of reader suite to read in IDEAS (.unv) files
MBReaderWriterSet.cpp : added register_factory for IDEAS format
The following header files have been added or modifed
MBGeomUtil.hpp
MBElemUtil.hpp
ReadIDEAS.hpp
Added: MOAB/trunk/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/MBElemUtil.cpp (rev 0)
+++ MOAB/trunk/MBElemUtil.cpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -0,0 +1,133 @@
+#include <iostream>
+#include <assert.h>
+
+#include "MBElemUtil.hpp"
+
+//
+// nat_coords_trilinear_hex
+// Given an MBEntityHandle defining a MBHEX element defined by a
+// trilinear basis and a set of xyz points in space, it finds the
+// cooresponding natural coordinates and returns them to the pointer
+// nat. The method uses an iterative technique, guessing initial
+// coordinates {0, 0, 0} and calculating new coordinates based on the
+// ones just calculated.
+//
+void MBElemUtil::nat_coords_trilinear_hex(MBCartVect hex[8],
+ MBCartVect xyz,
+ MBCartVect &ncoords,
+ double etol)
+{
+
+ MBCartVect nat(0.);
+ double A[8]; double B[8]; double C[8]; double D[8];
+ double Ax, By, Cz;
+ double err = 1.e37;
+
+ double xt, yt, zt;
+ double nxi, neta, nmu;
+ double pxi, peta, pmu;
+ double tmp;
+
+ // Iterative estimate of natural coordinates
+ while ( err > etol ) {
+
+ // Estimate the xi-coordinate
+ A[0] = (1 - nat[1]) * (1 - nat[2]);
+ A[1] = A[0];
+ A[2] = (1 + nat[1]) * (1 - nat[2]);
+ A[3] = A[2];
+ A[4] = (1 - nat[1]) * (1 + nat[2]);
+ A[5] = A[4];
+ A[6] = (1 + nat[1]) * (1 + nat[2]);
+ A[7] = A[6];
+
+ Ax = 0;
+ for (unsigned int i = 0; i < 8; i++) {
+ Ax = Ax + A[i]*hex[i][0];
+ }
+ nat[0] = (8*xyz[0] - Ax ) /
+ (-A[0]*hex[0][0] + A[1]*hex[1][0] + A[2]*hex[2][0] - A[3]*hex[3][0]
+ -A[4]*hex[4][0] + A[5]*hex[5][0] + A[6]*hex[6][0] - A[7]*hex[7][0]);
+
+ // Estimate the eta-coordinate
+ B[0] = (1 - nat[0]) * (1 - nat[2]);
+ B[1] = (1 + nat[0]) * (1 - nat[2]);
+ B[2] = B[1];
+ B[3] = B[0];
+ B[4] = (1 - nat[0]) * (1 + nat[2]);
+ B[5] = (1 + nat[0]) * (1 + nat[2]);
+ B[6] = B[5];
+ B[7] = B[4];
+
+ By = 0;
+ for (unsigned int i = 0; i < 8; i++) {
+ By = By + B[i]*hex[i][1];
+ }
+ nat[1] = (8*xyz[1] - By ) /
+ (-B[0]*hex[0][1] - B[1]*hex[1][1] + B[2]*hex[2][1] + B[3]*hex[3][1]
+ -B[4]*hex[4][1] - B[5]*hex[5][1] + B[6]*hex[6][1] + B[7]*hex[7][1]);
+
+ // Estimate the mu-coordinate
+ C[0] = (1 - nat[0]) * (1 - nat[1]);
+ C[1] = (1 + nat[0]) * (1 - nat[1]);
+ C[2] = (1 + nat[0]) * (1 + nat[1]);
+ C[3] = (1 - nat[0]) * (1 + nat[1]);
+ C[4] = C[0];
+ C[5] = C[1];
+ C[6] = C[2];
+ C[7] = C[3];
+
+ Cz = 0;
+ for (unsigned int i = 0; i < 8; i++) {
+ Cz = Cz + C[i]*hex[i][2];
+ }
+ nat[2] = (8*xyz[2] - Cz ) /
+ (-C[0]*hex[0][2] - C[1]*hex[1][2] - C[2]*hex[2][2] - C[3]*hex[3][2]
+ +C[4]*hex[4][2] + C[5]*hex[5][2] + C[6]*hex[6][2] + C[7]*hex[7][2]);
+
+
+ // Shortcut variables...
+ nxi = 1 - nat[0];
+ neta = 1 - nat[1];
+ nmu = 1 - nat[2];
+ pxi = 1 + nat[0];
+ peta = 1 + nat[1];
+ pmu = 1 + nat[2];
+ D[0] = nxi * neta * nmu;
+ D[1] = pxi * neta * nmu;
+ D[2] = pxi * peta * nmu;
+ D[3] = nxi * peta * nmu;
+ D[4] = nxi * neta * pmu;
+ D[5] = pxi * neta * pmu;
+ D[6] = pxi * peta * pmu;
+ D[7] = nxi * peta * pmu;
+
+ // Compute corresponding estimates for x, y, and z to check
+ xt = 0.125 * ( D[0] * hex[0][0] + D[1] * hex[1][0] +
+ D[2] * hex[2][0] + D[3] * hex[3][0] +
+ D[4] * hex[4][0] + D[5] * hex[5][0] +
+ D[6] * hex[6][0] + D[7] * hex[7][0] );
+
+ yt = 0.125 * ( D[0] * hex[0][1] + D[1] * hex[1][1] +
+ D[2] * hex[2][1] + D[3] * hex[3][1] +
+ D[4] * hex[4][1] + D[5] * hex[5][1] +
+ D[6] * hex[6][1] + D[7] * hex[7][1] );
+
+ zt = 0.125 * ( D[0] * hex[0][2] + D[1] * hex[1][2] +
+ D[2] * hex[2][2] + D[3] * hex[3][2] +
+ D[4] * hex[4][2] + D[5] * hex[5][2] +
+ D[6] * hex[6][2] + D[7] * hex[7][2] );
+
+ // Compute error
+ err = fabs(xt - xyz[0]);
+ tmp = fabs(yt - xyz[1]);
+ if (tmp > err) err = tmp;
+ tmp = fabs(zt - xyz[2]);
+ if (tmp > err) err = tmp;
+
+ }
+
+ ncoords = nat;
+ return;
+
+}
Added: MOAB/trunk/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/MBElemUtil.hpp (rev 0)
+++ MOAB/trunk/MBElemUtil.hpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -0,0 +1,12 @@
+#include "MBCore.hpp"
+#include "MBCartVect.hpp"
+
+class MBElemUtil {
+
+ public:
+
+ static void nat_coords_trilinear_hex(MBCartVect[8],
+ MBCartVect,
+ MBCartVect&,
+ double);
+};
Modified: MOAB/trunk/MBGeomUtil.cpp
===================================================================
--- MOAB/trunk/MBGeomUtil.cpp 2008-05-06 19:05:58 UTC (rev 1796)
+++ MOAB/trunk/MBGeomUtil.cpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -21,9 +21,11 @@
#include "MBCartVect.hpp"
#include "MBCN.hpp"
#include "MBGeomUtil.hpp"
+#include "MBElemUtil.hpp"
#include <cmath>
#include <algorithm>
#include <assert.h>
+#include <iostream>
#ifdef _MSC_VER
# include <float.h>
@@ -269,9 +271,12 @@
const MBCartVect& center,
const MBCartVect& dims )
{
+
switch (elem_type) {
case MBTRI:
return box_tri_overlap( elem_corners, center, dims );
+ case MBHEX:
+ return box_hex_overlap( elem_corners, center, dims );
case MBPOLYGON:
case MBPOLYHEDRON:
assert(false);
@@ -654,4 +659,195 @@
return closest % closest < tolerance * tolerance;
}
+bool box_hex_overlap( const MBCartVect hex_vertices[8],
+ const MBCartVect& box_center,
+ const MBCartVect& box_dims)
+{
+
+ // Mapping of vertices to each face on MBHEX
+ const unsigned int facids[] = {0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6,
+ 3, 0, 4, 7, 3, 2, 1, 0, 4, 5, 6, 7};
+ // Mapping of vertices to each edge on MBHEX
+ const unsigned int edgids[] = {0, 1, 1, 2, 2, 3, 3, 0,
+ 0, 4, 1, 5, 2, 6, 3, 7,
+ 4, 5, 5, 6, 6, 7, 7, 4};
+
+ const double eps = 1.e-10;
+
+ // Center the hex such that origin located at the box
+ MBCartVect hexv[8];
+ for (unsigned int i=0; i < 8; i++) {
+ hexv[i] = hex_vertices[i] - box_center;
+ }
+
+ // Check hex against box normals...
+ // Need to check both positive and negative faces
+ double t;
+ unsigned int positive = 0, negative = 0;
+ unsigned int out_plus = 0, out_minus = 0, inside = 0;
+ bool test;
+ // Negative face, loop over all vertices in hex
+ for (unsigned int i=0; i < 3; i++) {
+ // vxyz = box_center[i] + box_dims[i];
+ out_plus = 0; out_minus = 0; inside = 0;
+ // Test each vertex on hex...
+ for (unsigned int j=0; j < 8; j++) {
+
+ const double bxd = box_dims[i] + eps;
+
+ if (hexv[j][i] > bxd ) out_plus++;
+ else if (hexv[j][i] < -bxd ) out_minus++;
+ else inside++;
+
+ test = ((inside) || ( out_plus && out_minus ));
+ if (test) break;
+ }
+ if (!test) return false;
+ }
+
+ // Construct the vertices of the box...
+ MBCartVect boxv[8];
+
+ boxv[0][0] = -box_dims[0]; boxv[0][1] = -box_dims[1]; boxv[0][2] = -box_dims[2];
+ boxv[1][0] = box_dims[0]; boxv[1][1] = -box_dims[1]; boxv[1][2] = -box_dims[2];
+ boxv[2][0] = box_dims[0]; boxv[2][1] = box_dims[1]; boxv[2][2] = -box_dims[2];
+ boxv[3][0] = -box_dims[0]; boxv[3][1] = box_dims[1]; boxv[3][2] = -box_dims[2];
+ boxv[4][0] = -box_dims[0]; boxv[4][1] = -box_dims[1]; boxv[4][2] = box_dims[2];
+ boxv[5][0] = box_dims[0]; boxv[5][1] = -box_dims[1]; boxv[5][2] = box_dims[2];
+ boxv[6][0] = box_dims[0]; boxv[6][1] = box_dims[1]; boxv[6][2] = box_dims[2];
+ boxv[7][0] = -box_dims[0]; boxv[7][1] = box_dims[1]; boxv[7][2] = box_dims[2];
+
+ // Check box against hex normals ...
+ // Loop over each face of hex
+ for (unsigned int i = 0; i < 6; i++ ) {
+ // Compute the normal
+ const MBCartVect midpt1 = 0.5 * (( hexv[(facids[4*i+2])] + hexv[(facids[4*i+3])] )
+ - ( hexv[(facids[4*i])] + hexv[(facids[4*i+1])] ));
+ const MBCartVect midpt2 = 0.5 * (( hexv[(facids[4*i+3])] + hexv[(facids[4*i])] )
+ - ( hexv[(facids[4*i+1])] + hexv[(facids[4*i+2])] ));
+ const MBCartVect normal = midpt1 * midpt2;
+
+ // Loop over each vertex in the box
+ positive = 0; negative = 0;
+ for (unsigned int j=0; j < 8; j++) {
+ // Take dot product of the vector and normal
+ t = ( boxv[j] - hexv[(facids[4*i])] ) % normal;
+
+ // Do the comparison
+ if (t > eps ) positive++;
+ else negative++;
+ if (positive && negative) break;
+ }
+ // std::cout << positive << " " << negative << std::endl;
+ test = ( (positive && !negative) ? false : true );
+ if (!test) return false;
+ }
+
+ // Edge check ...
+ // For box, only need to check three edges due to orthogonality
+ int side1, side2;
+ MBCartVect edge_cross;
+
+ for (unsigned int d = 0; d < 3; d++ ) {
+ // const unsigned int d0 = d % 3;
+ const unsigned int d1 = (d + 1) % 3;
+ const unsigned int d2 = (d + 2) % 3;
+
+ for (unsigned int i = 0; i < 12; i++) {
+ // edge_cross[d0] = 0.;
+ edge_cross[d1] = hexv[(edgids[2*i])][d2] - hexv[(edgids[2*i+1])][d2];
+ edge_cross[d2] = hexv[(edgids[2*i+1])][d1] - hexv[(edgids[2*i])][d1];
+
+ if ((edge_cross[d1]*edge_cross[d1] + edge_cross[d2]*edge_cross[d2]) < eps) continue;
+
+ side1 = side2 = 0;
+
+ // Hex edge and box vertex tests
+ positive = 0; negative = 0; test = false;
+ for (unsigned int j = 0; j < 8; j++) {
+ // t = edge_cross % ( boxv[j] - hexv[(edgids[2*i])] );
+ t = edge_cross[d1] * ( boxv[j][d1] - hexv[(edgids[2*i])][d1] )
+ + edge_cross[d2] * ( boxv[j][d2] - hexv[(edgids[2*i])][d2] );
+ if (t > eps) positive++;
+ else if (t < -eps) negative++;
+
+ if (positive && negative) {
+ test = true;
+ break;
+ }
+ }
+ if (!test) side1 = ( (positive && !negative) ? +1 : -1 );
+ if (side1 == 0) continue;
+
+ // Hex edge and hex vertex tests
+ positive = 0; negative = 0; test = false;
+ for (unsigned int j = 0; j < 8; j++) {
+ // t = edge_cross % ( hexv[j] - hexv[(edgids[2*i])] );
+ t = edge_cross[d1] * ( boxv[j][d1] - hexv[(edgids[2*i])][d1] )
+ + edge_cross[d2] * ( boxv[j][d2] - hexv[(edgids[2*i])][d2] );
+ if (t > eps) positive++;
+ else if (t < -eps) negative++;
+
+ if (positive && negative) {
+ test = true;
+ break;
+ }
+ }
+ if (!test) side2 = ( (positive && !negative) ? +1 : -1 );
+ if (side2 == 0) continue;
+
+ // Test if on opposite sides
+ if (side1 * side2 < 0 ) return false;
+
+ }
+
+ }
+
+ // All possibilities exhausted, must overlap ...
+ return true;
+}
+
+
+bool point_in_trilinear_hex(MBCartVect hex[8],
+ MBCartVect xyz,
+ double etol)
+{
+
+ const double one = 1.000001;
+
+ MBCartVect nat(0.);
+ MBElemUtil::nat_coords_trilinear_hex(hex, xyz, nat, etol);
+
+ for (unsigned int i = 0; i < 3; i++) {
+ if ((nat[i] > one) || (nat[i] < -one)) return false;
+ }
+
+ return true;
+
+}
+
+
+bool point_in_trilinear_hex(MBCartVect hex[8],
+ MBCartVect xyz,
+ MBCartVect box_min, MBCartVect box_max,
+ double etol)
+{
+
+ const double one = 1.000001;
+
+ if ((xyz[0] < box_min[0]) || (xyz[0] > box_max[0])) return false;
+ if ((xyz[1] < box_min[1]) || (xyz[1] > box_max[1])) return false;
+ if ((xyz[2] < box_min[2]) || (xyz[2] > box_max[2])) return false;
+
+ MBCartVect nat(0.);
+ MBElemUtil::nat_coords_trilinear_hex(hex, xyz, nat, etol);
+
+ for (unsigned int i = 0; i < 3; i++) {
+ if ((nat[i] > one) || (nat[i] < -one)) return false;
+ }
+
+ return true;
+
+}
+
} // namespace MBGeoemtry
Modified: MOAB/trunk/MBGeomUtil.hpp
===================================================================
--- MOAB/trunk/MBGeomUtil.hpp 2008-05-06 19:05:58 UTC (rev 1796)
+++ MOAB/trunk/MBGeomUtil.hpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -226,6 +226,36 @@
MBCartVect& closest_out,
int& closest_topo );
+// Finds whether or not a box defined by the center and the half
+// width intersects a trilinear hex defined by its eight vertices.
+bool box_hex_overlap( const MBCartVect hexv[8],
+ const MBCartVect& box_center,
+ const MBCartVect& box_dims);
+
+//
+// point_in_trilinear_hex
+// Tests if a point in xyz space is within a hex element defined with
+// its eight vertex points forming a trilinear basis function. Computes
+// the natural coordinates with respect to the hex of the xyz point
+// and checks if each are between +/-1. If anyone is outside the range
+// the function returns false, otherwise it returns true.
+//
+bool point_in_trilinear_hex(MBCartVect hex[8],
+ MBCartVect xyz,
+ double etol);
+
+//
+// point_in_trilinear_hex
+// Tests if a point in xyz space is within a hex element defined with
+// its eight vertex points forming a trilinear basis function. Like the
+// test above except that it gets a bounding box as arguments to filter
+// the test for acceleration.
+//
+bool point_in_trilinear_hex(MBCartVect hex[8],
+ MBCartVect xyz,
+ MBCartVect box_min, MBCartVect box_max,
+ double etol);
+
} // namespace MBGeoemtry
#endif
Modified: MOAB/trunk/MBReaderWriterSet.cpp
===================================================================
--- MOAB/trunk/MBReaderWriterSet.cpp 2008-05-06 19:05:58 UTC (rev 1796)
+++ MOAB/trunk/MBReaderWriterSet.cpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -23,6 +23,7 @@
#include "ReadVtk.hpp"
#include "ReadSTL.hpp"
#include "ReadGmsh.hpp"
+#include "ReadIDEAS.hpp"
#include "Tqdcfr.hpp"
#include "WriteAns.hpp"
@@ -66,6 +67,8 @@
const char* exo_sufxs[] = { "exo", "exoII", "exo2", "g", "gen", NULL };
register_factory( ReadNCDF::factory, WriteNCDF::factory, "Exodus II", exo_sufxs, "EXODUS" );
#endif
+
+ register_factory( ReadIDEAS::factory, NULL, "IDEAS format", "unv", "UNV" );
register_factory( ReadVtk::factory, WriteVtk::factory, "Kitware VTK", "vtk", "VTK" );
Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am 2008-05-06 19:05:58 UTC (rev 1796)
+++ MOAB/trunk/Makefile.am 2008-05-06 20:54:34 UTC (rev 1797)
@@ -122,6 +122,7 @@
MBCN.cpp \
MBCNArrays.hpp \
MBCartVect.cpp \
+ MBElemUtil.cpp \
MBHandleUtils.cpp \
MBMatrix3.cpp \
MBMatrix3.hpp \
@@ -158,6 +159,8 @@
ReadSTL.hpp \
ReadVtk.cpp \
ReadVtk.hpp \
+ ReadIDEAS.cpp \
+ ReadIDEAS.hpp \
ScdElementData.cpp \
ScdElementData.hpp \
ScdVertexData.cpp \
@@ -216,6 +219,7 @@
MBCN_protos.h \
MBCartVect.hpp \
MBCore.hpp \
+ MBElemUtil.hpp \
MBEntityType.h \
MBEntityHandle.h \
MBError.hpp \
Added: MOAB/trunk/ReadIDEAS.cpp
===================================================================
--- MOAB/trunk/ReadIDEAS.cpp (rev 0)
+++ MOAB/trunk/ReadIDEAS.cpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -0,0 +1,195 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include "assert.h"
+
+#include "ReadIDEAS.hpp"
+#include "MBInterface.hpp"
+#include "MBInternals.hpp"
+#include "MBReadUtilIface.hpp"
+#include "FileTokenizer.hpp"
+#include "MBRange.hpp"
+
+MBReaderIface* ReadIDEAS::factory( MBInterface* iface )
+ { return new ReadIDEAS( iface ); }
+
+ReadIDEAS::ReadIDEAS(MBInterface* impl)
+ : MBI(impl)
+{
+ impl->query_interface("MBReadUtilIface", reinterpret_cast<void**>(&readMeshIface));
+}
+
+MBErrorCode ReadIDEAS::load_file(const char* fname,
+ MBEntityHandle& meshset,
+ const FileOptions& options,
+ const int* material_set_list,
+ int num_material_sets ) {
+
+ file.open( fname );
+
+ MBErrorCode rval;
+
+ rval = MBI->create_meshset(MESHSET_SET, mesh_handle);
+ if (MB_SUCCESS != rval) return rval;
+
+ char line[10000];
+ file.getline(line, 10000);
+ std::string s = line;
+ if (s.find("-1") > s.length()) return MB_FAILURE;
+
+ while (! file.eof() ) {
+ file.getline(line, 10000);
+ s = line;
+
+ unsigned int header_id = (unsigned int) strtol(line, NULL, 10);
+ switch (header_id) {
+ case VERTEX_LIST :
+ create_vertices();
+ break;
+ case MAKE_TETRAHEDRA :
+ create_tetrahedral_elements();
+ break;
+ default:
+ skip_header();
+ break;
+ }
+
+ }
+
+ meshset = mesh_handle;
+ file.close();
+ return MB_SUCCESS;
+
+}
+
+void ReadIDEAS::skip_header() {
+
+ // Go until finding a pair of -1 lines
+ char *ctmp;
+ char line[10000];
+ std::string s;
+
+ int end_of_block = 0;
+
+ long int il;
+
+ while (! file.eof() ) {
+ file.getline(line, 10000);
+
+ il = std::strtol(line, &ctmp, 10);
+ if (il == -1) {
+ s = ctmp+1;
+ if (s.empty()) end_of_block++;
+ }
+ else end_of_block = 0;
+
+ if (end_of_block >= 2) break;
+
+ }
+
+}
+
+
+void ReadIDEAS::create_vertices() {
+
+ // Read two lines: first has some data, second has coordinates
+ double *coords;
+ char line1[10000], line2[10000];
+ int il1, il2;
+ char *ctmp1, *ctmp2;
+ std::string s1, s2;
+
+ MBErrorCode rval;
+
+ int top_of_block = file.tellg();
+ unsigned int num_verts = 0;
+
+ while (! file.eof() ) {
+
+ file.getline(line1, 10000);
+ file.getline(line2, 10000);
+
+ // Check if we are at the end of the block
+ il1 = std::strtol(line1, &ctmp1, 10);
+ il2 = std::strtol(line2, &ctmp2, 10);
+ if ((il1 == -1) && (il2 == -1)) {
+ s1 = ctmp1+1;
+ s2 = ctmp2+1;
+ if ((s1.empty()) && (s2.empty())) break;
+ }
+ num_verts++;
+ }
+
+ file.seekg( top_of_block );
+ coords = new double [ 3*num_verts ];
+
+ for (unsigned int i = 0; i < num_verts; i++) {
+
+ file.getline(line1, 10000);
+ file.getline(line2, 10000);
+
+ // Get the doubles out of the 2nd line
+ coords[3*i ] = std::strtod(line2, &ctmp2);
+ coords[3*i+1] = std::strtod(ctmp2+1, &ctmp2);
+ coords[3*i+2] = std::strtod(ctmp2+1, NULL);
+
+ }
+
+ MBRange verts;
+ rval = MBI->create_vertices(coords, num_verts, verts);
+ assert( MB_SUCCESS == rval );
+
+ rval = MBI->add_entities( mesh_handle, verts );
+ assert( MB_SUCCESS == rval );
+
+ file.getline(line1, 10000);
+ file.getline(line2, 10000);
+
+}
+
+
+void ReadIDEAS::create_tetrahedral_elements() {
+
+ MBEntityHandle connect[4];
+ char line1[10000], line2[10000];
+ int il1, il2;
+ char *ctmp1, *ctmp2;
+ std::string s1, s2;
+
+ MBErrorCode rval;
+ MBEntityHandle handle;
+
+ MBRange verts;
+ rval = MBI->get_entities_by_type( mesh_handle, MBVERTEX, verts, true);
+ MBEntityHandle vstart = *(verts.begin());
+
+ while (! file.eof() ) {
+
+ file.getline(line1, 10000);
+ file.getline(line2, 10000);
+
+ // Check if we are at the end of the block
+ il1 = std::strtol(line1, &ctmp1, 10);
+ il2 = std::strtol(line2, &ctmp2, 10);
+ if ((il1 == -1) && (il2 == -1)) {
+ s1 = ctmp1+1;
+ s2 = ctmp2+1;
+ if ((s1.empty()) && (s2.empty())) break;
+ }
+
+ // Get the connectivity out of the 2nd line
+ connect[0] = vstart + strtol( line2, &ctmp2, 10) - 1;
+ connect[1] = vstart + strtol( ctmp2+1, &ctmp2, 10) - 1;
+ connect[2] = vstart + strtol( ctmp2+1, &ctmp2, 10) - 1;
+ connect[3] = vstart + strtol( ctmp2+1, &ctmp2, 10) - 1;
+
+ // Make the element
+ rval = MBI->create_element(MBTET, connect, 4, handle);
+ assert( MB_SUCCESS == rval );
+
+ rval = MBI->add_entities( mesh_handle, &handle, 1);
+ assert( MB_SUCCESS == rval );
+
+ }
+
+}
Added: MOAB/trunk/ReadIDEAS.hpp
===================================================================
--- MOAB/trunk/ReadIDEAS.hpp (rev 0)
+++ MOAB/trunk/ReadIDEAS.hpp 2008-05-06 20:54:34 UTC (rev 1797)
@@ -0,0 +1,51 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+
+#include "MBReaderIface.hpp"
+#include "MBInterface.hpp"
+
+#define VERTEX_LIST 2411
+#define MAKE_TETRAHEDRA 2412
+
+class MBReadUtilIface;
+
+class ReadIDEAS : public MBReaderIface
+{
+
+public:
+
+ static MBReaderIface* factory( MBInterface* );
+
+MBErrorCode load_file(const char* fname,
+ MBEntityHandle& meshset,
+ const FileOptions&,
+ const int* material_set_list,
+ int num_material_sets );
+
+ //! Constructor
+ ReadIDEAS(MBInterface* impl = NULL);
+
+ //! Destructor
+ virtual ~ReadIDEAS() {}
+
+protected:
+
+ void skip_header();
+ void create_vertices();
+ void create_tetrahedral_elements();
+
+private:
+
+ std::ifstream file;
+
+ // Read mesh interface
+ MBReadUtilIface* readMeshIface;
+
+ // MOAB Interface
+ MBInterface* MBI;
+
+ // Handle for the mesh
+ MBEntityHandle mesh_handle;
+
+};
More information about the moab-dev
mailing list