Adding files, uncertain status
This commit is contained in:
James E McClure
parent
9dd97d7f0f
commit
5a57449306
202
include/pmmc.h
202
include/pmmc.h
@ -3511,7 +3511,7 @@ inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList
|
|||||||
}
|
}
|
||||||
//--------------------------------------------------------------------------------------------------------
|
//--------------------------------------------------------------------------------------------------------
|
||||||
inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &CurveValues,
|
inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &CurveValues,
|
||||||
DTMutableList<Point> &Points, int npts)
|
DTMutableList<Point> &Points, int i, int j, int k, int npts)
|
||||||
{
|
{
|
||||||
int p;
|
int p;
|
||||||
Point A,B;
|
Point A,B;
|
||||||
@ -3535,7 +3535,10 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
|
|||||||
|
|
||||||
for (p=0; p<npts; p++){
|
for (p=0; p<npts; p++){
|
||||||
A = Points(p);
|
A = Points(p);
|
||||||
CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
}
|
}
|
||||||
|
|
||||||
integral = 0.0;
|
integral = 0.0;
|
||||||
@ -3543,9 +3546,11 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
|
|||||||
// Extract the line segment
|
// Extract the line segment
|
||||||
A = Points(p);
|
A = Points(p);
|
||||||
B = Points(p+1);
|
B = Points(p+1);
|
||||||
|
vA = CurveValues(p);
|
||||||
|
vB = CurveValues(p+1);
|
||||||
// Compute the length of the segment
|
// Compute the length of the segment
|
||||||
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
|
integral += 0.5*length*(vA + vB);
|
||||||
}
|
}
|
||||||
return integral;
|
return integral;
|
||||||
}
|
}
|
||||||
@ -3558,6 +3563,7 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
|
|||||||
int p;
|
int p;
|
||||||
Point A,B;
|
Point A,B;
|
||||||
double vA,vB;
|
double vA,vB;
|
||||||
|
double x,y,z;
|
||||||
double s,s1,s2,s3,temp;
|
double s,s1,s2,s3,temp;
|
||||||
double a,b,c,d,e,f,g,h;
|
double a,b,c,d,e,f,g,h;
|
||||||
double integral;
|
double integral;
|
||||||
@ -3600,9 +3606,12 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
|
|||||||
g = CubeValues(0,1,1)-a-c-d;
|
g = CubeValues(0,1,1)-a-c-d;
|
||||||
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
||||||
|
|
||||||
for (int p=0; p<npts; p++){
|
for (p=0; p<npts; p++){
|
||||||
A = Points(p);
|
A = Points(p);
|
||||||
CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
}
|
}
|
||||||
|
|
||||||
integral = 0.0;
|
integral = 0.0;
|
||||||
@ -3610,10 +3619,187 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
|
|||||||
// Extract the line segment
|
// Extract the line segment
|
||||||
A = Points(p);
|
A = Points(p);
|
||||||
B = Points(p+1);
|
B = Points(p+1);
|
||||||
|
vA = CurveValues(p);
|
||||||
|
vB = CurveValues(p+1);
|
||||||
// Compute the length of the segment
|
// Compute the length of the segment
|
||||||
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
|
integral += 0.5*length*(vA + vB);
|
||||||
}
|
}
|
||||||
|
|
||||||
return integral;
|
return integral;
|
||||||
}//--------------------------------------------------------------------------------------------------------
|
}
|
||||||
|
//--------------------------------------------------------------------------------------------------------
|
||||||
|
inline void pmmc_CubeSurfaceInterpVector(DoubleArray &Vec_x, DoubleArray &Vec_y, DoubleArray &Vec_z,
|
||||||
|
DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
|
||||||
|
DoubleArray &SurfaceVector, DoubleArray &VecAvg, int i, int j, int k, int npts, int ntris)
|
||||||
|
{
|
||||||
|
Point A,B,C;
|
||||||
|
int p;
|
||||||
|
double vA,vB,vC;
|
||||||
|
double x,y,z;
|
||||||
|
double s,s1,s2,s3,temp;
|
||||||
|
double a,b,c,d,e,f,g,h;
|
||||||
|
double integral;
|
||||||
|
|
||||||
|
// ................x component .............................
|
||||||
|
// Copy the curvature values for the cube
|
||||||
|
CubeValues(0,0,0) = Vec_x(i,j,k);
|
||||||
|
CubeValues(1,0,0) = Vec_x(i+1,j,k);
|
||||||
|
CubeValues(0,1,0) = Vec_x(i,j+1,k);
|
||||||
|
CubeValues(1,1,0) = Vec_x(i+1,j+1,k);
|
||||||
|
CubeValues(0,0,1) = Vec_x(i,j,k+1);
|
||||||
|
CubeValues(1,0,1) = Vec_x(i+1,j,k+1);
|
||||||
|
CubeValues(0,1,1) = Vec_x(i,j+1,k+1);
|
||||||
|
CubeValues(1,1,1) = Vec_x(i+1,j+1,k+1);
|
||||||
|
|
||||||
|
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
||||||
|
a = CubeValues(0,0,0);
|
||||||
|
b = CubeValues(1,0,0)-a;
|
||||||
|
c = CubeValues(0,1,0)-a;
|
||||||
|
d = CubeValues(0,0,1)-a;
|
||||||
|
e = CubeValues(1,1,0)-a-b-c;
|
||||||
|
f = CubeValues(1,0,1)-a-b-d;
|
||||||
|
g = CubeValues(0,1,1)-a-c-d;
|
||||||
|
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
||||||
|
|
||||||
|
for (p=0; p<npts; p++){
|
||||||
|
A = Points(p);
|
||||||
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
SurfaceVector(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ................y component .............................
|
||||||
|
// Copy the curvature values for the cube
|
||||||
|
CubeValues(0,0,0) = Vec_y(i,j,k);
|
||||||
|
CubeValues(1,0,0) = Vec_y(i+1,j,k);
|
||||||
|
CubeValues(0,1,0) = Vec_y(i,j+1,k);
|
||||||
|
CubeValues(1,1,0) = Vec_y(i+1,j+1,k);
|
||||||
|
CubeValues(0,0,1) = Vec_y(i,j,k+1);
|
||||||
|
CubeValues(1,0,1) = Vec_y(i+1,j,k+1);
|
||||||
|
CubeValues(0,1,1) = Vec_y(i,j+1,k+1);
|
||||||
|
CubeValues(1,1,1) = Vec_y(i+1,j+1,k+1);
|
||||||
|
|
||||||
|
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
||||||
|
a = CubeValues(0,0,0);
|
||||||
|
b = CubeValues(1,0,0)-a;
|
||||||
|
c = CubeValues(0,1,0)-a;
|
||||||
|
d = CubeValues(0,0,1)-a;
|
||||||
|
e = CubeValues(1,1,0)-a-b-c;
|
||||||
|
f = CubeValues(1,0,1)-a-b-d;
|
||||||
|
g = CubeValues(0,1,1)-a-c-d;
|
||||||
|
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
||||||
|
|
||||||
|
for (p=0; p<npts; p++){
|
||||||
|
A = Points(p);
|
||||||
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
SurfaceVector(npts+p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
|
}
|
||||||
|
|
||||||
|
// ................z component .............................
|
||||||
|
// Copy the curvature values for the cube
|
||||||
|
CubeValues(0,0,0) = Vec_z(i,j,k);
|
||||||
|
CubeValues(1,0,0) = Vec_z(i+1,j,k);
|
||||||
|
CubeValues(0,1,0) = Vec_z(i,j+1,k);
|
||||||
|
CubeValues(1,1,0) = Vec_z(i+1,j+1,k);
|
||||||
|
CubeValues(0,0,1) = Vec_z(i,j,k+1);
|
||||||
|
CubeValues(1,0,1) = Vec_z(i+1,j,k+1);
|
||||||
|
CubeValues(0,1,1) = Vec_z(i,j+1,k+1);
|
||||||
|
CubeValues(1,1,1) = Vec_z(i+1,j+1,k+1);
|
||||||
|
|
||||||
|
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
||||||
|
a = CubeValues(0,0,0);
|
||||||
|
b = CubeValues(1,0,0)-a;
|
||||||
|
c = CubeValues(0,1,0)-a;
|
||||||
|
d = CubeValues(0,0,1)-a;
|
||||||
|
e = CubeValues(1,1,0)-a-b-c;
|
||||||
|
f = CubeValues(1,0,1)-a-b-d;
|
||||||
|
g = CubeValues(0,1,1)-a-c-d;
|
||||||
|
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
||||||
|
|
||||||
|
for (p=0; p<npts; p++){
|
||||||
|
A = Points(p);
|
||||||
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
SurfaceVector(2*npts+p) = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
|
}
|
||||||
|
//.............................................................................
|
||||||
|
|
||||||
|
for (int r=0; r<ntris; r++){
|
||||||
|
A = Points(Triangles(0,r));
|
||||||
|
B = Points(Triangles(1,r));
|
||||||
|
C = Points(Triangles(2,r));
|
||||||
|
// Compute length of sides (assume dx=dy=dz)
|
||||||
|
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
|
||||||
|
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
|
||||||
|
s = 0.5*(s1+s2+s3);
|
||||||
|
temp = s*(s-s1)*(s-s2)*(s-s3);
|
||||||
|
if (temp > 0.0){
|
||||||
|
// Increment the averaged values
|
||||||
|
// x component
|
||||||
|
vA = SurfaceVector(Triangles(0,r));
|
||||||
|
vB = SurfaceVector(Triangles(1,r));
|
||||||
|
vC = SurfaceVector(Triangles(2,r));
|
||||||
|
VecAvg(0) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
|
||||||
|
// y component
|
||||||
|
vA = SurfaceVector(npts+Triangles(0,r));
|
||||||
|
vB = SurfaceVector(npts+Triangles(1,r));
|
||||||
|
vC = SurfaceVector(npts+Triangles(2,r));
|
||||||
|
VecAvg(1) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
|
||||||
|
// z component
|
||||||
|
vA = SurfaceVector(2*npts+Triangles(0,r));
|
||||||
|
vB = SurfaceVector(2*npts+Triangles(1,r));
|
||||||
|
vC = SurfaceVector(2*npts+Triangles(2,r));
|
||||||
|
VecAvg(2) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
//--------------------------------------------------------------------------------------------------------
|
||||||
|
inline double pmmc_CubeCurveInterpVector(DoubleArray &CubeValues, DoubleArray &CurveValues,
|
||||||
|
DTMutableList<Point> &Points, int i, int j, int k, int npts)
|
||||||
|
{
|
||||||
|
int p;
|
||||||
|
Point A,B;
|
||||||
|
double vA,vB;
|
||||||
|
double x,y,z;
|
||||||
|
double s,s1,s2,s3,temp;
|
||||||
|
double a,b,c,d,e,f,g,h;
|
||||||
|
double integral;
|
||||||
|
double length;
|
||||||
|
|
||||||
|
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
||||||
|
// Evaluate the coefficients
|
||||||
|
a = CubeValues(0,0,0);
|
||||||
|
b = CubeValues(1,0,0)-a;
|
||||||
|
c = CubeValues(0,1,0)-a;
|
||||||
|
d = CubeValues(0,0,1)-a;
|
||||||
|
e = CubeValues(1,1,0)-a-b-c;
|
||||||
|
f = CubeValues(1,0,1)-a-b-d;
|
||||||
|
g = CubeValues(0,1,1)-a-c-d;
|
||||||
|
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
|
||||||
|
|
||||||
|
for (p=0; p<npts; p++){
|
||||||
|
A = Points(p);
|
||||||
|
x = A.x-1.0*i;
|
||||||
|
y = A.y-1.0*j;
|
||||||
|
z = A.z-1.0*k;
|
||||||
|
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
|
||||||
|
}
|
||||||
|
|
||||||
|
integral = 0.0;
|
||||||
|
for (p=0; p < npts-1; p++){
|
||||||
|
// Extract the line segment
|
||||||
|
A = Points(p);
|
||||||
|
B = Points(p+1);
|
||||||
|
vA = CurveValues(p);
|
||||||
|
vB = CurveValues(p+1);
|
||||||
|
// Compute the length of the segment
|
||||||
|
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
integral += 0.5*length*(vA + vB);
|
||||||
|
}
|
||||||
|
return integral;
|
||||||
|
}
|
||||||
|
|||||||
283
tests/TestContactAngle.cpp
Normal file
283
tests/TestContactAngle.cpp
Normal file
@ -0,0 +1,283 @@
|
|||||||
|
#include <iostream>
|
||||||
|
#include <math.h>
|
||||||
|
#include "pmmc.h"
|
||||||
|
//#include "PointList.h"
|
||||||
|
//#include "Array.h"
|
||||||
|
|
||||||
|
#define RADIUS 15
|
||||||
|
#define HEIGHT 15.5
|
||||||
|
#define N 60
|
||||||
|
#define PI 3.14159
|
||||||
|
|
||||||
|
int main (int argc, char *argv[])
|
||||||
|
{
|
||||||
|
|
||||||
|
// printf("Radius = %s \n,"RADIUS);
|
||||||
|
int SIZE = N*N*N;
|
||||||
|
int Nx,Ny,Nz;
|
||||||
|
Nx = Ny = Nz = N;
|
||||||
|
int i,j,k,p,q,r;
|
||||||
|
|
||||||
|
// double *Solid; // cylinder
|
||||||
|
// double *Phase; // region of the cylinder
|
||||||
|
// Solid = new double [SIZE];
|
||||||
|
// Phase = new double [SIZE];
|
||||||
|
DoubleArray SignDist(Nx,Ny,Nz);
|
||||||
|
DoubleArray Phase(Nx,Ny,Nz);
|
||||||
|
double fluid_isovalue = 0.0;
|
||||||
|
double solid_isovalue = 0.0;
|
||||||
|
|
||||||
|
|
||||||
|
/* ****************************************************************
|
||||||
|
VARIABLES FOR THE PMMC ALGORITHM
|
||||||
|
****************************************************************** */
|
||||||
|
//...........................................................................
|
||||||
|
// Averaging variables
|
||||||
|
//...........................................................................
|
||||||
|
double awn,ans,aws,lwns,nwp_volume;
|
||||||
|
double As;
|
||||||
|
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
|
||||||
|
double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;
|
||||||
|
double As_global;
|
||||||
|
// bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
|
||||||
|
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||||
|
// int count_in=0,count_out=0;
|
||||||
|
// int nodx,nody,nodz;
|
||||||
|
// initialize lists for vertices for surfaces, common line
|
||||||
|
DTMutableList<Point> nw_pts(20);
|
||||||
|
DTMutableList<Point> ns_pts(20);
|
||||||
|
DTMutableList<Point> ws_pts(20);
|
||||||
|
DTMutableList<Point> nws_pts(20);
|
||||||
|
// initialize triangle lists for surfaces
|
||||||
|
IntArray nw_tris(3,20);
|
||||||
|
IntArray ns_tris(3,20);
|
||||||
|
IntArray ws_tris(3,20);
|
||||||
|
// initialize list for line segments
|
||||||
|
IntArray nws_seg(2,20);
|
||||||
|
|
||||||
|
DTMutableList<Point> tmp(20);
|
||||||
|
// IntArray store;
|
||||||
|
|
||||||
|
int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||||
|
int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||||
|
|
||||||
|
double s,s1,s2,s3; // Triangle sides (lengths)
|
||||||
|
Point A,B,C,P;
|
||||||
|
// double area;
|
||||||
|
|
||||||
|
// Initialize arrays for local solid surface
|
||||||
|
DTMutableList<Point> local_sol_pts(20);
|
||||||
|
int n_local_sol_pts = 0;
|
||||||
|
IntArray local_sol_tris(3,18);
|
||||||
|
int n_local_sol_tris;
|
||||||
|
DoubleArray values(20);
|
||||||
|
DTMutableList<Point> local_nws_pts(20);
|
||||||
|
int n_local_nws_pts;
|
||||||
|
|
||||||
|
int c;
|
||||||
|
//...........................................................................
|
||||||
|
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||||
|
IntArray cubeList(3,ncubes);
|
||||||
|
pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
|
||||||
|
//...........................................................................
|
||||||
|
double Cx,Cy,Cz;
|
||||||
|
double dist1,dist2;
|
||||||
|
Cx = Cy = Cz = N*0.51;
|
||||||
|
for (k=0; k<N; k++){
|
||||||
|
for (j=0; j<N; j++){
|
||||||
|
for (i=0; i<N; i++){
|
||||||
|
dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS;
|
||||||
|
dist2 = fabs(Cz-k)-HEIGHT;
|
||||||
|
// printf("distances = %f, %f \n",dist1,dist2);
|
||||||
|
//Solid.data[k*Nx*Ny+j*Nx+i] = dist1;
|
||||||
|
//Phase[k*Nx*Ny+j*Nx+i] = dist2;
|
||||||
|
SignDist(i,j,k) = -dist1;
|
||||||
|
Phase(i,j,k) = dist2;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
FILE *STRIS;
|
||||||
|
STRIS = fopen("solid-triangles.out","w");
|
||||||
|
|
||||||
|
FILE *WN_TRIS;
|
||||||
|
WN_TRIS = fopen("wn-tris.out","w");
|
||||||
|
|
||||||
|
FILE *NS_TRIS;
|
||||||
|
NS_TRIS = fopen("ns-tris.out","w");
|
||||||
|
|
||||||
|
FILE *WS_TRIS;
|
||||||
|
WS_TRIS = fopen("ws-tris.out","w");
|
||||||
|
|
||||||
|
FILE *WNS_PTS;
|
||||||
|
WNS_PTS = fopen("wns-pts.out","w");
|
||||||
|
|
||||||
|
// End of the loop to set the values
|
||||||
|
awn = aws = ans = lwns = 0.0;
|
||||||
|
nwp_volume = 0.0;
|
||||||
|
As = 0.0;
|
||||||
|
|
||||||
|
for (c=0;c<ncubes;c++){
|
||||||
|
// Get cube from the list
|
||||||
|
i = cubeList(0,c);
|
||||||
|
j = cubeList(1,c);
|
||||||
|
k = cubeList(2,c);
|
||||||
|
|
||||||
|
for (p=0;p<8;p++){
|
||||||
|
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0
|
||||||
|
&& SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
|
||||||
|
nwp_volume += 0.125;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Run PMMC
|
||||||
|
n_local_sol_tris = 0;
|
||||||
|
n_local_sol_pts = 0;
|
||||||
|
n_local_nws_pts = 0;
|
||||||
|
|
||||||
|
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||||
|
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||||
|
|
||||||
|
// Construct the interfaces and common curve
|
||||||
|
pmmc_ConstructLocalCube(SignDist, Phase, 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);
|
||||||
|
|
||||||
|
|
||||||
|
//*******************************************************************
|
||||||
|
// Compute the Interfacial Areas, Common Line length for blob p
|
||||||
|
// nw surface
|
||||||
|
double temp;
|
||||||
|
for (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));
|
||||||
|
// Compute length of sides (assume dx=dy=dz)
|
||||||
|
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
|
||||||
|
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
|
||||||
|
s = 0.5*(s1+s2+s3);
|
||||||
|
temp = s*(s-s1)*(s-s2)*(s-s3);
|
||||||
|
if (temp > 0.0) awn += sqrt(temp);
|
||||||
|
|
||||||
|
}
|
||||||
|
for (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));
|
||||||
|
// Compute length of sides (assume dx=dy=dz)
|
||||||
|
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
|
||||||
|
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
|
||||||
|
s = 0.5*(s1+s2+s3);
|
||||||
|
//ans=ans+sqrt(s*(s-s1)*(s-s2)*(s-s3));
|
||||||
|
temp = s*(s-s1)*(s-s2)*(s-s3);
|
||||||
|
if (temp > 0.0) ans += sqrt(temp);
|
||||||
|
}
|
||||||
|
for (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));
|
||||||
|
// Compute length of sides (assume dx=dy=dz)
|
||||||
|
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
|
||||||
|
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
|
||||||
|
s = 0.5*(s1+s2+s3);
|
||||||
|
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
|
||||||
|
temp = s*(s-s1)*(s-s2)*(s-s3);
|
||||||
|
if (temp > 0.0) aws += sqrt(temp);
|
||||||
|
}
|
||||||
|
for (r=0;r<n_local_sol_tris;r++){
|
||||||
|
A = local_sol_pts(local_sol_tris(0,r));
|
||||||
|
B = local_sol_pts(local_sol_tris(1,r));
|
||||||
|
C = local_sol_pts(local_sol_tris(2,r));
|
||||||
|
// Compute length of sides (assume dx=dy=dz)
|
||||||
|
s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
|
||||||
|
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
|
||||||
|
s = 0.5*(s1+s2+s3);
|
||||||
|
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
|
||||||
|
temp = s*(s-s1)*(s-s2)*(s-s3);
|
||||||
|
if (temp > 0.0) As += sqrt(temp);
|
||||||
|
}
|
||||||
|
for (p=0; p < n_local_nws_pts-1; p++){
|
||||||
|
// Extract the line segment
|
||||||
|
A = local_nws_pts(p);
|
||||||
|
B = local_nws_pts(p+1);
|
||||||
|
// Compute the length of the segment
|
||||||
|
s = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
|
||||||
|
// Add the length to the common line
|
||||||
|
lwns += s;
|
||||||
|
}
|
||||||
|
//.......................................................................................
|
||||||
|
// Write the triangle lists to text file
|
||||||
|
for (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));
|
||||||
|
fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||||
|
}
|
||||||
|
for (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));
|
||||||
|
fprintf(WS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||||
|
}
|
||||||
|
for (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));
|
||||||
|
fprintf(NS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||||
|
}
|
||||||
|
for (p=0; p < n_nws_pts; p++){
|
||||||
|
P = nws_pts(p);
|
||||||
|
fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
|
||||||
|
}
|
||||||
|
|
||||||
|
//.......................................................................................
|
||||||
|
|
||||||
|
for (r=0;r<n_local_sol_tris;r++){
|
||||||
|
A = local_sol_pts(local_sol_tris(0,r));
|
||||||
|
B = local_sol_pts(local_sol_tris(1,r));
|
||||||
|
C = local_sol_pts(local_sol_tris(2,r));
|
||||||
|
fprintf(STRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
|
||||||
|
}
|
||||||
|
//*******************************************************************
|
||||||
|
// Reset the triangle counts to zero
|
||||||
|
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||||
|
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||||
|
// n_ns_tris_beg = 0;//n_ns_tris;
|
||||||
|
// n_ws_tris_beg = 0;//n_ws_tris;
|
||||||
|
// n_nws_seg_beg = n_nws_seg;
|
||||||
|
//*******************************************************************
|
||||||
|
}
|
||||||
|
fclose(WN_TRIS);
|
||||||
|
fclose(NS_TRIS);
|
||||||
|
fclose(WS_TRIS);
|
||||||
|
fclose(WNS_PTS);
|
||||||
|
fclose(STRIS);
|
||||||
|
|
||||||
|
printf("-------------------------------- \n");
|
||||||
|
printf("NWP volume = %f \n", nwp_volume);
|
||||||
|
printf("Area wn = %f, Analytical = %f \n", awn,2*PI*RADIUS*RADIUS);
|
||||||
|
printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
|
||||||
|
printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
|
||||||
|
printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
|
||||||
|
printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
|
||||||
|
printf("-------------------------------- \n");
|
||||||
|
//.........................................................................
|
||||||
|
|
||||||
|
/* FILE *PHASE;
|
||||||
|
PHASE = fopen("Phase.in","wb");
|
||||||
|
fwrite(Phase,8,SIZE,PHASE);
|
||||||
|
fclose(PHASE);
|
||||||
|
|
||||||
|
FILE *SOLID;
|
||||||
|
SOLID = fopen("Distance.in","wb");
|
||||||
|
fwrite(Solid,8,SIZE,SOLID);
|
||||||
|
fclose(SOLID);
|
||||||
|
*/
|
||||||
|
}
|
||||||
@ -145,7 +145,6 @@ int main (int argc, char *argv[])
|
|||||||
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_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,
|
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);
|
i, j, k, Nx, Ny, Nz);
|
||||||
|
|
||||||
|
|
||||||
//*******************************************************************
|
//*******************************************************************
|
||||||
// Compute the Interfacial Areas, Common Line length for blob p
|
// Compute the Interfacial Areas, Common Line length for blob p
|
||||||
|
|||||||
@ -141,7 +141,7 @@ int main(int argc, char **argv)
|
|||||||
|
|
||||||
printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
|
printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
|
||||||
|
|
||||||
FILE *CURVATURE;
|
/* FILE *CURVATURE;
|
||||||
CURVATURE = fopen("Curvature.dat","wb");
|
CURVATURE = fopen("Curvature.dat","wb");
|
||||||
fwrite(MeanCurvature.data,8,N,CURVATURE);
|
fwrite(MeanCurvature.data,8,N,CURVATURE);
|
||||||
fclose(CURVATURE);
|
fclose(CURVATURE);
|
||||||
@ -150,5 +150,5 @@ int main(int argc, char **argv)
|
|||||||
DISTANCE = fopen("SignDist.dat","wb");
|
DISTANCE = fopen("SignDist.dat","wb");
|
||||||
fwrite(Phase.data,8,N,DISTANCE);
|
fwrite(Phase.data,8,N,DISTANCE);
|
||||||
fclose(DISTANCE);
|
fclose(DISTANCE);
|
||||||
|
*/
|
||||||
}
|
}
|
||||||
|
|||||||
145
tests/TestSphereOrientation.cpp
Normal file
145
tests/TestSphereOrientation.cpp
Normal file
@ -0,0 +1,145 @@
|
|||||||
|
#include <iostream>
|
||||||
|
#include <math.h>
|
||||||
|
#include "pmmc.h"
|
||||||
|
#include "Domain.h"
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
/*
|
||||||
|
* Compare the measured and analytical curvature for a sphere
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
|
||||||
|
int main(int argc, char **argv)
|
||||||
|
{
|
||||||
|
int i,j,k;
|
||||||
|
int Nx,Ny,Nz,N;
|
||||||
|
double Lx,Ly,Lz;
|
||||||
|
double fluid_isovalue=0.0;
|
||||||
|
double solid_isovalue=0.0;
|
||||||
|
|
||||||
|
Lx = Ly = Lz = 1.0;
|
||||||
|
Nx = Ny = Nz = 102;
|
||||||
|
N = Nx*Ny*Nz;
|
||||||
|
|
||||||
|
//...........................................................................
|
||||||
|
// Set up the cube list
|
||||||
|
//...........................................................................
|
||||||
|
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||||
|
IntArray cubeList(3,ncubes);
|
||||||
|
pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
|
||||||
|
//...........................................................................
|
||||||
|
|
||||||
|
//****************************************************************************************
|
||||||
|
// Create the structures needed to carry out the PMMC algorithm
|
||||||
|
//****************************************************************************************
|
||||||
|
DTMutableList<Point> nw_pts(20);
|
||||||
|
DTMutableList<Point> ns_pts(20);
|
||||||
|
DTMutableList<Point> ws_pts(20);
|
||||||
|
DTMutableList<Point> nws_pts(20);
|
||||||
|
// initialize triangle lists for surfaces
|
||||||
|
IntArray nw_tris(3,20);
|
||||||
|
IntArray ns_tris(3,20);
|
||||||
|
IntArray ws_tris(3,20);
|
||||||
|
// initialize list for line segments
|
||||||
|
IntArray nws_seg(2,20);
|
||||||
|
|
||||||
|
DTMutableList<Point> tmp(20);
|
||||||
|
|
||||||
|
DoubleArray wn_curvature(20);
|
||||||
|
|
||||||
|
int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||||
|
int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||||
|
|
||||||
|
// Initialize arrays for local solid surface
|
||||||
|
DTMutableList<Point> local_sol_pts(20);
|
||||||
|
int n_local_sol_pts = 0;
|
||||||
|
IntArray local_sol_tris(3,18);
|
||||||
|
int n_local_sol_tris;
|
||||||
|
DoubleArray values(20);
|
||||||
|
DTMutableList<Point> local_nws_pts(20);
|
||||||
|
int n_local_nws_pts;
|
||||||
|
//****************************************************************************************
|
||||||
|
|
||||||
|
int nspheres=1;
|
||||||
|
// Set up the signed distance function for a sphere
|
||||||
|
double *cx,*cy,*cz,*rad;
|
||||||
|
cx = new double[nspheres];
|
||||||
|
cy = new double[nspheres];
|
||||||
|
cz = new double[nspheres];
|
||||||
|
rad = new double[nspheres];
|
||||||
|
|
||||||
|
// Sphere in the middle of the domain
|
||||||
|
cx[0] = cy[0] = cz[0] = 0.5;
|
||||||
|
rad[0] = 0.3;
|
||||||
|
|
||||||
|
DoubleArray SignDist(Nx,Ny,Nz);
|
||||||
|
DoubleArray Phase(Nx,Ny,Nz);
|
||||||
|
DoubleArray Phase_x(Nx,Ny,Nz);
|
||||||
|
DoubleArray Phase_y(Nx,Ny,Nz);
|
||||||
|
DoubleArray Phase_z(Nx,Ny,Nz);
|
||||||
|
DoubleArray CubeValues(2,2,2);
|
||||||
|
|
||||||
|
// Compute the signed distance function
|
||||||
|
SignedDistance(Phase.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
|
||||||
|
|
||||||
|
for (k=0; k<Nz; k++){
|
||||||
|
for (j=0; j<Ny; j++){
|
||||||
|
for (i=0; i<Nx; i++){
|
||||||
|
SignDist(i,j,k) = 100.0;
|
||||||
|
Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SignedDistance(SignDist.data,0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
|
||||||
|
|
||||||
|
double Gxx_sum,Gyy_sum,Gzz_sum,Gxy_sum,Gxz_sum,Gyz_sum;
|
||||||
|
double wn_curvature_sum = 0.0;
|
||||||
|
double wn_area_sum = 0.0;
|
||||||
|
|
||||||
|
pmmc_MeshGradient(Phase, Phase_x, Phase_y, Phase_z, Nx, Ny, Nz);
|
||||||
|
|
||||||
|
for (int c=0;c<ncubes;c++){
|
||||||
|
|
||||||
|
// Get cube from the list
|
||||||
|
i = cubeList(0,c);
|
||||||
|
j = cubeList(1,c);
|
||||||
|
k = cubeList(2,c);
|
||||||
|
|
||||||
|
// Run PMMC
|
||||||
|
n_local_sol_tris = 0;
|
||||||
|
n_local_sol_pts = 0;
|
||||||
|
n_local_nws_pts = 0;
|
||||||
|
|
||||||
|
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
|
||||||
|
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
|
||||||
|
|
||||||
|
// Construct the interfaces and common curve
|
||||||
|
pmmc_ConstructLocalCube(SignDist, Phase, 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);
|
||||||
|
|
||||||
|
// Copy the curvature values for the cube
|
||||||
|
CubeValues(0,0,0) = Phase_x(i,j,k);
|
||||||
|
CubeValues(1,0,0) = Phase_x(i+1,j,k);
|
||||||
|
CubeValues(0,1,0) = Phase_x(i,j+1,k);
|
||||||
|
CubeValues(1,1,0) = Phase_x(i+1,j+1,k);
|
||||||
|
CubeValues(0,0,1) = Phase_x(i,j,k+1);
|
||||||
|
CubeValues(1,0,1) = Phase_x(i+1,j,k+1);
|
||||||
|
CubeValues(0,1,1) = Phase_x(i,j+1,k+1);
|
||||||
|
CubeValues(1,1,1) = Phase_x(i+1,j+1,k+1);
|
||||||
|
|
||||||
|
// Interpolate the curvature onto the surface
|
||||||
|
// wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, nw_pts, nw_tris,
|
||||||
|
// wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
|
||||||
|
|
||||||
|
wn_area_sum += pmmc_CubeSurfaceArea(nw_pts, nw_tris, n_nw_tris);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
printf("Gxx = %f, Analytical = 1/3 \n", wn_curvature_sum/wn_area_sum);
|
||||||
|
|
||||||
|
}
|
||||||
Loading…
Reference in New Issue
Block a user