stl writer for dcel
This commit is contained in:
JamesEMcclure
parent
7423610f6d
commit
4fd74a78a7
@ -1,7 +1,5 @@
|
|||||||
#include "analysis/dcel.h"
|
#include "analysis/dcel.h"
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
DECL::DECL(){
|
DECL::DECL(){
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -369,243 +367,37 @@ double DECL::EdgeAngle(int edge)
|
|||||||
return angle;
|
return angle;
|
||||||
}
|
}
|
||||||
|
|
||||||
void Isosurface(DoubleArray &A, const double &v)
|
void iso_surface(const Array<double>&Field, const double isovalue)
|
||||||
{
|
{
|
||||||
NULL_USE( v );
|
DECL object;
|
||||||
|
int e1,e2,e3;
|
||||||
Point P,Q;
|
FILE *TRIANGLES;
|
||||||
Point PlaceHolder;
|
TRIANGLES = fopen("isosurface.stl","w");
|
||||||
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
fprintf(TRIANGLES,"solid isosurface\n");
|
||||||
|
int Nx = Field.size(0);
|
||||||
int TriangleCount;
|
int Ny = Field.size(1);
|
||||||
int VertexCount;
|
int Nz = Field.size(2);
|
||||||
int CubeIndex;
|
for (int k=1; k<Nz-1; k++){
|
||||||
|
for (int j=1; j<Ny-1; j++){
|
||||||
Point VertexList[12];
|
for (int i=1; i<Nx-1; i++){
|
||||||
Point NewVertexList[12];
|
object.LocalIsosurface(Field,isovalue,i,j,k);
|
||||||
int LocalRemap[12];
|
for (int idx=0; idx<object.TriangleCount; idx++){
|
||||||
|
e1 = object.Face(idx);
|
||||||
Point cellvertices[20];
|
e2 = object.halfedge.next(e1);
|
||||||
std::array<std::array<int,3>,20> Triangles;
|
e3 = object.halfedge.next(e2);
|
||||||
Triangles.fill( { 0 } );
|
auto P1 = object.vertex.coords(object.halfedge.v1(e1));
|
||||||
|
auto P2 = object.vertex.coords(object.halfedge.v1(e2));
|
||||||
// Values from array 'A' at the cube corners
|
auto P3 = object.vertex.coords(object.halfedge.v1(e3));
|
||||||
double CubeValues[8];
|
P1.x += 1.0*i; P1.y += 1.0*j; P1.z +=1.0*k;
|
||||||
|
P2.x += 1.0*i; P2.y += 1.0*j; P2.z +=1.0*k;
|
||||||
int Nx = A.size(0);
|
P3.x += 1.0*i; P3.y += 1.0*j; P3.z +=1.0*k;
|
||||||
int Ny = A.size(1);
|
fprintf(TRIANGLES,"vertex %f %f %f\n",P1.x,P1.y,P1.z);
|
||||||
int Nz = A.size(2);
|
fprintf(TRIANGLES,"vertex %f %f %f\n",P2.x,P2.y,P2.z);
|
||||||
|
fprintf(TRIANGLES,"vertex %f %f %f\n",P3.x,P3.y,P3.z);
|
||||||
// Points corresponding to cube corners
|
|
||||||
C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
|
|
||||||
C1.x = 1.0; C1.y = 0.0; C1.z = 0.0;
|
|
||||||
C2.x = 1.0; C2.y = 1.0; C2.z = 0.0;
|
|
||||||
C3.x = 0.0; C3.y = 1.0; C3.z = 0.0;
|
|
||||||
C4.x = 0.0; C4.y = 0.0; C4.z = 1.0;
|
|
||||||
C5.x = 1.0; C5.y = 0.0; C5.z = 1.0;
|
|
||||||
C6.x = 1.0; C6.y = 1.0; C6.z = 1.0;
|
|
||||||
C7.x = 0.0; C7.y = 1.0; C7.z = 1.0;
|
|
||||||
|
|
||||||
std::vector<std::array<int,6>> HalfEdge;
|
|
||||||
for (int k=1; k<Nz-1; k++){
|
|
||||||
for (int j=1; j<Ny-1; j++){
|
|
||||||
for (int i=1; i<Nx-1; i++){
|
|
||||||
// Set the corner values for this cube
|
|
||||||
CubeValues[0] = A(i,j,k);
|
|
||||||
CubeValues[1] = A(i+1,j,k);
|
|
||||||
CubeValues[2] = A(i+1,j+1,k);
|
|
||||||
CubeValues[3] = A(i,j+1,k);
|
|
||||||
CubeValues[4] = A(i,j,k+1);
|
|
||||||
CubeValues[5] = A(i+1,j,k+1);
|
|
||||||
CubeValues[6] = A(i+1,j+1,k+1);
|
|
||||||
CubeValues[7] = A(i,j+1,k+1);
|
|
||||||
|
|
||||||
//Determine the index into the edge table which
|
|
||||||
//tells us which vertices are inside of the surface
|
|
||||||
CubeIndex = 0;
|
|
||||||
if (CubeValues[0] < 0.0f) CubeIndex |= 1;
|
|
||||||
if (CubeValues[1] < 0.0f) CubeIndex |= 2;
|
|
||||||
if (CubeValues[2] < 0.0f) CubeIndex |= 4;
|
|
||||||
if (CubeValues[3] < 0.0f) CubeIndex |= 8;
|
|
||||||
if (CubeValues[4] < 0.0f) CubeIndex |= 16;
|
|
||||||
if (CubeValues[5] < 0.0f) CubeIndex |= 32;
|
|
||||||
if (CubeValues[6] < 0.0f) CubeIndex |= 64;
|
|
||||||
if (CubeValues[7] < 0.0f) CubeIndex |= 128;
|
|
||||||
|
|
||||||
//Find the vertices where the surface intersects the cube
|
|
||||||
if (edgeTable[CubeIndex] & 1){
|
|
||||||
P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]);
|
|
||||||
VertexList[0] = P;
|
|
||||||
Q = C0;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 2){
|
|
||||||
P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]);
|
|
||||||
VertexList[1] = P;
|
|
||||||
Q = C1;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 4){
|
|
||||||
P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]);
|
|
||||||
VertexList[2] = P;
|
|
||||||
Q = C2;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 8){
|
|
||||||
P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]);
|
|
||||||
VertexList[3] = P;
|
|
||||||
Q = C3;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 16){
|
|
||||||
P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]);
|
|
||||||
VertexList[4] = P;
|
|
||||||
Q = C4;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 32){
|
|
||||||
P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]);
|
|
||||||
VertexList[5] = P;
|
|
||||||
Q = C5;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 64){
|
|
||||||
P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]);
|
|
||||||
VertexList[6] = P;
|
|
||||||
Q = C6;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 128){
|
|
||||||
P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]);
|
|
||||||
VertexList[7] = P;
|
|
||||||
Q = C7;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 256){
|
|
||||||
P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]);
|
|
||||||
VertexList[8] = P;
|
|
||||||
Q = C0;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 512){
|
|
||||||
P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]);
|
|
||||||
VertexList[9] = P;
|
|
||||||
Q = C1;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 1024){
|
|
||||||
P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]);
|
|
||||||
VertexList[10] = P;
|
|
||||||
Q = C2;
|
|
||||||
}
|
|
||||||
if (edgeTable[CubeIndex] & 2048){
|
|
||||||
P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]);
|
|
||||||
VertexList[11] = P;
|
|
||||||
Q = C3;
|
|
||||||
}
|
|
||||||
|
|
||||||
VertexCount=0;
|
|
||||||
for (int idx=0;idx<12;idx++)
|
|
||||||
LocalRemap[idx] = -1;
|
|
||||||
|
|
||||||
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++)
|
|
||||||
{
|
|
||||||
if(LocalRemap[triTable[CubeIndex][idx]] == -1)
|
|
||||||
{
|
|
||||||
NewVertexList[VertexCount] = VertexList[triTable[CubeIndex][idx]];
|
|
||||||
LocalRemap[triTable[CubeIndex][idx]] = VertexCount;
|
|
||||||
VertexCount++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
for (int idx=0;idx<VertexCount;idx++) {
|
|
||||||
P = NewVertexList[idx];
|
|
||||||
//P.x += i;
|
|
||||||
//P.y += j;
|
|
||||||
//P.z += k;
|
|
||||||
cellvertices[idx] = P;
|
|
||||||
}
|
|
||||||
|
|
||||||
TriangleCount = 0;
|
|
||||||
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
|
|
||||||
Triangles[TriangleCount][0] = LocalRemap[triTable[CubeIndex][idx+0]];
|
|
||||||
Triangles[TriangleCount][1] = LocalRemap[triTable[CubeIndex][idx+1]];
|
|
||||||
Triangles[TriangleCount][2] = LocalRemap[triTable[CubeIndex][idx+2]];
|
|
||||||
TriangleCount++;
|
|
||||||
}
|
|
||||||
int nTris = TriangleCount;
|
|
||||||
|
|
||||||
// Now add the local values to the DECL data structure
|
|
||||||
HalfEdge.resize(nTris*3);
|
|
||||||
int idx_edge=0;
|
|
||||||
for (int idx=0; idx<TriangleCount; idx++){
|
|
||||||
int V1 = Triangles[idx][0];
|
|
||||||
int V2 = Triangles[idx][1];
|
|
||||||
int V3 = Triangles[idx][2];
|
|
||||||
// first edge: V1->V2
|
|
||||||
HalfEdge[idx_edge][0] = V1; // first vertex
|
|
||||||
HalfEdge[idx_edge][1] = V2; // second vertex
|
|
||||||
HalfEdge[idx_edge][2] = idx; // triangle
|
|
||||||
HalfEdge[idx_edge][3] = -1; // twin
|
|
||||||
HalfEdge[idx_edge][4] = idx_edge+2; // previous edge
|
|
||||||
HalfEdge[idx_edge][5] = idx_edge+1; // next edge
|
|
||||||
idx_edge++;
|
|
||||||
// second edge: V2->V3
|
|
||||||
HalfEdge[idx_edge][0] = V2; // first vertex
|
|
||||||
HalfEdge[idx_edge][1] = V3; // second vertex
|
|
||||||
HalfEdge[idx_edge][2] = idx; // triangle
|
|
||||||
HalfEdge[idx_edge][3] = -1; // twin
|
|
||||||
HalfEdge[idx_edge][4] = idx_edge-1; // previous edge
|
|
||||||
HalfEdge[idx_edge][5] = idx_edge+1; // next edge
|
|
||||||
idx_edge++;
|
|
||||||
// third edge: V3->V1
|
|
||||||
HalfEdge[idx_edge][0] = V3; // first vertex
|
|
||||||
HalfEdge[idx_edge][1] = V1; // second vertex
|
|
||||||
HalfEdge[idx_edge][2] = idx; // triangle
|
|
||||||
HalfEdge[idx_edge][3] = -1; // twin
|
|
||||||
HalfEdge[idx_edge][4] = idx_edge-1; // previous edge
|
|
||||||
HalfEdge[idx_edge][5] = idx_edge-2; // next edge
|
|
||||||
idx_edge++;
|
|
||||||
}
|
|
||||||
int EdgeCount=idx_edge;
|
|
||||||
for (int idx=0; idx<EdgeCount; idx++){
|
|
||||||
int V1=HalfEdge[idx][0];
|
|
||||||
int V2=HalfEdge[idx][1];
|
|
||||||
// Find all the twins within the cube
|
|
||||||
for (int jdx=0; idx<EdgeCount; jdx++){
|
|
||||||
if (HalfEdge[jdx][1] == V1 && HalfEdge[jdx][0] == V2){
|
|
||||||
// this is the pair
|
|
||||||
HalfEdge[idx][3] = jdx;
|
|
||||||
HalfEdge[jdx][3] = idx;
|
|
||||||
}
|
|
||||||
if (HalfEdge[jdx][1] == V2 && HalfEdge[jdx][0] == V1 && !(idx==jdx)){
|
|
||||||
std::printf("WARNING: half edges with identical orientation! \n");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// Use "ghost" twins if edge is on a cube face
|
|
||||||
P = cellvertices[V1];
|
|
||||||
Q = cellvertices[V2];
|
|
||||||
if (P.x == 0.0 && Q.x == 0.0) HalfEdge[idx_edge][3] = -1; // ghost twin for x=0 face
|
|
||||||
if (P.x == 1.0 && Q.x == 1.0) HalfEdge[idx_edge][3] = -2; // ghost twin for x=1 face
|
|
||||||
if (P.y == 0.0 && Q.y == 0.0) HalfEdge[idx_edge][3] = -3; // ghost twin for y=0 face
|
|
||||||
if (P.y == 1.0 && Q.y == 1.0) HalfEdge[idx_edge][3] = -4; // ghost twin for y=1 face
|
|
||||||
if (P.z == 0.0 && Q.z == 0.0) HalfEdge[idx_edge][3] = -5; // ghost twin for z=0 face
|
|
||||||
if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6; // ghost twin for z=1 face
|
|
||||||
}
|
|
||||||
// Find all the angles
|
|
||||||
/*for (int idx=0; idx<EdgeCount; idx++){
|
|
||||||
int V1=HalfEdge[idx][0];
|
|
||||||
int V2=HalfEdge[idx][1];
|
|
||||||
int T1= HalfEdge[idx_edge][2];
|
|
||||||
int twin=HalfEdge[idx_edge][3];
|
|
||||||
if (twin == -1){
|
|
||||||
|
|
||||||
}
|
|
||||||
}*/
|
|
||||||
|
|
||||||
// Map vertices to global coordinates
|
|
||||||
for (int idx=0;idx<VertexCount;idx++) {
|
|
||||||
P = cellvertices[idx];
|
|
||||||
P.x += i;
|
|
||||||
P.y += j;
|
|
||||||
P.z += k;
|
|
||||||
cellvertices[idx] = P;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fclose(TRIANGLES);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -1,3 +1,6 @@
|
|||||||
|
#ifndef DCEL_INC
|
||||||
|
#define DCEL_INC
|
||||||
|
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include "analysis/pmmc.h"
|
#include "analysis/pmmc.h"
|
||||||
|
|
||||||
@ -79,3 +82,7 @@ public:
|
|||||||
private:
|
private:
|
||||||
std::vector<int> FaceData;
|
std::vector<int> FaceData;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
void iso_surface(const Array<double>&Field, const double isovalue);
|
||||||
|
|
||||||
|
#endif
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user