working on DECL isosurface tools
This commit is contained in:
James E McClure
parent
39d7fa0430
commit
db1fe9e327
@ -33,38 +33,51 @@ Point Vertex::coords(unsigned long int idx){
|
|||||||
return P;
|
return P;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
Halfedge::Halfedge(){
|
||||||
|
}
|
||||||
|
|
||||||
|
Halfedge::~Halfedge(){
|
||||||
|
|
||||||
|
}
|
||||||
unsigned long int Halfedge::v1(unsigned long int edge){
|
unsigned long int Halfedge::v1(unsigned long int edge){
|
||||||
return HalfEdge(0,edge);
|
return data(0,edge);
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned long int Halfedge::v2(unsigned long int edge){
|
unsigned long int Halfedge::v2(unsigned long int edge){
|
||||||
return HalfEdge(1,edge);
|
return data(1,edge);
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned long int Halfedge::face(unsigned long int edge){
|
unsigned long int Halfedge::face(unsigned long int edge){
|
||||||
return HalfEdge(2,edge);
|
return data(2,edge);
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned long int Halfedge::twin(unsigned long int edge){
|
unsigned long int Halfedge::twin(unsigned long int edge){
|
||||||
return HalfEdge(3,edge);
|
return data(3,edge);
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned long int Halfedge::prev(unsigned long int edge){
|
unsigned long int Halfedge::prev(unsigned long int edge){
|
||||||
return HalfEdge(4,edge);
|
return data(4,edge);
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned long int Halfedge::next(unsigned long int edge){
|
unsigned long int Halfedge::next(unsigned long int edge){
|
||||||
return HalfEdge(5,edge);
|
return data(5,edge);
|
||||||
|
}
|
||||||
|
|
||||||
|
DECL::DECL(){
|
||||||
|
}
|
||||||
|
|
||||||
|
DECL::~DECL(){
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int k){
|
void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int k){
|
||||||
Point P,Q;
|
Point P,Q;
|
||||||
Point PlaceHolder;
|
Point PlaceHolder;
|
||||||
double temp;
|
|
||||||
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
||||||
|
|
||||||
int TriangleCount;
|
int TriangleCount;
|
||||||
int NewVertexCount;
|
int NewVertexCount;
|
||||||
int CubeIndex;
|
int CubeIndex;
|
||||||
int nTris, nVert;
|
int nTris, nVert;
|
||||||
|
|
||||||
@ -78,10 +91,6 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int
|
|||||||
// Values from array 'A' at the cube corners
|
// Values from array 'A' at the cube corners
|
||||||
double CubeValues[8];
|
double CubeValues[8];
|
||||||
|
|
||||||
int Nx = A.size(0);
|
|
||||||
int Ny = A.size(1);
|
|
||||||
int Nz = A.size(2);
|
|
||||||
|
|
||||||
// Points corresponding to cube corners
|
// Points corresponding to cube corners
|
||||||
C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
|
C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
|
||||||
C1.x = 1.0; C1.y = 0.0; C1.z = 0.0;
|
C1.x = 1.0; C1.y = 0.0; C1.z = 0.0;
|
||||||
@ -92,188 +101,186 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int
|
|||||||
C6.x = 1.0; C6.y = 1.0; C6.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;
|
C7.x = 0.0; C7.y = 1.0; C7.z = 1.0;
|
||||||
|
|
||||||
CubeValues[0] = A(i,j,k) - value;
|
CubeValues[0] = A(i,j,k) - value;
|
||||||
CubeValues[1] = A(i+1,j,k) - value;
|
CubeValues[1] = A(i+1,j,k) - value;
|
||||||
CubeValues[2] = A(i+1,j+1,k) - value;
|
CubeValues[2] = A(i+1,j+1,k) - value;
|
||||||
CubeValues[3] = A(i,j+1,k) - value;
|
CubeValues[3] = A(i,j+1,k) - value;
|
||||||
CubeValues[4] = A(i,j,k+1) - value;
|
CubeValues[4] = A(i,j,k+1) - value;
|
||||||
CubeValues[5] = A(i+1,j,k+1) - value;
|
CubeValues[5] = A(i+1,j,k+1) - value;
|
||||||
CubeValues[6] = A(i+1,j+1,k+1) - value;
|
CubeValues[6] = A(i+1,j+1,k+1) - value;
|
||||||
CubeValues[7] = A(i,j+1,k+1) -value;
|
CubeValues[7] = A(i,j+1,k+1) -value;
|
||||||
|
|
||||||
//Determine the index into the edge table which
|
//Determine the index into the edge table which
|
||||||
//tells us which vertices are inside of the surface
|
//tells us which vertices are inside of the surface
|
||||||
CubeIndex = 0;
|
CubeIndex = 0;
|
||||||
if (CubeValues[0] < 0.0f) CubeIndex |= 1;
|
if (CubeValues[0] < 0.0f) CubeIndex |= 1;
|
||||||
if (CubeValues[1] < 0.0f) CubeIndex |= 2;
|
if (CubeValues[1] < 0.0f) CubeIndex |= 2;
|
||||||
if (CubeValues[2] < 0.0f) CubeIndex |= 4;
|
if (CubeValues[2] < 0.0f) CubeIndex |= 4;
|
||||||
if (CubeValues[3] < 0.0f) CubeIndex |= 8;
|
if (CubeValues[3] < 0.0f) CubeIndex |= 8;
|
||||||
if (CubeValues[4] < 0.0f) CubeIndex |= 16;
|
if (CubeValues[4] < 0.0f) CubeIndex |= 16;
|
||||||
if (CubeValues[5] < 0.0f) CubeIndex |= 32;
|
if (CubeValues[5] < 0.0f) CubeIndex |= 32;
|
||||||
if (CubeValues[6] < 0.0f) CubeIndex |= 64;
|
if (CubeValues[6] < 0.0f) CubeIndex |= 64;
|
||||||
if (CubeValues[7] < 0.0f) CubeIndex |= 128;
|
if (CubeValues[7] < 0.0f) CubeIndex |= 128;
|
||||||
|
|
||||||
//Find the vertices where the surface intersects the cube
|
//Find the vertices where the surface intersects the cube
|
||||||
if (edgeTable[CubeIndex] & 1){
|
if (edgeTable[CubeIndex] & 1){
|
||||||
P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]);
|
P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]);
|
||||||
VertexList[0] = P;
|
VertexList[0] = P;
|
||||||
Q = C0;
|
Q = C0;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 2){
|
if (edgeTable[CubeIndex] & 2){
|
||||||
P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]);
|
P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]);
|
||||||
VertexList[1] = P;
|
VertexList[1] = P;
|
||||||
Q = C1;
|
Q = C1;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 4){
|
if (edgeTable[CubeIndex] & 4){
|
||||||
P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]);
|
P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]);
|
||||||
VertexList[2] = P;
|
VertexList[2] = P;
|
||||||
Q = C2;
|
Q = C2;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 8){
|
if (edgeTable[CubeIndex] & 8){
|
||||||
P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]);
|
P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]);
|
||||||
VertexList[3] = P;
|
VertexList[3] = P;
|
||||||
Q = C3;
|
Q = C3;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 16){
|
if (edgeTable[CubeIndex] & 16){
|
||||||
P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]);
|
P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]);
|
||||||
VertexList[4] = P;
|
VertexList[4] = P;
|
||||||
Q = C4;
|
Q = C4;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 32){
|
if (edgeTable[CubeIndex] & 32){
|
||||||
P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]);
|
P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]);
|
||||||
VertexList[5] = P;
|
VertexList[5] = P;
|
||||||
Q = C5;
|
Q = C5;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 64){
|
if (edgeTable[CubeIndex] & 64){
|
||||||
P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]);
|
P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]);
|
||||||
VertexList[6] = P;
|
VertexList[6] = P;
|
||||||
Q = C6;
|
Q = C6;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 128){
|
if (edgeTable[CubeIndex] & 128){
|
||||||
P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]);
|
P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]);
|
||||||
VertexList[7] = P;
|
VertexList[7] = P;
|
||||||
Q = C7;
|
Q = C7;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 256){
|
if (edgeTable[CubeIndex] & 256){
|
||||||
P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]);
|
P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]);
|
||||||
VertexList[8] = P;
|
VertexList[8] = P;
|
||||||
Q = C0;
|
Q = C0;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 512){
|
if (edgeTable[CubeIndex] & 512){
|
||||||
P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]);
|
P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]);
|
||||||
VertexList[9] = P;
|
VertexList[9] = P;
|
||||||
Q = C1;
|
Q = C1;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 1024){
|
if (edgeTable[CubeIndex] & 1024){
|
||||||
P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]);
|
P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]);
|
||||||
VertexList[10] = P;
|
VertexList[10] = P;
|
||||||
Q = C2;
|
Q = C2;
|
||||||
}
|
}
|
||||||
if (edgeTable[CubeIndex] & 2048){
|
if (edgeTable[CubeIndex] & 2048){
|
||||||
P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]);
|
P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]);
|
||||||
VertexList[11] = P;
|
VertexList[11] = P;
|
||||||
Q = C3;
|
Q = C3;
|
||||||
}
|
}
|
||||||
|
|
||||||
NewVertexCount=0;
|
NewVertexCount=0;
|
||||||
for (int idx=0;idx<12;idx++)
|
for (int idx=0;idx<12;idx++)
|
||||||
LocalRemap[idx] = -1;
|
LocalRemap[idx] = -1;
|
||||||
|
|
||||||
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++)
|
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++)
|
||||||
{
|
{
|
||||||
if(LocalRemap[triTable[CubeIndex][idx]] == -1)
|
if(LocalRemap[triTable[CubeIndex][idx]] == -1)
|
||||||
{
|
{
|
||||||
NewVertexList[NewVertexCount] = VertexList[triTable[CubeIndex][idx]];
|
NewVertexList[NewVertexCount] = VertexList[triTable[CubeIndex][idx]];
|
||||||
LocalRemap[triTable[CubeIndex][idx]] = NewVertexCount;
|
LocalRemap[triTable[CubeIndex][idx]] = NewVertexCount;
|
||||||
NewVertexCount++;
|
NewVertexCount++;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
for (int idx=0;idx<NewVertexCount;idx++) {
|
for (int idx=0;idx<NewVertexCount;idx++) {
|
||||||
P = NewVertexList[idx];
|
P = NewVertexList[idx];
|
||||||
//P.x += i;
|
//P.x += i;
|
||||||
//P.y += j;
|
//P.y += j;
|
||||||
//P.z += k;
|
//P.z += k;
|
||||||
cellvertices(idx) = P;
|
cellvertices(idx) = P;
|
||||||
}
|
}
|
||||||
nVert = NewVertexCount;
|
nVert = NewVertexCount;
|
||||||
|
|
||||||
TriangleCount = 0;
|
TriangleCount = 0;
|
||||||
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
|
for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
|
||||||
Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
|
Triangles(0,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+0]];
|
||||||
Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
|
Triangles(1,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+1]];
|
||||||
Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
|
Triangles(2,TriangleCount) = LocalRemap[triTable[CubeIndex][idx+2]];
|
||||||
TriangleCount++;
|
TriangleCount++;
|
||||||
}
|
}
|
||||||
nTris = TriangleCount;
|
nTris = TriangleCount;
|
||||||
|
|
||||||
// Now add the local values to the DECL data structure
|
|
||||||
IntArray HalfEdge(6,nTris*3);
|
|
||||||
DoubleArray EdgeAngles(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);
|
|
||||||
// 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
|
|
||||||
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
|
|
||||||
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
|
|
||||||
idx_edge++;
|
|
||||||
}
|
|
||||||
int EdgeCount=idx_edge;
|
|
||||||
for (int idx=0; idx<EdgeCount; idx++){
|
|
||||||
int V1=HalfEdge(0,idx);
|
|
||||||
int V2=HalfEdge(1,idx);
|
|
||||||
// Find all the twins within the cube
|
|
||||||
for (int jdx=0; idx<EdgeCount; jdx++){
|
|
||||||
if (HalfEdge(1,jdx) == V1 && HalfEdge(0,jdx) == V2){
|
|
||||||
// this is the pair
|
|
||||||
HalfEdge(3,idx) = jdx;
|
|
||||||
HalfEdge(3,jdx) = idx;
|
|
||||||
}
|
|
||||||
if (HalfEdge(1,jdx) == V2 && HalfEdge(0,jdx) == 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
|
|
||||||
}
|
|
||||||
// Map vertices to global coordinates
|
|
||||||
for (int idx=0;idx<NewVertexCount;idx++) {
|
|
||||||
P = cellvertices(idx);
|
|
||||||
P.x += i;
|
|
||||||
P.y += j;
|
|
||||||
P.z += k;
|
|
||||||
cellvertices(idx) = P;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
// Now add the local values to the DECL data structure
|
||||||
|
halfedge.data.resize(6,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);
|
||||||
|
// first edge: V1->V2
|
||||||
|
halfedge.data(0,idx_edge) = V1; // first vertex
|
||||||
|
halfedge.data(1,idx_edge) = V2; // second vertex
|
||||||
|
halfedge.data(2,idx_edge) = idx; // triangle
|
||||||
|
halfedge.data(3,idx_edge) = -1; // twin
|
||||||
|
halfedge.data(4,idx_edge) = idx_edge+2; // previous edge
|
||||||
|
halfedge.data(5,idx_edge) = idx_edge+1; // next edge
|
||||||
|
idx_edge++;
|
||||||
|
// second edge: V2->V3
|
||||||
|
halfedge.data(0,idx_edge) = V2; // first vertex
|
||||||
|
halfedge.data(1,idx_edge) = V3; // second vertex
|
||||||
|
halfedge.data(2,idx_edge) = idx; // triangle
|
||||||
|
halfedge.data(3,idx_edge) = -1; // twin
|
||||||
|
halfedge.data(4,idx_edge) = idx_edge-1; // previous edge
|
||||||
|
halfedge.data(5,idx_edge) = idx_edge+1; // next edge
|
||||||
|
idx_edge++;
|
||||||
|
// third edge: V3->V1
|
||||||
|
halfedge.data(0,idx_edge) = V3; // first vertex
|
||||||
|
halfedge.data(1,idx_edge) = V1; // second vertex
|
||||||
|
halfedge.data(2,idx_edge) = idx; // triangle
|
||||||
|
halfedge.data(3,idx_edge) = -1; // twin
|
||||||
|
halfedge.data(4,idx_edge) = idx_edge-1; // previous edge
|
||||||
|
halfedge.data(5,idx_edge) = idx_edge-2; // next edge
|
||||||
|
idx_edge++;
|
||||||
|
}
|
||||||
|
int EdgeCount=idx_edge;
|
||||||
|
for (int idx=0; idx<EdgeCount; idx++){
|
||||||
|
int V1=halfedge.data(0,idx);
|
||||||
|
int V2=halfedge.data(1,idx);
|
||||||
|
// Find all the twins within the cube
|
||||||
|
for (int jdx=0; idx<EdgeCount; jdx++){
|
||||||
|
if (halfedge.data(1,jdx) == V1 && halfedge.data(0,jdx) == V2){
|
||||||
|
// this is the pair
|
||||||
|
halfedge.data(3,idx) = jdx;
|
||||||
|
halfedge.data(3,jdx) = idx;
|
||||||
|
}
|
||||||
|
if (halfedge.data(1,jdx) == V2 && halfedge.data(0,jdx) == 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.data(3,idx_edge) = -1; // ghost twin for x=0 face
|
||||||
|
if (P.x == 1.0 && Q.x == 1.0) halfedge.data(3,idx_edge) = -2; // ghost twin for x=1 face
|
||||||
|
if (P.y == 0.0 && Q.y == 0.0) halfedge.data(3,idx_edge) = -3; // ghost twin for y=0 face
|
||||||
|
if (P.y == 1.0 && Q.y == 1.0) halfedge.data(3,idx_edge) = -4; // ghost twin for y=1 face
|
||||||
|
if (P.z == 0.0 && Q.z == 0.0) halfedge.data(3,idx_edge) = -5; // ghost twin for z=0 face
|
||||||
|
if (P.z == 1.0 && Q.z == 1.0) halfedge.data(3,idx_edge) = -6; // ghost twin for z=1 face
|
||||||
|
}
|
||||||
|
// Map vertices to global coordinates
|
||||||
|
for (int idx=0;idx<NewVertexCount;idx++) {
|
||||||
|
P = cellvertices(idx);
|
||||||
|
P.x += i;
|
||||||
|
P.y += j;
|
||||||
|
P.z += k;
|
||||||
|
cellvertices(idx) = P;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
Point DECL::TriNormal(int edge)
|
Point DECL::TriNormal(int edge)
|
||||||
@ -281,9 +288,15 @@ Point DECL::TriNormal(int edge)
|
|||||||
Point P,Q;
|
Point P,Q;
|
||||||
double ux,uy,uz,vx,vy,vz;
|
double ux,uy,uz,vx,vy,vz;
|
||||||
double nx,ny,nz,len;
|
double nx,ny,nz,len;
|
||||||
if (edge == -1) P.x = 1.0; P.y = 0.0; P.z = 0.0; // x cube face
|
if (edge == -1){
|
||||||
else if (edge == -2) P.x = 0.0; P.y = 1.0; P.z = 0.0; // y cube face
|
P.x = 1.0; P.y = 0.0; P.z = 0.0; // x cube face
|
||||||
else if (edge == -3) P.x = 0.0; P.y = 0.0; P.z = 1.0; // z cube face
|
}
|
||||||
|
else if (edge == -2){
|
||||||
|
P.x = 0.0; P.y = 1.0; P.z = 0.0; // y cube face
|
||||||
|
}
|
||||||
|
else if (edge == -3){
|
||||||
|
P.x = 0.0; P.y = 0.0; P.z = 1.0; // z cube face
|
||||||
|
}
|
||||||
else{
|
else{
|
||||||
// coordinates for first edge
|
// coordinates for first edge
|
||||||
P = vertex.coords(halfedge.v1(edge));
|
P = vertex.coords(halfedge.v1(edge));
|
||||||
@ -325,7 +338,7 @@ void Isosurface(DoubleArray &A, const double &v)
|
|||||||
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
||||||
|
|
||||||
int TriangleCount;
|
int TriangleCount;
|
||||||
int NewVertexCount;
|
int NewVertexCount;
|
||||||
int CubeIndex;
|
int CubeIndex;
|
||||||
int nTris, nVert;
|
int nTris, nVert;
|
||||||
|
|
||||||
|
|||||||
@ -32,8 +32,8 @@ public:
|
|||||||
unsigned long int next(unsigned long int edge);
|
unsigned long int next(unsigned long int edge);
|
||||||
unsigned long int prev(unsigned long int edge);
|
unsigned long int prev(unsigned long int edge);
|
||||||
|
|
||||||
|
Array<unsigned long int> data;
|
||||||
private:
|
private:
|
||||||
Array<unsigned long int> HalfEdge;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
// DECL
|
// DECL
|
||||||
@ -45,7 +45,7 @@ public:
|
|||||||
unsigned long int face();
|
unsigned long int face();
|
||||||
Vertex vertex;
|
Vertex vertex;
|
||||||
Halfedge halfedge;
|
Halfedge halfedge;
|
||||||
void AddCube(); // need a function to add new faces based on marching cubes surface
|
void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k);
|
||||||
|
|
||||||
double origin(int edge);
|
double origin(int edge);
|
||||||
double EdgeAngle(int edge);
|
double EdgeAngle(int edge);
|
||||||
|
|||||||
@ -16,23 +16,23 @@
|
|||||||
|
|
||||||
AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)
|
AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)
|
||||||
{
|
{
|
||||||
lhs = static_cast<AnalysisType>(
|
lhs = static_cast<AnalysisType>(
|
||||||
static_cast<std::underlying_type<AnalysisType>::type>(lhs) |
|
static_cast<std::underlying_type<AnalysisType>::type>(lhs) |
|
||||||
static_cast<std::underlying_type<AnalysisType>::type>(rhs)
|
static_cast<std::underlying_type<AnalysisType>::type>(rhs)
|
||||||
);
|
);
|
||||||
return lhs;
|
return lhs;
|
||||||
}
|
}
|
||||||
bool matches( AnalysisType x, AnalysisType y )
|
bool matches( AnalysisType x, AnalysisType y )
|
||||||
{
|
{
|
||||||
return static_cast<std::underlying_type<AnalysisType>::type>(x) &
|
return static_cast<std::underlying_type<AnalysisType>::type>(x) &
|
||||||
static_cast<std::underlying_type<AnalysisType>::type>(y) != 0;
|
static_cast<std::underlying_type<AnalysisType>::type>(y) != 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
template<class TYPE>
|
template<class TYPE>
|
||||||
void DeleteArray( const TYPE *p )
|
void DeleteArray( const TYPE *p )
|
||||||
{
|
{
|
||||||
delete [] p;
|
delete [] p;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -40,32 +40,32 @@ void DeleteArray( const TYPE *p )
|
|||||||
class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
|
class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cphi_, std::shared_ptr<double> cfq_, int N_ ):
|
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cphi_, std::shared_ptr<double> cfq_, int N_ ):
|
||||||
filename(filename_), cphi(cphi_), cfq(cfq_), N(N_) {}
|
filename(filename_), cphi(cphi_), cfq(cfq_), N(N_) {}
|
||||||
virtual void run() {
|
virtual void run() {
|
||||||
PROFILE_START("Save Checkpoint",1);
|
PROFILE_START("Save Checkpoint",1);
|
||||||
double value;
|
double value;
|
||||||
ofstream File(filename,ios::binary);
|
ofstream File(filename,ios::binary);
|
||||||
for (int n=0; n<N; n++){
|
for (int n=0; n<N; n++){
|
||||||
// Write the two density values
|
// Write the two density values
|
||||||
value = cphi.get()[n];
|
value = cphi.get()[n];
|
||||||
File.write((char*) &value, sizeof(value));
|
File.write((char*) &value, sizeof(value));
|
||||||
// Write the distributions
|
// Write the distributions
|
||||||
for (int q=0; q<19; q++){
|
for (int q=0; q<19; q++){
|
||||||
value = cfq.get()[q*N+n];
|
value = cfq.get()[q*N+n];
|
||||||
File.write((char*) &value, sizeof(value));
|
File.write((char*) &value, sizeof(value));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
File.close();
|
File.close();
|
||||||
PROFILE_STOP("Save Checkpoint",1);
|
PROFILE_STOP("Save Checkpoint",1);
|
||||||
};
|
};
|
||||||
private:
|
private:
|
||||||
WriteRestartWorkItem();
|
WriteRestartWorkItem();
|
||||||
const char* filename;
|
const char* filename;
|
||||||
std::shared_ptr<double> cfq,cphi;
|
std::shared_ptr<double> cfq,cphi;
|
||||||
// const DoubleArray& phase;
|
// const DoubleArray& phase;
|
||||||
//const DoubleArray& dist;
|
//const DoubleArray& dist;
|
||||||
const int N;
|
const int N;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
@ -74,78 +74,78 @@ static const std::string id_map_filename = "lbpm_id_map.txt";
|
|||||||
class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet<void>
|
class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet<void>
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
||||||
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
||||||
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ):
|
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ):
|
||||||
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
||||||
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
|
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
~BlobIdentificationWorkItem1() { }
|
~BlobIdentificationWorkItem1() { }
|
||||||
virtual void run() {
|
virtual void run() {
|
||||||
// Compute the global blob id and compare to the previous version
|
// Compute the global blob id and compare to the previous version
|
||||||
PROFILE_START("Identify blobs",1);
|
PROFILE_START("Identify blobs",1);
|
||||||
double vF = 0.0;
|
double vF = 0.0;
|
||||||
double vS = -1.0; // one voxel buffer region around solid
|
double vS = -1.0; // one voxel buffer region around solid
|
||||||
IntArray& ids = new_index->second;
|
IntArray& ids = new_index->second;
|
||||||
new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm);
|
new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm);
|
||||||
PROFILE_STOP("Identify blobs",1);
|
PROFILE_STOP("Identify blobs",1);
|
||||||
}
|
}
|
||||||
private:
|
private:
|
||||||
BlobIdentificationWorkItem1();
|
BlobIdentificationWorkItem1();
|
||||||
int timestep;
|
int timestep;
|
||||||
int Nx, Ny, Nz;
|
int Nx, Ny, Nz;
|
||||||
const RankInfoStruct& rank_info;
|
const RankInfoStruct& rank_info;
|
||||||
std::shared_ptr<const DoubleArray> phase;
|
std::shared_ptr<const DoubleArray> phase;
|
||||||
const DoubleArray& dist;
|
const DoubleArray& dist;
|
||||||
BlobIDstruct last_id, new_index, new_id;
|
BlobIDstruct last_id, new_index, new_id;
|
||||||
BlobIDList new_list;
|
BlobIDList new_list;
|
||||||
runAnalysis::commWrapper comm;
|
runAnalysis::commWrapper comm;
|
||||||
};
|
};
|
||||||
class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet<void>
|
class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet<void>
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
||||||
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
||||||
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ):
|
BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ):
|
||||||
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
||||||
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
|
phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
~BlobIdentificationWorkItem2() { }
|
~BlobIdentificationWorkItem2() { }
|
||||||
virtual void run() {
|
virtual void run() {
|
||||||
// Compute the global blob id and compare to the previous version
|
// Compute the global blob id and compare to the previous version
|
||||||
PROFILE_START("Identify blobs maps",1);
|
PROFILE_START("Identify blobs maps",1);
|
||||||
const IntArray& ids = new_index->second;
|
const IntArray& ids = new_index->second;
|
||||||
static int max_id = -1;
|
static int max_id = -1;
|
||||||
new_id->first = new_index->first;
|
new_id->first = new_index->first;
|
||||||
new_id->second = new_index->second;
|
new_id->second = new_index->second;
|
||||||
if ( last_id.get()!=NULL ) {
|
if ( last_id.get()!=NULL ) {
|
||||||
// Compute the timestep-timestep map
|
// Compute the timestep-timestep map
|
||||||
const IntArray& old_ids = last_id->second;
|
const IntArray& old_ids = last_id->second;
|
||||||
ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm);
|
ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm);
|
||||||
// Renumber the current timestep's ids
|
// Renumber the current timestep's ids
|
||||||
getNewIDs(map,max_id,*new_list);
|
getNewIDs(map,max_id,*new_list);
|
||||||
renumberIDs(*new_list,new_id->second);
|
renumberIDs(*new_list,new_id->second);
|
||||||
writeIDMap(map,timestep,id_map_filename);
|
writeIDMap(map,timestep,id_map_filename);
|
||||||
} else {
|
} else {
|
||||||
max_id = -1;
|
max_id = -1;
|
||||||
ID_map_struct map(new_id->first);
|
ID_map_struct map(new_id->first);
|
||||||
getNewIDs(map,max_id,*new_list);
|
getNewIDs(map,max_id,*new_list);
|
||||||
writeIDMap(map,timestep,id_map_filename);
|
writeIDMap(map,timestep,id_map_filename);
|
||||||
}
|
}
|
||||||
PROFILE_STOP("Identify blobs maps",1);
|
PROFILE_STOP("Identify blobs maps",1);
|
||||||
}
|
}
|
||||||
private:
|
private:
|
||||||
BlobIdentificationWorkItem2();
|
BlobIdentificationWorkItem2();
|
||||||
int timestep;
|
int timestep;
|
||||||
int Nx, Ny, Nz;
|
int Nx, Ny, Nz;
|
||||||
const RankInfoStruct& rank_info;
|
const RankInfoStruct& rank_info;
|
||||||
std::shared_ptr<const DoubleArray> phase;
|
std::shared_ptr<const DoubleArray> phase;
|
||||||
const DoubleArray& dist;
|
const DoubleArray& dist;
|
||||||
BlobIDstruct last_id, new_index, new_id;
|
BlobIDstruct last_id, new_index, new_id;
|
||||||
BlobIDList new_list;
|
BlobIDList new_list;
|
||||||
runAnalysis::commWrapper comm;
|
runAnalysis::commWrapper comm;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
@ -153,36 +153,36 @@ private:
|
|||||||
class WriteVisWorkItem: public ThreadPool::WorkItemRet<void>
|
class WriteVisWorkItem: public ThreadPool::WorkItemRet<void>
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
|
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
|
||||||
TwoPhase& Avgerages_, fillHalo<double>& fillData_, runAnalysis::commWrapper&& comm_ ):
|
TwoPhase& Avgerages_, fillHalo<double>& fillData_, runAnalysis::commWrapper&& comm_ ):
|
||||||
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_))
|
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_))
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
~WriteVisWorkItem() { }
|
~WriteVisWorkItem() { }
|
||||||
virtual void run() {
|
virtual void run() {
|
||||||
PROFILE_START("Save Vis",1);
|
PROFILE_START("Save Vis",1);
|
||||||
ASSERT(visData[0].vars[0]->name=="phase");
|
ASSERT(visData[0].vars[0]->name=="phase");
|
||||||
ASSERT(visData[0].vars[1]->name=="Pressure");
|
ASSERT(visData[0].vars[1]->name=="Pressure");
|
||||||
ASSERT(visData[0].vars[2]->name=="SignDist");
|
ASSERT(visData[0].vars[2]->name=="SignDist");
|
||||||
ASSERT(visData[0].vars[3]->name=="BlobID");
|
ASSERT(visData[0].vars[3]->name=="BlobID");
|
||||||
Array<double>& PhaseData = visData[0].vars[0]->data;
|
Array<double>& PhaseData = visData[0].vars[0]->data;
|
||||||
Array<double>& PressData = visData[0].vars[1]->data;
|
Array<double>& PressData = visData[0].vars[1]->data;
|
||||||
Array<double>& SignData = visData[0].vars[2]->data;
|
Array<double>& SignData = visData[0].vars[2]->data;
|
||||||
Array<double>& BlobData = visData[0].vars[3]->data;
|
Array<double>& BlobData = visData[0].vars[3]->data;
|
||||||
fillData.copy(Averages.SDn,PhaseData);
|
fillData.copy(Averages.SDn,PhaseData);
|
||||||
fillData.copy(Averages.Press,PressData);
|
fillData.copy(Averages.Press,PressData);
|
||||||
fillData.copy(Averages.SDs,SignData);
|
fillData.copy(Averages.SDs,SignData);
|
||||||
fillData.copy(Averages.Label_NWP,BlobData);
|
fillData.copy(Averages.Label_NWP,BlobData);
|
||||||
IO::writeData( timestep, visData, comm.comm );
|
IO::writeData( timestep, visData, comm.comm );
|
||||||
PROFILE_STOP("Save Vis",1);
|
PROFILE_STOP("Save Vis",1);
|
||||||
};
|
};
|
||||||
private:
|
private:
|
||||||
WriteVisWorkItem();
|
WriteVisWorkItem();
|
||||||
int timestep;
|
int timestep;
|
||||||
std::vector<IO::MeshDataStruct>& visData;
|
std::vector<IO::MeshDataStruct>& visData;
|
||||||
TwoPhase& Averages;
|
TwoPhase& Averages;
|
||||||
fillHalo<double>& fillData;
|
fillHalo<double>& fillData;
|
||||||
runAnalysis::commWrapper comm;
|
runAnalysis::commWrapper comm;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
@ -191,46 +191,46 @@ private:
|
|||||||
class AnalysisWorkItem: public ThreadPool::WorkItemRet<void>
|
class AnalysisWorkItem: public ThreadPool::WorkItemRet<void>
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_,
|
AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_,
|
||||||
BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
|
BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
|
||||||
type(type_), timestep(timestep_), Averages(Averages_),
|
type(type_), timestep(timestep_), Averages(Averages_),
|
||||||
blob_ids(ids), id_list(id_list_), beta(beta_) { }
|
blob_ids(ids), id_list(id_list_), beta(beta_) { }
|
||||||
~AnalysisWorkItem() { }
|
~AnalysisWorkItem() { }
|
||||||
virtual void run() {
|
virtual void run() {
|
||||||
Averages.NumberComponents_NWP = blob_ids->first;
|
Averages.NumberComponents_NWP = blob_ids->first;
|
||||||
Averages.Label_NWP = blob_ids->second;
|
Averages.Label_NWP = blob_ids->second;
|
||||||
Averages.Label_NWP_map = *id_list;
|
Averages.Label_NWP_map = *id_list;
|
||||||
Averages.NumberComponents_WP = 1;
|
Averages.NumberComponents_WP = 1;
|
||||||
Averages.Label_WP.fill(0.0);
|
Averages.Label_WP.fill(0.0);
|
||||||
if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
|
if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
|
||||||
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
|
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
|
||||||
}
|
}
|
||||||
if ( matches(type,AnalysisType::ComputeAverages) ) {
|
if ( matches(type,AnalysisType::ComputeAverages) ) {
|
||||||
PROFILE_START("Compute dist",1);
|
PROFILE_START("Compute dist",1);
|
||||||
Averages.Initialize();
|
Averages.Initialize();
|
||||||
Averages.ComputeDelPhi();
|
Averages.ComputeDelPhi();
|
||||||
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
|
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
|
||||||
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
|
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
|
||||||
Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
|
Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
|
||||||
Averages.UpdateMeshValues();
|
Averages.UpdateMeshValues();
|
||||||
Averages.ComputeLocal();
|
Averages.ComputeLocal();
|
||||||
Averages.Reduce();
|
Averages.Reduce();
|
||||||
Averages.PrintAll(timestep);
|
Averages.PrintAll(timestep);
|
||||||
Averages.Initialize();
|
Averages.Initialize();
|
||||||
Averages.ComponentAverages();
|
Averages.ComponentAverages();
|
||||||
Averages.SortBlobs();
|
Averages.SortBlobs();
|
||||||
Averages.PrintComponents(timestep);
|
Averages.PrintComponents(timestep);
|
||||||
PROFILE_STOP("Compute dist",1);
|
PROFILE_STOP("Compute dist",1);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
private:
|
private:
|
||||||
AnalysisWorkItem();
|
AnalysisWorkItem();
|
||||||
AnalysisType type;
|
AnalysisType type;
|
||||||
int timestep;
|
int timestep;
|
||||||
TwoPhase& Averages;
|
TwoPhase& Averages;
|
||||||
BlobIDstruct blob_ids;
|
BlobIDstruct blob_ids;
|
||||||
BlobIDList id_list;
|
BlobIDList id_list;
|
||||||
double beta;
|
double beta;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
@ -238,44 +238,44 @@ private:
|
|||||||
* MPI comm wrapper for use with analysis *
|
* MPI comm wrapper for use with analysis *
|
||||||
******************************************************************/
|
******************************************************************/
|
||||||
runAnalysis::commWrapper::commWrapper( int tag_, MPI_Comm comm_, runAnalysis* analysis_ ):
|
runAnalysis::commWrapper::commWrapper( int tag_, MPI_Comm comm_, runAnalysis* analysis_ ):
|
||||||
comm(comm_),
|
comm(comm_),
|
||||||
tag(tag_),
|
tag(tag_),
|
||||||
analysis(analysis_)
|
analysis(analysis_)
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ):
|
runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ):
|
||||||
comm(rhs.comm),
|
comm(rhs.comm),
|
||||||
tag(rhs.tag),
|
tag(rhs.tag),
|
||||||
analysis(rhs.analysis)
|
analysis(rhs.analysis)
|
||||||
{
|
{
|
||||||
rhs.tag = -1;
|
rhs.tag = -1;
|
||||||
}
|
}
|
||||||
runAnalysis::commWrapper::~commWrapper()
|
runAnalysis::commWrapper::~commWrapper()
|
||||||
{
|
{
|
||||||
if ( tag == -1 )
|
if ( tag == -1 )
|
||||||
return;
|
return;
|
||||||
MPI_Barrier( comm );
|
MPI_Barrier( comm );
|
||||||
analysis->d_comm_used[tag] = false;
|
analysis->d_comm_used[tag] = false;
|
||||||
}
|
}
|
||||||
runAnalysis::commWrapper runAnalysis::getComm( )
|
runAnalysis::commWrapper runAnalysis::getComm( )
|
||||||
{
|
{
|
||||||
// Get a tag from root
|
// Get a tag from root
|
||||||
int tag = -1;
|
int tag = -1;
|
||||||
if ( d_rank == 0 ) {
|
if ( d_rank == 0 ) {
|
||||||
for (int i=0; i<1024; i++) {
|
for (int i=0; i<1024; i++) {
|
||||||
if ( !d_comm_used[i] ) {
|
if ( !d_comm_used[i] ) {
|
||||||
tag = i;
|
tag = i;
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
if ( tag == -1 )
|
if ( tag == -1 )
|
||||||
ERROR("Unable to get comm");
|
ERROR("Unable to get comm");
|
||||||
}
|
}
|
||||||
MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm );
|
MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm );
|
||||||
d_comm_used[tag] = true;
|
d_comm_used[tag] = true;
|
||||||
if ( d_comms[tag] == MPI_COMM_NULL )
|
if ( d_comms[tag] == MPI_COMM_NULL )
|
||||||
MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] );
|
MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] );
|
||||||
return commWrapper(tag,d_comms[tag],this);
|
return commWrapper(tag,d_comms[tag],this);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -283,29 +283,29 @@ runAnalysis::commWrapper runAnalysis::getComm( )
|
|||||||
* Constructor/Destructors *
|
* Constructor/Destructors *
|
||||||
******************************************************************/
|
******************************************************************/
|
||||||
runAnalysis::runAnalysis( std::shared_ptr<Database> db,
|
runAnalysis::runAnalysis( std::shared_ptr<Database> db,
|
||||||
const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
|
const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
|
||||||
int Np, bool Regular, double beta, IntArray Map ):
|
int Np, bool Regular, double beta, IntArray Map ):
|
||||||
d_Np( Np ),
|
d_Np( Np ),
|
||||||
d_beta( beta ),
|
d_beta( beta ),
|
||||||
d_regular ( Regular),
|
d_regular ( Regular),
|
||||||
d_rank_info( rank_info ),
|
d_rank_info( rank_info ),
|
||||||
d_Map( Map ),
|
d_Map( Map ),
|
||||||
d_ScaLBL_Comm( ScaLBL_Comm),
|
d_ScaLBL_Comm( ScaLBL_Comm),
|
||||||
d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1)
|
d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1)
|
||||||
{
|
{
|
||||||
|
|
||||||
char rankString[20];
|
char rankString[20];
|
||||||
sprintf(rankString,"%05d",Dm->rank());
|
sprintf(rankString,"%05d",Dm->rank());
|
||||||
d_N[0] = Dm->Nx;
|
d_N[0] = Dm->Nx;
|
||||||
d_N[1] = Dm->Ny;
|
d_N[1] = Dm->Ny;
|
||||||
d_N[2] = Dm->Nz;
|
d_N[2] = Dm->Nz;
|
||||||
d_restart_interval = db->getScalar<int>( "restart_interval" );
|
d_restart_interval = db->getScalar<int>( "restart_interval" );
|
||||||
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
|
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
|
||||||
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
|
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
|
||||||
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
|
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
|
||||||
auto restart_file = db->getScalar<std::string>( "restart_file" );
|
auto restart_file = db->getScalar<std::string>( "restart_file" );
|
||||||
d_restartFile = restart_file + "." + rankString;
|
d_restartFile = restart_file + "." + rankString;
|
||||||
d_rank = MPI_WORLD_RANK();
|
d_rank = MPI_WORLD_RANK();
|
||||||
writeIDMap(ID_map_struct(),0,id_map_filename);
|
writeIDMap(ID_map_struct(),0,id_map_filename);
|
||||||
// Initialize IO for silo
|
// Initialize IO for silo
|
||||||
IO::initialize("","silo","false");
|
IO::initialize("","silo","false");
|
||||||
@ -337,41 +337,41 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
|
|||||||
BlobIDVar->dim = 1;
|
BlobIDVar->dim = 1;
|
||||||
BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
|
BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
|
||||||
d_meshData[0].vars.push_back(BlobIDVar);
|
d_meshData[0].vars.push_back(BlobIDVar);
|
||||||
// Initialize the comms
|
// Initialize the comms
|
||||||
MPI_Comm_dup(MPI_COMM_WORLD,&d_comm);
|
MPI_Comm_dup(MPI_COMM_WORLD,&d_comm);
|
||||||
for (int i=0; i<1024; i++) {
|
for (int i=0; i<1024; i++) {
|
||||||
d_comms[i] = MPI_COMM_NULL;
|
d_comms[i] = MPI_COMM_NULL;
|
||||||
d_comm_used[i] = false;
|
d_comm_used[i] = false;
|
||||||
}
|
}
|
||||||
// Initialize the threads
|
// Initialize the threads
|
||||||
int N_threads = db->getWithDefault<int>( "N_threads", 4 );
|
int N_threads = db->getWithDefault<int>( "N_threads", 4 );
|
||||||
auto method = db->getWithDefault<std::string>( "load_balance", "default" );
|
auto method = db->getWithDefault<std::string>( "load_balance", "default" );
|
||||||
createThreads( method, N_threads );
|
createThreads( method, N_threads );
|
||||||
}
|
}
|
||||||
runAnalysis::~runAnalysis( )
|
runAnalysis::~runAnalysis( )
|
||||||
{
|
{
|
||||||
// Finish processing analysis
|
// Finish processing analysis
|
||||||
finish();
|
finish();
|
||||||
// Clear internal data
|
// Clear internal data
|
||||||
MPI_Comm_free( &d_comm );
|
MPI_Comm_free( &d_comm );
|
||||||
for (int i=0; i<1024; i++) {
|
for (int i=0; i<1024; i++) {
|
||||||
if ( d_comms[i] != MPI_COMM_NULL )
|
if ( d_comms[i] != MPI_COMM_NULL )
|
||||||
MPI_Comm_free(&d_comms[i]);
|
MPI_Comm_free(&d_comms[i]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
void runAnalysis::finish( )
|
void runAnalysis::finish( )
|
||||||
{
|
{
|
||||||
PROFILE_START("finish");
|
PROFILE_START("finish");
|
||||||
// Wait for the work items to finish
|
// Wait for the work items to finish
|
||||||
d_tpool.wait_pool_finished();
|
d_tpool.wait_pool_finished();
|
||||||
// Clear the wait ids
|
// Clear the wait ids
|
||||||
d_wait_blobID.reset();
|
d_wait_blobID.reset();
|
||||||
d_wait_analysis.reset();
|
d_wait_analysis.reset();
|
||||||
d_wait_vis.reset();
|
d_wait_vis.reset();
|
||||||
d_wait_restart.reset();
|
d_wait_restart.reset();
|
||||||
// Syncronize
|
// Syncronize
|
||||||
MPI_Barrier( d_comm );
|
MPI_Barrier( d_comm );
|
||||||
PROFILE_STOP("finish");
|
PROFILE_STOP("finish");
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -380,50 +380,50 @@ void runAnalysis::finish( )
|
|||||||
******************************************************************/
|
******************************************************************/
|
||||||
void print( const std::vector<int>& ids )
|
void print( const std::vector<int>& ids )
|
||||||
{
|
{
|
||||||
if ( ids.empty() )
|
if ( ids.empty() )
|
||||||
return;
|
return;
|
||||||
printf("%i",ids[0]);
|
printf("%i",ids[0]);
|
||||||
for (size_t i=1; i<ids.size(); i++)
|
for (size_t i=1; i<ids.size(); i++)
|
||||||
printf(", %i",ids[i]);
|
printf(", %i",ids[i]);
|
||||||
printf("\n");
|
printf("\n");
|
||||||
}
|
}
|
||||||
void runAnalysis::createThreads( const std::string& method, int N_threads )
|
void runAnalysis::createThreads( const std::string& method, int N_threads )
|
||||||
{
|
{
|
||||||
// Check if we are not using analysis threads
|
// Check if we are not using analysis threads
|
||||||
if ( method == "none" )
|
if ( method == "none" )
|
||||||
return;
|
return;
|
||||||
// Check if we have thread support
|
// Check if we have thread support
|
||||||
int thread_support;
|
int thread_support;
|
||||||
MPI_Query_thread( &thread_support );
|
MPI_Query_thread( &thread_support );
|
||||||
if ( thread_support < MPI_THREAD_MULTIPLE ) {
|
if ( thread_support < MPI_THREAD_MULTIPLE ) {
|
||||||
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
|
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
// Create the threads
|
// Create the threads
|
||||||
const auto cores = d_tpool.getProcessAffinity();
|
const auto cores = d_tpool.getProcessAffinity();
|
||||||
if ( cores.empty() ) {
|
if ( cores.empty() ) {
|
||||||
// We were not able to get the cores for the process
|
// We were not able to get the cores for the process
|
||||||
d_tpool.setNumThreads( N_threads );
|
d_tpool.setNumThreads( N_threads );
|
||||||
} else if ( method == "default" ) {
|
} else if ( method == "default" ) {
|
||||||
// Create the given number of threads, but let the OS manage affinities
|
// Create the given number of threads, but let the OS manage affinities
|
||||||
d_tpool.setNumThreads( N_threads );
|
d_tpool.setNumThreads( N_threads );
|
||||||
} else if ( method == "independent" ) {
|
} else if ( method == "independent" ) {
|
||||||
int N = cores.size() - 1;
|
int N = cores.size() - 1;
|
||||||
d_tpool.setNumThreads( N );
|
d_tpool.setNumThreads( N );
|
||||||
d_tpool.setThreadAffinity( { cores[0] } );
|
d_tpool.setThreadAffinity( { cores[0] } );
|
||||||
for ( int i=0; i<N; i++)
|
for ( int i=0; i<N; i++)
|
||||||
d_tpool.setThreadAffinity( i, { cores[i+1] } );
|
d_tpool.setThreadAffinity( i, { cores[i+1] } );
|
||||||
}
|
}
|
||||||
// Print the current affinities
|
// Print the current affinities
|
||||||
if ( d_rank == 0 ) {
|
if ( d_rank == 0 ) {
|
||||||
printf("Affinities - rank 0:\n");
|
printf("Affinities - rank 0:\n");
|
||||||
printf("Main: ");
|
printf("Main: ");
|
||||||
print(d_tpool.getProcessAffinity());
|
print(d_tpool.getProcessAffinity());
|
||||||
for (int i=0; i<d_tpool.getNumThreads(); i++) {
|
for (int i=0; i<d_tpool.getNumThreads(); i++) {
|
||||||
printf("Thread %i: ",i+1);
|
printf("Thread %i: ",i+1);
|
||||||
print(d_tpool.getThreadAffinity(i));
|
print(d_tpool.getThreadAffinity(i));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -432,17 +432,17 @@ void runAnalysis::createThreads( const std::string& method, int N_threads )
|
|||||||
******************************************************************/
|
******************************************************************/
|
||||||
AnalysisType runAnalysis::computeAnalysisType( int timestep )
|
AnalysisType runAnalysis::computeAnalysisType( int timestep )
|
||||||
{
|
{
|
||||||
AnalysisType type = AnalysisType::AnalyzeNone;
|
AnalysisType type = AnalysisType::AnalyzeNone;
|
||||||
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
|
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
|
||||||
// Copy the phase indicator field for the earlier timestep
|
// Copy the phase indicator field for the earlier timestep
|
||||||
// printf("Copy phase indicator,timestep=%i\n",timestep);
|
// printf("Copy phase indicator,timestep=%i\n",timestep);
|
||||||
type |= AnalysisType::CopyPhaseIndicator;
|
type |= AnalysisType::CopyPhaseIndicator;
|
||||||
}
|
}
|
||||||
if ( timestep%d_blobid_interval == 0 ) {
|
if ( timestep%d_blobid_interval == 0 ) {
|
||||||
// Identify blobs and update global ids in time
|
// Identify blobs and update global ids in time
|
||||||
type |= AnalysisType::IdentifyBlobs;
|
type |= AnalysisType::IdentifyBlobs;
|
||||||
}
|
}
|
||||||
/*#ifdef USE_CUDA
|
/*#ifdef USE_CUDA
|
||||||
if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
|
if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
|
||||||
// Keep a few blob identifications queued up to keep the processors busy,
|
// Keep a few blob identifications queued up to keep the processors busy,
|
||||||
// allowing us to track the blobs as fast as possible
|
// allowing us to track the blobs as fast as possible
|
||||||
@ -450,28 +450,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
|
|||||||
type |= AnalysisType::IdentifyBlobs;
|
type |= AnalysisType::IdentifyBlobs;
|
||||||
}
|
}
|
||||||
#endif */
|
#endif */
|
||||||
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
|
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
|
||||||
// Copy the averages to the CPU (and identify blobs)
|
// Copy the averages to the CPU (and identify blobs)
|
||||||
//printf("Copy sim state, timestep=%i \n",timestep);
|
//printf("Copy sim state, timestep=%i \n",timestep);
|
||||||
type |= AnalysisType::CopySimState;
|
type |= AnalysisType::CopySimState;
|
||||||
type |= AnalysisType::IdentifyBlobs;
|
type |= AnalysisType::IdentifyBlobs;
|
||||||
}
|
}
|
||||||
if ( timestep%d_analysis_interval == 0 ) {
|
if ( timestep%d_analysis_interval == 0 ) {
|
||||||
// Run the analysis
|
// Run the analysis
|
||||||
//printf("Compute averages, timestep=%i \n",timestep);
|
//printf("Compute averages, timestep=%i \n",timestep);
|
||||||
type |= AnalysisType::ComputeAverages;
|
type |= AnalysisType::ComputeAverages;
|
||||||
}
|
}
|
||||||
if (timestep%d_restart_interval == 0) {
|
if (timestep%d_restart_interval == 0) {
|
||||||
// Write the restart file
|
// Write the restart file
|
||||||
type |= AnalysisType::CreateRestart;
|
type |= AnalysisType::CreateRestart;
|
||||||
}
|
}
|
||||||
if (timestep%d_visualization_interval == 0) {
|
if (timestep%d_visualization_interval == 0) {
|
||||||
// Write the visualization data
|
// Write the visualization data
|
||||||
type |= AnalysisType::WriteVis;
|
type |= AnalysisType::WriteVis;
|
||||||
type |= AnalysisType::CopySimState;
|
type |= AnalysisType::CopySimState;
|
||||||
type |= AnalysisType::IdentifyBlobs;
|
type |= AnalysisType::IdentifyBlobs;
|
||||||
}
|
}
|
||||||
return type;
|
return type;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -480,28 +480,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
|
|||||||
* Run the analysis *
|
* Run the analysis *
|
||||||
******************************************************************/
|
******************************************************************/
|
||||||
void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
|
void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
|
||||||
double *Pressure, double *Velocity, double *fq, double *Den)
|
double *Pressure, double *Velocity, double *fq, double *Den)
|
||||||
{
|
{
|
||||||
int N = d_N[0]*d_N[1]*d_N[2];
|
int N = d_N[0]*d_N[1]*d_N[2];
|
||||||
|
|
||||||
// Check which analysis steps we need to perform
|
// Check which analysis steps we need to perform
|
||||||
auto type = computeAnalysisType( timestep );
|
auto type = computeAnalysisType( timestep );
|
||||||
if ( type == AnalysisType::AnalyzeNone )
|
if ( type == AnalysisType::AnalyzeNone )
|
||||||
return;
|
return;
|
||||||
|
|
||||||
// Check how may queued items we have
|
// Check how may queued items we have
|
||||||
if ( d_tpool.N_queued() > 20 ) {
|
if ( d_tpool.N_queued() > 20 ) {
|
||||||
std::cerr << "Analysis queue is getting behind, waiting ...\n";
|
std::cerr << "Analysis queue is getting behind, waiting ...\n";
|
||||||
finish();
|
finish();
|
||||||
}
|
}
|
||||||
|
|
||||||
PROFILE_START("run");
|
PROFILE_START("run");
|
||||||
|
|
||||||
// Copy the appropriate variables to the host (so we can spawn new threads)
|
// Copy the appropriate variables to the host (so we can spawn new threads)
|
||||||
ScaLBL_DeviceBarrier();
|
ScaLBL_DeviceBarrier();
|
||||||
PROFILE_START("Copy data to host",1);
|
PROFILE_START("Copy data to host",1);
|
||||||
std::shared_ptr<DoubleArray> phase;
|
std::shared_ptr<DoubleArray> phase;
|
||||||
/* if ( matches(type,AnalysisType::CopyPhaseIndicator) ||
|
/* if ( matches(type,AnalysisType::CopyPhaseIndicator) ||
|
||||||
matches(type,AnalysisType::ComputeAverages) ||
|
matches(type,AnalysisType::ComputeAverages) ||
|
||||||
matches(type,AnalysisType::CopySimState) ||
|
matches(type,AnalysisType::CopySimState) ||
|
||||||
matches(type,AnalysisType::IdentifyBlobs) )
|
matches(type,AnalysisType::IdentifyBlobs) )
|
||||||
@ -531,113 +531,113 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
|
|||||||
}
|
}
|
||||||
delete [] TmpDat;
|
delete [] TmpDat;
|
||||||
}
|
}
|
||||||
*/
|
*/
|
||||||
//if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
|
//if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
|
||||||
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
|
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
|
||||||
if (d_regular)
|
if (d_regular)
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus);
|
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus);
|
||||||
else
|
else
|
||||||
ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
|
ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
|
||||||
//memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
|
//memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
|
||||||
}
|
}
|
||||||
if ( timestep%d_analysis_interval == 0 ) {
|
if ( timestep%d_analysis_interval == 0 ) {
|
||||||
if (d_regular)
|
if (d_regular)
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus);
|
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus);
|
||||||
else
|
else
|
||||||
ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
|
ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
|
||||||
//memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
|
//memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
|
||||||
}
|
}
|
||||||
//if ( matches(type,AnalysisType::CopySimState) ) {
|
//if ( matches(type,AnalysisType::CopySimState) ) {
|
||||||
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
|
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
|
||||||
// Copy the members of Averages to the cpu (phase was copied above)
|
// Copy the members of Averages to the cpu (phase was copied above)
|
||||||
PROFILE_START("Copy-Pressure",1);
|
PROFILE_START("Copy-Pressure",1);
|
||||||
ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
|
ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
|
||||||
ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
|
ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
|
||||||
ScaLBL_DeviceBarrier();
|
ScaLBL_DeviceBarrier();
|
||||||
PROFILE_STOP("Copy-Pressure",1);
|
PROFILE_STOP("Copy-Pressure",1);
|
||||||
PROFILE_START("Copy-Wait",1);
|
PROFILE_START("Copy-Wait",1);
|
||||||
PROFILE_STOP("Copy-Wait",1);
|
PROFILE_STOP("Copy-Wait",1);
|
||||||
PROFILE_START("Copy-State",1);
|
PROFILE_START("Copy-State",1);
|
||||||
//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
|
//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
|
||||||
if (d_regular)
|
if (d_regular)
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
|
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
|
||||||
else
|
else
|
||||||
ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
|
ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
|
||||||
// copy other variables
|
// copy other variables
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
|
d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
|
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
|
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
|
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
|
||||||
PROFILE_STOP("Copy-State",1);
|
PROFILE_STOP("Copy-State",1);
|
||||||
}
|
}
|
||||||
std::shared_ptr<double> cfq,cPhi;
|
std::shared_ptr<double> cfq,cPhi;
|
||||||
//if ( matches(type,AnalysisType::CreateRestart) ) {
|
//if ( matches(type,AnalysisType::CreateRestart) ) {
|
||||||
if (timestep%d_restart_interval==0){
|
if (timestep%d_restart_interval==0){
|
||||||
// Copy restart data to the CPU
|
// Copy restart data to the CPU
|
||||||
cPhi = std::shared_ptr<double>(new double[d_Np],DeleteArray<double>);
|
cPhi = std::shared_ptr<double>(new double[d_Np],DeleteArray<double>);
|
||||||
cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
|
cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
|
||||||
ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
|
ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
|
||||||
ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double));
|
ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double));
|
||||||
}
|
}
|
||||||
PROFILE_STOP("Copy data to host",1);
|
PROFILE_STOP("Copy data to host",1);
|
||||||
|
|
||||||
// Spawn threads to do blob identification work
|
// Spawn threads to do blob identification work
|
||||||
if ( matches(type,AnalysisType::IdentifyBlobs) ) {
|
if ( matches(type,AnalysisType::IdentifyBlobs) ) {
|
||||||
phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
|
phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
|
||||||
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase);
|
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase);
|
||||||
|
|
||||||
BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
|
BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
|
||||||
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
|
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
|
||||||
BlobIDList new_list(new std::vector<BlobIDType>());
|
BlobIDList new_list(new std::vector<BlobIDType>());
|
||||||
auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
|
auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
|
||||||
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
|
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
|
||||||
auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
|
auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
|
||||||
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
|
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
|
||||||
work1->add_dependency(d_wait_blobID);
|
work1->add_dependency(d_wait_blobID);
|
||||||
work2->add_dependency(d_tpool.add_work(work1));
|
work2->add_dependency(d_tpool.add_work(work1));
|
||||||
d_wait_blobID = d_tpool.add_work(work2);
|
d_wait_blobID = d_tpool.add_work(work2);
|
||||||
d_last_index = new_index;
|
d_last_index = new_index;
|
||||||
d_last_ids = new_ids;
|
d_last_ids = new_ids;
|
||||||
d_last_id_map = new_list;
|
d_last_id_map = new_list;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Spawn threads to do the analysis work
|
// Spawn threads to do the analysis work
|
||||||
//if (timestep%d_restart_interval==0){
|
//if (timestep%d_restart_interval==0){
|
||||||
// if ( matches(type,AnalysisType::ComputeAverages) ) {
|
// if ( matches(type,AnalysisType::ComputeAverages) ) {
|
||||||
if ( timestep%d_analysis_interval == 0 ) {
|
if ( timestep%d_analysis_interval == 0 ) {
|
||||||
auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
|
auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
|
||||||
work->add_dependency(d_wait_blobID);
|
work->add_dependency(d_wait_blobID);
|
||||||
work->add_dependency(d_wait_analysis);
|
work->add_dependency(d_wait_analysis);
|
||||||
work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying
|
work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying
|
||||||
d_wait_analysis = d_tpool.add_work(work);
|
d_wait_analysis = d_tpool.add_work(work);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Spawn a thread to write the restart file
|
// Spawn a thread to write the restart file
|
||||||
// if ( matches(type,AnalysisType::CreateRestart) ) {
|
// if ( matches(type,AnalysisType::CreateRestart) ) {
|
||||||
if (timestep%d_restart_interval==0){
|
if (timestep%d_restart_interval==0){
|
||||||
|
|
||||||
if (d_rank==0) {
|
if (d_rank==0) {
|
||||||
FILE *Rst = fopen("Restart.txt","w");
|
FILE *Rst = fopen("Restart.txt","w");
|
||||||
fprintf(Rst,"%i\n",timestep+4);
|
fprintf(Rst,"%i\n",timestep+4);
|
||||||
fclose(Rst);
|
fclose(Rst);
|
||||||
}
|
}
|
||||||
// Write the restart file (using a seperate thread)
|
// Write the restart file (using a seperate thread)
|
||||||
auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cPhi,cfq,d_Np);
|
auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cPhi,cfq,d_Np);
|
||||||
work->add_dependency(d_wait_restart);
|
work->add_dependency(d_wait_restart);
|
||||||
d_wait_restart = d_tpool.add_work(work);
|
d_wait_restart = d_tpool.add_work(work);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Save the results for visualization
|
// Save the results for visualization
|
||||||
// if ( matches(type,AnalysisType::CreateRestart) ) {
|
// if ( matches(type,AnalysisType::CreateRestart) ) {
|
||||||
if (timestep%d_restart_interval==0){
|
if (timestep%d_restart_interval==0){
|
||||||
// Write the vis files
|
// Write the vis files
|
||||||
auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
|
auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
|
||||||
work->add_dependency(d_wait_blobID);
|
work->add_dependency(d_wait_blobID);
|
||||||
work->add_dependency(d_wait_analysis);
|
work->add_dependency(d_wait_analysis);
|
||||||
work->add_dependency(d_wait_vis);
|
work->add_dependency(d_wait_vis);
|
||||||
d_wait_vis = d_tpool.add_work(work);
|
d_wait_vis = d_tpool.add_work(work);
|
||||||
}
|
}
|
||||||
PROFILE_STOP("run");
|
PROFILE_STOP("run");
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user