Commit 826c4386 authored by Edward Vigmond's avatar Edward Vigmond
Browse files

Now use libredblack to determine the region surfaces. Still slow but

tractable.
parent a0b10b75
......@@ -410,25 +410,32 @@ MultiPoint ** MultiPoint::isosurf( DATA_TYPE *dat, DATA_TYPE val, int &npoly,
/**
* return alist of lists of nodes defining the surface of a volume element
*
* the node table is filled in by row with the number of nodes and then
* the global node numbers
*
* \param ve volume element
* \param ft table to fill in with the face info
* \param ns number of surfaces
* \param nn number of nodes per surface
* \param nl table of local nodes for each surface
* \param ve volume element
*
* \return the number of faces
*/
vector<vector<int> > *
VolElement:: make_surf_nodelist( const int ns, int nn, const int **nl, int ve )
int
VolElement:: make_surf_nodelist( int ve, int **ft, const int ns, int nn, const int **nl )
{
vector<vector<int> >* vv = new vector<vector<int> >;
for( int s=0; s<ns; s++ ) {
vector<int> nlv;
for( int i=0; i<nn; i++ ) {
int lnode = *(reinterpret_cast<int*>(nl)+s*nn+i);
if( lnode !=-1 )
nlv.push_back(_node[ve*_ptsPerObj+lnode]);
int i;
int *rptr = reinterpret_cast<int*>(nl)+s*nn;
for( i=0; i<nn; i++ ) {
int lnode = *rptr++;
if( lnode == -1 )
break;
ft[s][i+1] = _node[ve*_ptsPerObj+lnode];
}
vv->push_back(nlv);
ft[s][0] = i;
}
return vv;
return ns;
}
......@@ -11,6 +11,9 @@
gzFile openFile( const char*, const char* );
#define MAX_NUM_SURF 6
#define MAX_NUM_SURF_NODES 4
class DrawingObj
{
public:
......@@ -201,12 +204,12 @@ class VolElement : public MultiPoint
void add( int *n, int r=-1 );
virtual void draw_out_face( int )=0;
virtual SurfaceElement *cut( char*, GLfloat*, Interpolator<DATA_TYPE>*&, int e=0 )=0;
virtual vector<vector<int> >* surfaces(int v=0) = 0;
virtual int surfaces(int **, int v=0) = 0;
protected:
int* _region; //!< region for each element
Connection* _edges; //!< egdes for drawing elements
SurfaceElement* planecut( char*, GLfloat*, Interpolator<DATA_TYPE>*&, int, const int [][2], int e );
vector<vector<int> >* make_surf_nodelist(const int,int, const int **, int);
int make_surf_nodelist(int, int**, const int, int, const int **);
};
......@@ -221,7 +224,7 @@ class Tetrahedral : public VolElement
virtual void draw_out_face( int );
virtual SurfaceElement* cut(char*,GLfloat*,Interpolator<DATA_TYPE>*&,int=0);
const int* iso_polys( unsigned int );
virtual vector<vector<int> >* surfaces(int v=0);
virtual int surfaces(int **, int v=0);
};
class Prism : public VolElement
......@@ -235,7 +238,7 @@ class Prism : public VolElement
virtual void draw_out_face( int );
virtual SurfaceElement* cut(char *,GLfloat*,Interpolator<DATA_TYPE>*&,int);
const int* iso_polys(unsigned int);
virtual vector<vector<int> >* surfaces(int v=0);
virtual int surfaces(int **, int v=0);
};
......@@ -250,7 +253,7 @@ class Hexahedron : public VolElement
virtual void draw_out_face( int );
virtual SurfaceElement* cut(char *,GLfloat*,Interpolator<DATA_TYPE>*&,int);
const int* iso_polys( unsigned int );
virtual vector<vector<int> >* surfaces(int v=0);
virtual int surfaces(int **, int v=0);
};
class Pyramid : public VolElement
......@@ -264,6 +267,6 @@ class Pyramid : public VolElement
virtual void draw_out_face( int );
virtual SurfaceElement* cut(char*, GLfloat*, Interpolator<DATA_TYPE>*&,int);
const int* iso_polys( unsigned int );
virtual vector<vector<int> >* surfaces(int v=0);
virtual int surfaces(int **, int v=0);
};
#endif
......@@ -522,13 +522,15 @@ const int* Hexahedron::iso_polys(unsigned int index)
/*
* return list of lists of nodes defining the bounding surfaces
*
* \param ft face table to be filled in
* \param v element number
*
* \param pointer to a vector of vectors
*/
vector<vector<int> > *
Hexahedron::surfaces( int v )
int
Hexahedron::surfaces( int **ft, int v )
{
return make_surf_nodelist( hex_num_surf, int(hex_surface_table[1]-hex_surface_table[0]),
(const int**)hex_surface_table, v );
return make_surf_nodelist( v, ft, hex_num_surf,
int(hex_surface_table[1]-hex_surface_table[0]),
(const int**)hex_surface_table );
}
......@@ -3,7 +3,7 @@ HOSTMACHINE := $(shell uname)
FLTK_INC := $(shell fltk-config --use-gl --cxxflags)
FLTK_LD_FLAGS := $(shell fltk-config --use-gl --ldflags)
COMMON_LIBS = -lpng -lpthread -lm -lz
COMMON_LIBS = -lredblack -lpng -lpthread -lm -lz
COMMON_INC = -I. -O0 -g -DOBJ_CLASS -D_REENTRANT
ifeq ($(HOSTMACHINE),Linux)
......
......@@ -2,25 +2,50 @@
#include <set>
#include <map>
#include <vector>
#include<algorithm>
#include <string>
#include "DataOpacity.h"
#include "gl2ps.h"
#include <redblack.h>
#include "VecData.h"
struct eleFaceCmp
struct Face {
int nsort[MAX_NUM_SURF_NODES]; //!< sorted nodes
int norig[MAX_NUM_SURF_NODES]; //!< nodes in original order
int nnode; //!< number of nodes
};
int intcmp( const void *a, const void *b ){ return *(int *)a-*(int *)b; }
/** compare 2 faces for sorting */
int cmpface( const void *a, const void *b, const void *c )
{
bool operator()(const vector<int>&a, const vector<int>&b)
{
if( a.size()<b.size() ) return true;
if( a.size()>b.size() ) return false;
for( int i=0; i<a.size()-2; i++ ) // -2 since element \# and face are tacked on
if( a[i] < b[i] ) return true;
return false;
}
const Face *A = (const Face *)a;
const Face *B = (const Face *)b;
if( A->nnode != B->nnode ) return A->nnode - B->nnode;
for( int i=0; i< A->nnode; i++ )
if( A->nsort[i] != B->nsort[i] ) return A->nsort[i] - B->nsort[i];
return 0;
};
/** make a face from a node list
*
* \param f new face
* \param n number of nodes in surface face
* \param orig nodes in original order
*/
void
make_face( Face *f, int n, int* orig )
{
f->nnode = n;
memcpy( f->norig, orig, n*sizeof(int) );
memcpy( f->nsort, orig, n*sizeof(int) );
qsort( f->nsort, n, sizeof(int), intcmp );
}
/** return a new unique region label
*/
int Model::new_region_label()
......@@ -290,71 +315,76 @@ int Model::add_surface_from_elem( const char *fn )
/** find bounding surfaces for all the regions
*
* We do this by adding the faces of all elements into a set and removing a face if
* We do this by adding the faces of all elements into a tree and removing a face if
* it appears twice
*/
int Model::add_region_surfaces()
{
// copy all faces from elements in the reions into a set
// eliminating any face which appears twice
for( int r=0; r<_numReg; r++ ) {
RRegion *reg = _region[r];
set<vector<int>,eleFaceCmp> faces;
for( int e=0; e<_numVol; e++ ){
if( reg->ele_member(e) ) {
vector<vector<int> >*fnl = _vol[e]->surfaces();
for( int i=0; i<fnl->size(); i++ ) {
// sort the nodes but add the info so we can get the face back later
sort( fnl->at(i).begin(), fnl->at(i).end() );
fnl->at(i).push_back(e); // element
fnl->at(i).push_back(i); // face within element
set<vector<int>,eleFaceCmp>::iterator it;
it = faces.find( (*fnl)[i] );
if( it != faces.end() )
faces.erase(it);
else
faces.insert( (*fnl)[i] );
// set up matrix to hold element faces
int **faces = new int *[MAX_NUM_SURF];
for( int i=0; i<MAX_NUM_SURF; i++ ) faces[i] = new int[MAX_NUM_SURF_NODES+1];
for( int r=0; r<_numReg; r++ ) {
struct rbtree *facetree = rbinit( cmpface, NULL ); // tree to hold faces
Face *newface = new Face;
for( int e=0; e<_numVol; e++ )
if( _region[r]->ele_member(e) ) {
int ns = _vol[e]->surfaces( faces );
for( int i=0; i<ns; i++ ) {
Face *f1;
make_face( newface, faces[i][0], faces[i]+1 );
if ( (f1=(Face *)rbsearch( newface, facetree )) != newface ) { //in list
rbdelete( f1, facetree );
free( f1 );
} else
newface = new Face; // added to list, get memory for a new face
}
delete fnl;
}
}
// convert the left over faces into a surface
if( faces.size() ) {
// count the number of faces left in the list
delete newface;
RBLIST* rbl = rbopenlist( facetree );
int numface=0;
while( rbreadlist( rbl ) != NULL ) numface ++;
rbcloselist( rbl );
if( numface ) { // convert the left over faces into a surface
rbl = rbopenlist( facetree );
_surface.push_back( new Surfaces( &pt ) );
_surface.back()->num( faces.size() );
set<vector<int>,eleFaceCmp>::iterator fit = faces.begin();
_surface.back()->num( numface );
for( int e=0; e<faces.size(); e++, fit++ ){
int findex = fit->back();
int ele = fit->at(fit->size()-2);
vector<vector<int> >*fnl = _vol[ele]->surfaces();
int numnodes = fnl->at(findex).size();
int nl[numnodes];
for( int e=0; e<numface; e++ ){
for( int i=0; i<numnodes; i++ ) nl[i] = fnl->at(findex)[i];
newface = (Face *)rbreadlist( rbl );
if( numnodes==3 ) {
if( newface->nnode == 3 ) {
_surface.back()->ele(e) = new Triangle( &pt );
} else if( numnodes==4 ) {
} else if( newface->nnode==4 ) {
_surface.back()->ele(e) = new Quadrilateral( &pt );
} else
} else
assert(0);
_surface.back()->ele(e)->define( nl );
_surface.back()->ele(e)->define( newface->norig );
_surface.back()->ele(e)->compute_normals(0,0);
delete fnl;
delete newface; // mo longer needed
}
_surface.back()->determine_vert_norms( pt );
rbcloselist( rbl );
rbdestroy( facetree );
}
}
for( int i=0; i<MAX_NUM_SURF; i++ ) delete[] faces[i];
delete[] faces;
}
/** add a surface by reading in a .tri file, also try reading a normal file
*
/** add a surface by reading in a .tri file, also try reading a normal file
*
* \param fn file containing tri's
*
* \return the \#surfaces, -1 if not successful
......
......@@ -298,13 +298,15 @@ Prism::iso_polys(unsigned int index)
/*
* return list of lists of nodes defining the bounding surfaces
*
* \param ft face table to be filled in
* \param v element number
*
* \param pointer to a vector of vectors
*/
vector<vector<int> > *
Prism::surfaces( int v )
int
Prism::surfaces( int **ft, int v )
{
return make_surf_nodelist( prism_num_surf, int(prism_surface_table[1]-prism_surface_table[0]),
(const int **)prism_surface_table, v );
return make_surf_nodelist( v, ft, prism_num_surf,
int(prism_surface_table[1]-prism_surface_table[0]),
(const int **)prism_surface_table );
}
......@@ -245,13 +245,16 @@ const int* Pyramid::iso_polys(unsigned int index)
/*
* return list of lists of nodes defining the bounding surfaces
*
* \param v element number
* \param ft face table
* \param v element number
*
* \param pointer to a vector of vectors
*
* \pre the size of ft is big enough
*/
vector<vector<int> > *
Pyramid::surfaces( int v )
int
Pyramid::surfaces( int **ft, int v )
{
return make_surf_nodelist( pyr_num_surf, int(pyr_surface_table[1]-pyr_surface_table[0]),
(const int **)pyr_surface_table, v );
return make_surf_nodelist( v, ft, pyr_num_surf,
int(pyr_surface_table[1]-pyr_surface_table[0]), (const int **)pyr_surface_table );
}
......@@ -25,7 +25,7 @@ const int tetra_iso_table[][10] = {
};
const int tet_num_surf =4;
const int tetra_surface_table[][3] = {
{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {2, 3, 1}
{0, 1, 2}, {0, 3, 1}, {0, 2, 3}, {2, 1, 3}
};
......@@ -194,11 +194,12 @@ Tetrahedral::iso_polys(unsigned int index)
*
* \param pointer to a vector of vectors
*/
vector<vector<int> > *
Tetrahedral::surfaces( int v )
int
Tetrahedral::surfaces( int **nl, int v )
{
return make_surf_nodelist( tet_num_surf, int(tetra_surface_table[1]-tetra_surface_table[0]),
(const int **)tetra_surface_table, v );
return make_surf_nodelist( v, nl, tet_num_surf,
int(tetra_surface_table[1]-tetra_surface_table[0]),
(const int **)tetra_surface_table );
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment