OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2018-09-28 22:38:40 -04:00
commit 3acf143305
6 changed files with 141 additions and 325 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
analysis/Minkowski.cpp
View File

@ -32,13 +32,6 @@
#define PI 3.14159265359
// Constructor
Minkowski::Minkowski():
kstart(0), kfinish(0), isovalue(0), Volume(0),
LOGFILE(NULL), Dm(NULL), Vi(0), Vi_global(0)
{
}
Minkowski::Minkowski(std::shared_ptr <Domain> dm):
kstart(0), kfinish(0), isovalue(0), Volume(0),
LOGFILE(NULL), Dm(dm), Vi(0), Vi_global(0)
@ -61,7 +54,7 @@ Minkowski::Minkowski(std::shared_ptr <Domain> dm):
// Destructor
Minkowski::~Minkowski()
{
// if ( LOGFILE!=NULL ) { fclose(LOGFILE); }
if ( LOGFILE!=NULL ) { fclose(LOGFILE); }
}
double Minkowski::V(){

4
analysis/Minkowski.h
View File

@ -60,8 +60,8 @@ public:
double J();
double X();
//...........................................................................
Minkowski();
//..........................................................................
Minkowski(){};//NULL CONSTRUCTOR
Minkowski(std::shared_ptr <Domain> Dm);
~Minkowski();
void ComputeScalar(const DoubleArray& Field, const double isovalue);

158
analysis/TwoPhase.cpp
View File

@ -103,6 +103,9 @@ TwoPhase::TwoPhase(std::shared_ptr <Domain> dm):
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0;
TempID = new char[Nx*Ny*Nz];
wet_morph = std::shared_ptr<Minkowski>(new Minkowski(Dm));
nonwet_morph = std::shared_ptr<Minkowski>(new Minkowski(Dm));
// Global arrays
PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0);
@ -180,7 +183,9 @@ TwoPhase::TwoPhase(std::shared_ptr <Domain> dm):
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn "); //trimmed curvature,
fprintf(TIMELOG,"wwndnw wwnsdnwn Jwnwwndnw "); //kinematic quantities,
fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
fprintf(TIMELOG,"Vw Aw Jw Xw "); //miknowski measures,
fprintf(TIMELOG,"Vn An Jn Xn\n"); //miknowski measures,
// fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
}
NWPLOG = fopen("components.NWP.tcat","a+");
@ -209,7 +214,9 @@ TwoPhase::TwoPhase(std::shared_ptr <Domain> dm):
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn "); //trimmed curvature,
fprintf(TIMELOG,"wwndnw wwnsdnwn Jwnwwndnw "); //kinematic quantities,
fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
fprintf(TIMELOG,"Vw Aw Jw Xw "); //miknowski measures,
fprintf(TIMELOG,"Vn An Jn Xn\n"); //miknowski measures,
// fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
}
}
@ -572,7 +579,7 @@ void TwoPhase::ComputeLocal()
}
}
}
/*
Array <char> phase_label(Nx,Ny,Nz);
Array <double> phase_distance(Nx,Ny,Nz);
// Analyze the wetting fluid
@ -597,8 +604,8 @@ void TwoPhase::ComputeLocal()
}
}
CalcDist(phase_distance,phase_label,*Dm);
wet_morph.ComputeScalar(phase_distance,0.f);
printf("generating distance at rank=%i \n",Dm->rank());
wet_morph->ComputeScalar(phase_distance,0.f);
//printf("generating distance at rank=%i \n",Dm->rank());
// Analyze the wetting fluid
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
@ -620,12 +627,12 @@ void TwoPhase::ComputeLocal()
}
}
}
printf("calculate distance at rank=%i \n",Dm->rank());
//printf("calculate distance at rank=%i \n",Dm->rank());
CalcDist(phase_distance,phase_label,*Dm);
printf("morphological analysis at rank=%i \n",Dm->rank());
nonwet_morph.ComputeScalar(phase_distance,0.f);
printf("rank=%i completed \n",Dm->rank());
*/
//printf("morphological analysis at rank=%i \n",Dm->rank());
nonwet_morph->ComputeScalar(phase_distance,0.f);
//printf("rank=%i completed \n",Dm->rank());
}
@ -1076,129 +1083,6 @@ void TwoPhase::ComponentAverages()
}
}
void TwoPhase::WriteSurfaces(int logcount)
{
int i,j,k;
int ncubes=(Nx-1)*(Ny-1)*(Nz-1);
Point P,A,B,C;
std::shared_ptr<IO::TriList> wn_mesh( new IO::TriList() );
wn_mesh->A.reserve(8*ncubes);
wn_mesh->B.reserve(8*ncubes);
wn_mesh->C.reserve(8*ncubes);
std::shared_ptr<IO::TriList> ns_mesh( new IO::TriList() );
ns_mesh->A.reserve(8*ncubes);
ns_mesh->B.reserve(8*ncubes);
ns_mesh->C.reserve(8*ncubes);
std::shared_ptr<IO::TriList> ws_mesh( new IO::TriList() );
ws_mesh->A.reserve(8*ncubes);
ws_mesh->B.reserve(8*ncubes);
ws_mesh->C.reserve(8*ncubes);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){ // Get cube from the list
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris,
local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
i, j, k, Nx, Ny, Nz);
//.......................................................................................
// Write the triangle lists to text file
for (int r=0;r<n_nw_tris;r++){
A = nw_pts(nw_tris(0,r));
B = nw_pts(nw_tris(1,r));
C = nw_pts(nw_tris(2,r));
// compare the trianlge orientation against the color gradient
// Orientation of the triangle
double tri_normal_x = (A.y-B.y)*(B.z-C.z) - (A.z-B.z)*(B.y-C.y);
double tri_normal_y = (A.z-B.z)*(B.x-C.x) - (A.x-B.x)*(B.z-C.z);
double tri_normal_z = (A.x-B.x)*(B.y-C.y) - (A.y-B.y)*(B.x-C.x);
double normal_x = SDn_x(i,j,k);
double normal_y = SDn_y(i,j,k);
double normal_z = SDn_z(i,j,k);
// If the normals don't point in the same direction, flip the orientation of the triangle
// Right hand rule for triangle orientation is used to determine rendering for most software
if (normal_x*tri_normal_x + normal_y*tri_normal_y + normal_z*tri_normal_z < 0.0){
P = A;
A = C;
C = P;
}
// Remap the points
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
wn_mesh->A.push_back(A);
wn_mesh->B.push_back(B);
wn_mesh->C.push_back(C);
}
for (int r=0;r<n_ws_tris;r++){
A = ws_pts(ws_tris(0,r));
B = ws_pts(ws_tris(1,r));
C = ws_pts(ws_tris(2,r));
// Remap the points
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
ws_mesh->A.push_back(A);
ws_mesh->B.push_back(B);
ws_mesh->C.push_back(C);
}
for (int r=0;r<n_ns_tris;r++){
A = ns_pts(ns_tris(0,r));
B = ns_pts(ns_tris(1,r));
C = ns_pts(ns_tris(2,r));
// Remap the points
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
ns_mesh->A.push_back(A);
ns_mesh->B.push_back(B);
ns_mesh->C.push_back(C);
}
}
}
}
std::vector<IO::MeshDataStruct> meshData(4);
meshData[0].meshName = "wn-tris";
meshData[0].mesh = wn_mesh;
meshData[1].meshName = "ws-tris";
meshData[1].mesh = ws_mesh;
meshData[2].meshName = "ns-tris";
meshData[2].mesh = ns_mesh;
IO::writeData( logcount, meshData, Dm->Comm );
}
void TwoPhase::Reduce()
{
int i;
@ -1327,7 +1211,9 @@ void TwoPhase::PrintAll(int timestep)
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g ",wwndnw_global, wwnsdnwn_global, Jwnwwndnw_global); // kinematic quantities
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",wet_morph->V(), wet_morph->A(), wet_morph->J(), wet_morph->X());
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",nonwet_morph->V(), nonwet_morph->A(), nonwet_morph->J(), nonwet_morph->X());
// fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fflush(TIMELOG);
}
else{
@ -1352,7 +1238,9 @@ void TwoPhase::PrintAll(int timestep)
Gws(0),Gws(1),Gws(2),Gws(3),Gws(4),Gws(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn, trJwn, trRwn); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g ",wwndnw, wwnsdnwn, Jwnwwndnw); // kinematic quantities
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler, Kn, Jn, An); // minkowski measures
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",wet_morph->Vi, wet_morph->Ai, wet_morph->Ji, wet_morph->Xi);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",nonwet_morph->Vi, nonwet_morph->Ai, nonwet_morph->Ji, nonwet_morph->Xi);
// fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler, Kn, Jn, An); // minkowski measures
fflush(TIMELOG);
}

5
analysis/TwoPhase.h
View File

@ -160,8 +160,8 @@ public:
DoubleArray Vel_z;
DoubleArray PhaseDistance;
Minkowski wet_morph;
Minkowski nonwet_morph;
std::shared_ptr<Minkowski> wet_morph;
std::shared_ptr<Minkowski> nonwet_morph;
// Container for averages;
DoubleArray ComponentAverages_WP;
@ -179,7 +179,6 @@ public:
void AssignComponentLabels();
void ComponentAverages();
void Reduce();
void WriteSurfaces(int logcount);
void NonDimensionalize(double D, double viscosity, double IFT);
void PrintAll(int timestep);
int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);

229
analysis/dcel.cpp
View File

@ -1,78 +1,5 @@
#include "analysis/dcel.h"
/*
Double connected edge list (DECL)
*/
Vertex::Vertex(){
size_ = 0;
vertex_data.resize(36);
}
Vertex::~Vertex(){
}
void Vertex::add(Point P){
vertex_data.push_back(P.x);
vertex_data.push_back(P.y);
vertex_data.push_back(P.z);
size_++;
}
void Vertex::assign(int idx, Point P){
vertex_data[3*idx] = P.x;
vertex_data[3*idx+1] = P.y;
vertex_data[3*idx+2] = P.z;
}
int Vertex::size(){
return size_;
}
Point Vertex::coords(int idx){
Point P;
P.x = vertex_data[3*idx];
P.y = vertex_data[3*idx+1];
P.z = vertex_data[3*idx+2];
return P;
}
Halfedge::Halfedge(){
size_=0;
}
Halfedge::~Halfedge(){
}
int Halfedge::v1(int edge){
return data(0,edge);
}
int Halfedge::v2(int edge){
return data(1,edge);
}
int Halfedge::face(int edge){
return data(2,edge);
}
int Halfedge::twin(int edge){
return data(3,edge);
}
int Halfedge::prev(int edge){
return data(4,edge);
}
int Halfedge::next(int edge){
return data(5,edge);
}
int Halfedge::size(){
return size_;
}
DECL::DECL(){
@ -85,24 +12,20 @@ DECL::~DECL(){
}
int DECL::Face(int index){
return FaceData(index);
return FaceData[index];
}
void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const int j, const int k){
void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, const int j, const int k){
Point P,Q;
Point PlaceHolder;
Point C0,C1,C2,C3,C4,C5,C6,C7;
int CubeIndex;
int nTris = 0;
int nVert =0;
Point VertexList[12];
Point NewVertexList[12];
int LocalRemap[12];
DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
IntArray Triangles = IntArray(3,20);
Point cellvertices[20];
std::array<std::array<int,3>,20> Triangles;
// Values from array 'A' at the cube corners
double CubeValues[8];
@ -129,7 +52,7 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
//Determine the index into the edge table which
//tells us which vertices are inside of the surface
CubeIndex = 0;
int CubeIndex = 0;
if (CubeValues[0] < 0.0f) CubeIndex |= 1;
if (CubeValues[1] < 0.0f) CubeIndex |= 2;
if (CubeValues[2] < 0.0f) CubeIndex |= 4;
@ -222,31 +145,30 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
//P.x += i;
//P.y += j;
//P.z += k;
cellvertices(idx) = P;
cellvertices[idx] = P;
}
nVert = VertexCount;
TriangleCount = 0;
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
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++;
}
nTris = TriangleCount;
int nTris = TriangleCount;
// Now add the local values to the DECL data structure
if (nTris>0){
FaceData.resize(TriangleCount);
//printf("Construct halfedge structure... \n");
//printf(" Construct %i triangles \n",nTris);
halfedge.data.resize(6,nTris*3);
halfedge.resize(nTris*3);
int idx_edge=0;
for (int idx=0; idx<TriangleCount; idx++){
int V1 = Triangles(0,idx);
int V2 = Triangles(1,idx);
int V3 = Triangles(2,idx);
FaceData(idx) = idx_edge;
int V1 = Triangles[idx][0];
int V2 = Triangles[idx][1];
int V3 = Triangles[idx][2];
FaceData[idx] = idx_edge;
// first edge: V1->V2
halfedge.data(0,idx_edge) = V1; // first vertex
halfedge.data(1,idx_edge) = V2; // second vertex
@ -290,8 +212,8 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
}
}
// Use "ghost" twins if edge is on a cube face
P = cellvertices(V1);
Q = cellvertices(V2);
P = cellvertices[V1];
Q = cellvertices[V2];
if (P.x == 0.0 && Q.x == 0.0) halfedge.data(3,idx) = -1; // ghost twin for x=0 face
if (P.x == 1.0 && Q.x == 1.0) halfedge.data(3,idx) = -4; // ghost twin for x=1 face
if (P.y == 0.0 && Q.y == 0.0) halfedge.data(3,idx) = -2; // ghost twin for y=0 face
@ -303,7 +225,7 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, const int i, const
// Map vertices to global coordinates
for (int idx=0;idx<VertexCount;idx++) {
P = cellvertices(idx);
P = cellvertices[idx];
P.x += i;
P.y += j;
P.z += k;
@ -358,8 +280,7 @@ Point DECL::TriNormal(int edge)
double DECL::EdgeAngle(int edge)
{
double angle;
double nx,ny,nz;
double dotprod,length,hypotenuse;
double dotprod;
Point P,Q,R; // triangle vertices
Point U,V,W; // normal vectors
int e2 = halfedge.next(edge);
@ -372,14 +293,14 @@ double DECL::EdgeAngle(int edge)
if (halfedge.twin(edge) < 0 ){
// compute edge normal in plane of cube face
W = P - Q; // edge tangent vector
length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z);
double length = sqrt(W.x*W.x+W.y*W.y+W.z*W.z);
W.x /= length;
W.y /= length;
W.z /= length;
// edge normal within the plane of the cube face
nx = W.y*V.z - W.z*V.y;
ny = W.z*V.x - W.x*V.z;
nz = W.x*V.y - W.y*V.x;
double nx = W.y*V.z - W.z*V.y;
double ny = W.z*V.x - W.x*V.z;
double nz = W.x*V.y - W.y*V.x;
length = sqrt(nx*nx+ny*ny+nz*nz);
// new value for V is this normal vector
V.x = nx/length; V.y = ny/length; V.z = nz/length;
@ -431,20 +352,19 @@ void Isosurface(DoubleArray &A, const double &v)
{
Point P,Q;
Point PlaceHolder;
double temp;
Point C0,C1,C2,C3,C4,C5,C6,C7;
int TriangleCount;
int VertexCount;
int CubeIndex;
int nTris, nVert;
Point VertexList[12];
Point NewVertexList[12];
int LocalRemap[12];
DTMutableList<Point> cellvertices = DTMutableList<Point>(20);
IntArray Triangles = IntArray(3,20);
Point cellvertices[20];
std::array<std::array<int,3>,20> Triangles;
Triangles.fill( { 0 } );
// Values from array 'A' at the cube corners
double CubeValues[8];
@ -463,6 +383,7 @@ void Isosurface(DoubleArray &A, const double &v)
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++){
@ -569,83 +490,81 @@ void Isosurface(DoubleArray &A, const double &v)
//P.x += i;
//P.y += j;
//P.z += k;
cellvertices(idx) = P;
cellvertices[idx] = P;
}
nVert = VertexCount;
TriangleCount = 0;
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
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++;
}
nTris = TriangleCount;
int nTris = TriangleCount;
// Now add the local values to the DECL data structure
IntArray HalfEdge(6,nTris*3);
DoubleArray EdgeAngles(nTris*3);
HalfEdge.resize(nTris*3);
int idx_edge=0;
for (int idx=0; idx<TriangleCount; idx++){
int V1 = Triangles(0,idx);
int V2 = Triangles(1,idx);
int V3 = Triangles(2,idx);
int V1 = Triangles[idx][0];
int V2 = Triangles[idx][1];
int V3 = Triangles[idx][2];
// first edge: V1->V2
HalfEdge(0,idx_edge) = V1; // first vertex
HalfEdge(1,idx_edge) = V2; // second vertex
HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge(3,idx_edge) = -1; // twin
HalfEdge(4,idx_edge) = idx_edge+2; // previous edge
HalfEdge(5,idx_edge) = idx_edge+1; // next edge
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(0,idx_edge) = V2; // first vertex
HalfEdge(1,idx_edge) = V3; // second vertex
HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge(3,idx_edge) = -1; // twin
HalfEdge(4,idx_edge) = idx_edge-1; // previous edge
HalfEdge(5,idx_edge) = idx_edge+1; // next edge
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(0,idx_edge) = V3; // first vertex
HalfEdge(1,idx_edge) = V1; // second vertex
HalfEdge(2,idx_edge) = idx; // triangle
HalfEdge(3,idx_edge) = -1; // twin
HalfEdge(4,idx_edge) = idx_edge-1; // previous edge
HalfEdge(5,idx_edge) = idx_edge-2; // next edge
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(0,idx);
int V2=HalfEdge(1,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(1,jdx) == V1 && HalfEdge(0,jdx) == V2){
if (HalfEdge[jdx][1] == V1 && HalfEdge[jdx][0] == V2){
// this is the pair
HalfEdge(3,idx) = jdx;
HalfEdge(3,jdx) = idx;
HalfEdge[idx][3] = jdx;
HalfEdge[jdx][3] = idx;
}
if (HalfEdge(1,jdx) == V2 && HalfEdge(0,jdx) == V1 && !(idx==jdx)){
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(3,idx_edge) = -1; // ghost twin for x=0 face
if (P.x == 1.0 && Q.x == 1.0) HalfEdge(3,idx_edge) = -2; // ghost twin for x=1 face
if (P.y == 0.0 && Q.y == 0.0) HalfEdge(3,idx_edge) = -3; // ghost twin for y=0 face
if (P.y == 1.0 && Q.y == 1.0) HalfEdge(3,idx_edge) = -4; // ghost twin for y=1 face
if (P.z == 0.0 && Q.z == 0.0) HalfEdge(3,idx_edge) = -5; // ghost twin for z=0 face
if (P.z == 1.0 && Q.z == 1.0) HalfEdge(3,idx_edge) = -6; // ghost twin for z=1 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(0,idx);
int V2=HalfEdge(1,idx);
int T1= HalfEdge(2,idx_edge);
int twin=HalfEdge(3,idx_edge);
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){
}
@ -653,11 +572,11 @@ void Isosurface(DoubleArray &A, const double &v)
// Map vertices to global coordinates
for (int idx=0;idx<VertexCount;idx++) {
P = cellvertices(idx);
P = cellvertices[idx];
P.x += i;
P.y += j;
P.z += k;
cellvertices(idx) = P;
cellvertices[idx] = P;
}
}
}

61
analysis/dcel.h
View File

@ -8,36 +8,53 @@ Doubly-connected edge list (DECL)
// Vertex structure
class Vertex{
public:
Vertex();
~Vertex();
void add(Point P);
void assign( int idx, Point P);
int size();
Point coords(int idx);
Vertex() { d_data.resize(12); }
~Vertex() = default;
Vertex( const Vertex& ) = delete;
Vertex operator=( const Vertex& ) = delete;
// Add/assign a point
inline void add( const Point& P ) { d_data.push_back( P ); }
inline void assign( int idx, const Point& P ) { d_data[idx] = P; }
// Get a point
inline Point& coords( int idx ) { return d_data[idx]; }
inline const Point& coords( int idx ) const { return d_data[idx]; }
int IncidentEdge();
// Return the number of points
inline int size() const { return d_data.size(); }
private:
std::vector<double> vertex_data;
int size_;
std::vector<Point> d_data;
};
// Halfedge structure
// Face
class Halfedge{
public:
Halfedge();
~Halfedge();
Halfedge() = default;
~Halfedge() = default;
Halfedge( const Halfedge& ) = delete;
Halfedge operator=( const Halfedge& ) = delete;
int v1(int edge);
int v2(int edge);
int twin(int edge);
int face(int edge);
int next(int edge);
int prev(int edge);
int size();
inline int v1(int edge) const { return d_data[edge][0]; }
inline int v2(int edge) const { return d_data[edge][1]; }
inline int face(int edge) const { return d_data[edge][2]; }
inline int twin(int edge) const { return d_data[edge][3]; }
inline int prev(int edge) const { return d_data[edge][4]; }
inline int next(int edge) const { return d_data[edge][5]; }
inline int size() const { return d_data.size(); }
inline void resize( int N ) { d_data.resize( N ); }
inline int& data( int i, int j ) { return d_data[j][i]; }
inline const int& data( int i, int j ) const { return d_data[j][i]; }
Array<int> data;
private:
int size_;
std::vector<std::array<int,6>> d_data;
};
// DECL
@ -49,15 +66,15 @@ public:
int face();
Vertex vertex;
Halfedge halfedge;
void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k);
void LocalIsosurface(const DoubleArray& A, double value, int i, int j, int k);
int Face(int index);
double origin(int edge);
double EdgeAngle(int edge);
Point TriNormal(int edge);
int TriangleCount;
int VertexCount;
int VertexCount;
private:
Array <int> FaceData;
std::vector<int> FaceData;
};