Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA
This commit is contained in:
James McClure
@@ -58,9 +58,9 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue)
|
||||
//int Nx = Field.size(0);
|
||||
//int Ny = Field.size(1);
|
||||
//int Nz = Field.size(2);
|
||||
for (int k=0; k<Nz-1; k++){
|
||||
for (int j=0; j<Ny-1; j++){
|
||||
for (int i=0; i<Nx-1; i++){
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
object.LocalIsosurface(Field,isovalue,i,j,k);
|
||||
for (int idx=0; idx<object.TriangleCount; idx++){
|
||||
e1 = object.Face(idx);
|
||||
@@ -137,6 +137,7 @@ void Minkowski::MeasureObject(){
|
||||
}
|
||||
}
|
||||
CalcDist(distance,id,*Dm);
|
||||
Eikonal(distance, id, *Dm, 10, {true, true, true});
|
||||
ComputeScalar(distance,0.0);
|
||||
|
||||
}
|
||||
|
||||
@@ -1,7 +1,5 @@
|
||||
#include "analysis/dcel.h"
|
||||
|
||||
|
||||
|
||||
DECL::DECL(){
|
||||
}
|
||||
|
||||
@@ -15,6 +13,25 @@ int DECL::Face(int index){
|
||||
return FaceData[index];
|
||||
}
|
||||
|
||||
void DECL::Write(){
|
||||
int e1,e2,e3;
|
||||
FILE *TRIANGLES;
|
||||
TRIANGLES = fopen("triangles.stl","w");
|
||||
fprintf(TRIANGLES,"solid \n");
|
||||
for (int idx=0; idx<TriangleCount; idx++){
|
||||
e1 = Face(idx);
|
||||
e2 = halfedge.next(e1);
|
||||
e3 = halfedge.next(e2);
|
||||
auto P1 = vertex.coords(halfedge.v1(e1));
|
||||
auto P2 = vertex.coords(halfedge.v1(e2));
|
||||
auto P3 = vertex.coords(halfedge.v1(e3));
|
||||
fprintf(TRIANGLES,"vertex %f %f %f\n",P1.x,P1.y,P1.z);
|
||||
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);
|
||||
}
|
||||
fclose(TRIANGLES);
|
||||
}
|
||||
|
||||
void DECL::LocalIsosurface(const DoubleArray& A, double value, const int i, const int j, const int k){
|
||||
Point P,Q;
|
||||
Point PlaceHolder;
|
||||
@@ -350,243 +367,43 @@ double DECL::EdgeAngle(int edge)
|
||||
return angle;
|
||||
}
|
||||
|
||||
void Isosurface(DoubleArray &A, const double &v)
|
||||
void iso_surface(const Array<double>&Field, const double isovalue)
|
||||
{
|
||||
NULL_USE( v );
|
||||
|
||||
Point P,Q;
|
||||
Point PlaceHolder;
|
||||
Point C0,C1,C2,C3,C4,C5,C6,C7;
|
||||
|
||||
int TriangleCount;
|
||||
int VertexCount;
|
||||
int CubeIndex;
|
||||
|
||||
Point VertexList[12];
|
||||
Point NewVertexList[12];
|
||||
int LocalRemap[12];
|
||||
|
||||
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];
|
||||
|
||||
int Nx = A.size(0);
|
||||
int Ny = A.size(1);
|
||||
int Nz = A.size(2);
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
DECL object;
|
||||
int e1,e2,e3;
|
||||
FILE *TRIANGLES;
|
||||
TRIANGLES = fopen("isosurface.stl","w");
|
||||
fprintf(TRIANGLES,"solid isosurface\n");
|
||||
int Nx = Field.size(0);
|
||||
int Ny = Field.size(1);
|
||||
int Nz = Field.size(2);
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
object.LocalIsosurface(Field,isovalue,i,j,k);
|
||||
for (int idx=0; idx<object.TriangleCount; idx++){
|
||||
e1 = object.Face(idx);
|
||||
e2 = object.halfedge.next(e1);
|
||||
e3 = object.halfedge.next(e2);
|
||||
auto P1 = object.vertex.coords(object.halfedge.v1(e1));
|
||||
auto P2 = object.vertex.coords(object.halfedge.v1(e2));
|
||||
auto P3 = object.vertex.coords(object.halfedge.v1(e3));
|
||||
auto Normal = object.TriNormal(e1);
|
||||
// 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;
|
||||
//P3.x += 1.0*i; P3.y += 1.0*j; P3.z +=1.0*k;
|
||||
fprintf(TRIANGLES,"facet normal %f %f %f\n",Normal.x,Normal.y,Normal.z);
|
||||
fprintf(TRIANGLES," outer loop\n");
|
||||
fprintf(TRIANGLES," vertex %f %f %f\n",P1.x,P1.y,P1.z);
|
||||
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);
|
||||
fprintf(TRIANGLES," endloop\n");
|
||||
fprintf(TRIANGLES,"endfacet\n");
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
fprintf(TRIANGLES,"endsolid isosurface\n");
|
||||
fclose(TRIANGLES);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,3 +1,6 @@
|
||||
#ifndef DCEL_INC
|
||||
#define DCEL_INC
|
||||
|
||||
#include <vector>
|
||||
#include "analysis/pmmc.h"
|
||||
|
||||
@@ -67,6 +70,7 @@ public:
|
||||
Vertex vertex;
|
||||
Halfedge halfedge;
|
||||
void LocalIsosurface(const DoubleArray& A, double value, int i, int j, int k);
|
||||
void Write();
|
||||
int Face(int index);
|
||||
|
||||
double origin(int edge);
|
||||
@@ -78,3 +82,7 @@ public:
|
||||
private:
|
||||
std::vector<int> FaceData;
|
||||
};
|
||||
|
||||
void iso_surface(const Array<double>&Field, const double isovalue);
|
||||
|
||||
#endif
|
||||
|
||||
@@ -182,6 +182,148 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
|
||||
}
|
||||
}
|
||||
|
||||
double Eikonal(DoubleArray &Distance, const Array<char> &ID, Domain &Dm, int timesteps, const std::array<bool,3>& periodic){
|
||||
|
||||
/*
|
||||
* This routine converts the data in the Distance array to a signed distance
|
||||
* by solving the equation df/dt = sign(1-|grad f|), where Distance provides
|
||||
* the values of f on the mesh associated with domain Dm
|
||||
* It has been tested with segmented data initialized to values [-1,1]
|
||||
* and will converge toward the signed distance to the surface bounding the associated phases
|
||||
*
|
||||
* Reference:
|
||||
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics229
|
||||
*/
|
||||
|
||||
int i,j,k;
|
||||
double dt=0.1;
|
||||
double Dx,Dy,Dz;
|
||||
double Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
|
||||
double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
|
||||
double sign,norm;
|
||||
double LocalVar,GlobalVar,LocalMax,GlobalMax;
|
||||
|
||||
int xdim,ydim,zdim;
|
||||
xdim=Dm.Nx-2;
|
||||
ydim=Dm.Ny-2;
|
||||
zdim=Dm.Nz-2;
|
||||
//fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
|
||||
fillHalo<double> fillData( Dm.Comm, Dm.rank_info, {xdim, ydim, zdim}, {1,1,1}, 50, 1, {true,true,true}, periodic );
|
||||
|
||||
// Arrays to store the second derivatives
|
||||
DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
|
||||
DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
|
||||
DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
|
||||
|
||||
int count = 0;
|
||||
while (count < timesteps){
|
||||
|
||||
// Communicate the halo of values
|
||||
fillData.fill(Distance);
|
||||
|
||||
// Compute second order derivatives
|
||||
for (k=1;k<Dm.Nz-1;k++){
|
||||
for (j=1;j<Dm.Ny-1;j++){
|
||||
for (i=1;i<Dm.Nx-1;i++){
|
||||
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
|
||||
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
|
||||
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
|
||||
}
|
||||
}
|
||||
}
|
||||
fillData.fill(Dxx);
|
||||
fillData.fill(Dyy);
|
||||
fillData.fill(Dzz);
|
||||
|
||||
LocalMax=LocalVar=0.0;
|
||||
// Execute the next timestep
|
||||
for (k=1;k<Dm.Nz-1;k++){
|
||||
for (j=1;j<Dm.Ny-1;j++){
|
||||
for (i=1;i<Dm.Nx-1;i++){
|
||||
|
||||
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
|
||||
|
||||
sign = -1;
|
||||
if (ID(i,j,k) == 1) sign = 1;
|
||||
|
||||
// local second derivative terms
|
||||
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
|
||||
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
|
||||
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
|
||||
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
|
||||
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
|
||||
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
|
||||
|
||||
/* //............Compute upwind derivatives ...................
|
||||
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
|
||||
Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
|
||||
Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
|
||||
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
|
||||
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
|
||||
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
|
||||
*/
|
||||
Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
|
||||
Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
|
||||
Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
|
||||
|
||||
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
|
||||
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
|
||||
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
|
||||
|
||||
// Compute upwind derivatives for Godunov Hamiltonian
|
||||
if (sign < 0.0){
|
||||
if (Dxp + Dxm > 0.f) Dx = Dxp*Dxp;
|
||||
else Dx = Dxm*Dxm;
|
||||
|
||||
if (Dyp + Dym > 0.f) Dy = Dyp*Dyp;
|
||||
else Dy = Dym*Dym;
|
||||
|
||||
if (Dzp + Dzm > 0.f) Dz = Dzp*Dzp;
|
||||
else Dz = Dzm*Dzm;
|
||||
}
|
||||
else{
|
||||
|
||||
if (Dxp + Dxm < 0.f) Dx = Dxp*Dxp;
|
||||
else Dx = Dxm*Dxm;
|
||||
|
||||
if (Dyp + Dym < 0.f) Dy = Dyp*Dyp;
|
||||
else Dy = Dym*Dym;
|
||||
|
||||
if (Dzp + Dzm < 0.f) Dz = Dzp*Dzp;
|
||||
else Dz = Dzm*Dzm;
|
||||
}
|
||||
|
||||
//Dx = max(Dxp*Dxp,Dxm*Dxm);
|
||||
//Dy = max(Dyp*Dyp,Dym*Dym);
|
||||
//Dz = max(Dzp*Dzp,Dzm*Dzm);
|
||||
|
||||
norm=sqrt(Dx + Dy + Dz);
|
||||
if (norm > 1.0) norm=1.0;
|
||||
|
||||
Distance(i,j,k) += dt*sign*(1.0 - norm);
|
||||
LocalVar += dt*sign*(1.0 - norm);
|
||||
|
||||
if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
|
||||
LocalMax = fabs(dt*sign*(1.0 - norm));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
|
||||
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm);
|
||||
GlobalVar /= Dm.Volume;
|
||||
count++;
|
||||
|
||||
if (count%50 == 0 && Dm.rank()==0 )
|
||||
printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
|
||||
|
||||
if (fabs(GlobalMax) < 1e-5){
|
||||
if (Dm.rank()==0) printf("Exiting with max tolerance of 1e-5 \n");
|
||||
count=timesteps;
|
||||
}
|
||||
}
|
||||
return GlobalVar;
|
||||
}
|
||||
|
||||
// Explicit instantiations
|
||||
template void CalcDist<float>( Array<float>&, const Array<char>&, const Domain&, const std::array<bool,3>&, const std::array<double,3>& );
|
||||
|
||||
@@ -16,6 +16,16 @@ struct Vec {
|
||||
};
|
||||
inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.z < r.x*r.x+r.y*r.y+r.z*r.z; }
|
||||
|
||||
inline double minmod(double &a, double &b){
|
||||
|
||||
double value;
|
||||
|
||||
value = a;
|
||||
if ( a*b < 0.0) value=0.0;
|
||||
else if (fabs(a) > fabs(b)) value = b;
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
/*!
|
||||
* @brief Calculate the distance using a simple method
|
||||
@@ -40,4 +50,16 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
|
||||
void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
|
||||
const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Calculate the distance based on solution of Eikonal equation
|
||||
* @details This routine calculates the signed distance to the nearest domain surface.
|
||||
* @param[out] Distance Distance function
|
||||
* @param[in] ID Domain id
|
||||
* @param[in] Dm Domain information
|
||||
* @param[in] timesteps number of timesteps to run for Eikonal solver
|
||||
* @param[in] periodic Directions that are periodic
|
||||
*/
|
||||
double Eikonal(DoubleArray &Distance, const Array<char> &ID, Domain &Dm, int timesteps, const std::array<bool,3>& periodic);
|
||||
|
||||
#endif
|
||||
|
||||
@@ -232,7 +232,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
|
||||
if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
|
||||
// Fill remaining variables
|
||||
N = Nx*Ny*Nz;
|
||||
Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
|
||||
Volume = nx*ny*nz*nproc[0]*nproc[1]*nproc[2]*1.0;
|
||||
|
||||
if (myrank==0) printf("voxel length = %f micron \n", voxel_length);
|
||||
|
||||
@@ -266,8 +266,8 @@ void Domain::Decomp( const std::string& Filename )
|
||||
outlet_layers_x = 0;
|
||||
outlet_layers_y = 0;
|
||||
outlet_layers_z = 0;
|
||||
inlet_layers_phase=1;
|
||||
outlet_layers_phase=2;
|
||||
inlet_layers_phase=1;
|
||||
outlet_layers_phase=2;
|
||||
checkerSize = 32;
|
||||
|
||||
// Read domain parameters
|
||||
|
||||
Reference in New Issue
Block a user