Two phase LBM with basic averaging
This commit is contained in:
James McClure
parent
1e8e27d530
commit
066b84acb2
2
.gitignore
vendored
2
.gitignore
vendored
|
|
@ -1,6 +1,8 @@
|
|||
*~
|
||||
*.o
|
||||
*.o*
|
||||
*.e*
|
||||
*.dat
|
||||
hybrid
|
||||
bin
|
||||
lib
|
||||
|
|
|
|||
41
Makefile
41
Makefile
|
|
@ -1,21 +1,39 @@
|
|||
CUDA_FLAGS=-arch sm_20
|
||||
CUDA_LIB=/opt/apps/cuda/5.0.35/lib64
|
||||
CUDA_FLAGS=-arch sm_35
|
||||
|
||||
bin/Color-WIA:gpu/lb2_Color_wia_mpi.o lib/libcuColor.a lib/libcuD3Q19.a lib/libcuD3Q7.a lib/libcuExtras.a
|
||||
mkdir -p bin
|
||||
CC -O3 -o bin/Color-WIA gpu/lb2_Color_wia_mpi.o -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib -Iinclude
|
||||
|
||||
|
||||
bin/Color-WIA-CBUB:gpu/lb2_Color_wia_mpi_cbub.o lib/libcuColor.a lib/libcuD3Q19.a lib/libcuD3Q7.a lib/libcuExtras.a
|
||||
mkdir -p bin
|
||||
CC -O3 -D CBUB -o bin/Color-WIA-CBUB gpu/lb2_Color_wia_mpi_cbub.o -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib -Iinclude
|
||||
|
||||
|
||||
bin/ColorLBM:gpu/lb2_Color_mpi.cpp lib/libcuColor.a lib/libcuD3Q19.a lib/libcuD3Q7.a lib/libcuExtras.a
|
||||
mkdir -p bin
|
||||
mpicxx -O3 -o bin/ColorLBM gpu/lb2_Color_mpi.cpp -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib
|
||||
|
||||
bin/gpu-Color:gpu/lb2_Color.cpp lib/libcuColor.a lib/libcuD3Q19.a lib/libcuD3Q7.a lib/libcuExtras.a
|
||||
mkdir -p bin
|
||||
g++ -O3 -o bin/gpu-Color gpu/lb2_Color.cpp -lcudart -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib -L$(CUDA_LIB)
|
||||
CC -O3 -o bin/ColorLBM gpu/lb2_Color_mpi.cpp -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib
|
||||
|
||||
#bin/gpuMRT:gpu/lb1_MRT.cu lib/libcuMRT.a lib/libcuD3Q19.a
|
||||
# mkdir -p bin
|
||||
# nvcc -O3 -o bin/gpuMRT $(CUDA_FLAGS) gpu/lb1_MRT.cu -lcuMRT -lcuD3Q19 -Llib
|
||||
|
||||
bin/gpuColor:gpu/lb2_Color.cu lib/libcuColor.a lib/libcuD3Q19.a lib/libcuExtras.a lib/libcuD3Q7.a
|
||||
mkdir -p bin
|
||||
nvcc -o bin/gpuColor $(CUDA_FLAGS) gpu/lb2_Color.cu -lcuColor -lcuD3Q19 -lcuD3Q7 -lcuExtras -Llib
|
||||
#bin/gpuColor:gpu/lb2_Color.cu lib/libcuColor.a lib/libcuD3Q19.a
|
||||
# mkdir -p bin
|
||||
# nvcc -o bin/gpuColor $(CUDA_FLAGS) gpu/lb2_Color.cu -lcuColor -lcuD3Q19 -Llib
|
||||
|
||||
gpu/lb2_Color_wia_mpi.o:gpu/lb2_Color_wia_mpi.cpp
|
||||
CC -c -o gpu/lb2_Color_wia_mpi.o gpu/lb2_Color_wia_mpi.cpp -Iinclude
|
||||
|
||||
|
||||
gpu/lb2_Color_wia_mpi_cbub.o:gpu/lb2_Color_wia_mpi.cpp
|
||||
CC -c -DCBUB -o gpu/lb2_Color_wia_mpi_cbub.o gpu/lb2_Color_wia_mpi.cpp -Iinclude
|
||||
|
||||
#pmmc/pmmc.o:pmmc/pmmc.cpp
|
||||
# CC -c -o pmmc/pmmc.o pmmc/pmmc.cpp -Iinclude
|
||||
|
||||
#include/Array.o:pmmc/Array.cpp
|
||||
# CC -c -o include/Array.o include/Array.cpp -Iinclude
|
||||
|
||||
lib/libcuExtras.a: gpu/CudaExtras.cu
|
||||
mkdir -p lib
|
||||
|
|
@ -38,5 +56,6 @@ lib/libcuColor.a: gpu/Color.cu
|
|||
nvcc -lib $(CUDA_FLAGS) gpu/Color.cu -o lib/libcuColor.a
|
||||
|
||||
clean:
|
||||
rm gpu/*.o
|
||||
rm bin/*
|
||||
rm lib/*
|
||||
rm lib/*
|
||||
|
|
|
|||
101
domain/Domain.h
101
domain/Domain.h
|
|
@ -144,6 +144,97 @@ inline void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *
|
|||
}
|
||||
}
|
||||
|
||||
inline void SignedDistance(double *Distance, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
|
||||
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
|
||||
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz)
|
||||
{
|
||||
// Use sphere lists to determine which nodes are in porespace
|
||||
// Write out binary file for nodes
|
||||
int N = Nx*Ny*Nz; // Domain size, including the halo
|
||||
double hx,hy,hz;
|
||||
double x,y,z;
|
||||
double cx,cy,cz,r;
|
||||
int imin,imax,jmin,jmax,kmin,kmax;
|
||||
int p,i,j,k,n;
|
||||
//............................................
|
||||
double min_x,min_y,min_z;
|
||||
double distance;
|
||||
//............................................
|
||||
// Lattice spacing for the entire domain
|
||||
// It should generally be true that hx=hy=hz
|
||||
// Otherwise, you will end up with ellipsoids
|
||||
hx = Lx/((Nx-2)*nprocx-1);
|
||||
hy = Ly/((Ny-2)*nprocy-1);
|
||||
hz = Lz/((Nz-2)*nprocz-1);
|
||||
//............................................
|
||||
// Get maximum and minimum for this domain
|
||||
// Halo is included !
|
||||
min_x = double(iproc*(Nx-2)-1)*hx;
|
||||
min_y = double(jproc*(Ny-2)-1)*hy;
|
||||
min_z = double(kproc*(Nz-2)-1)*hz;
|
||||
//............................................
|
||||
|
||||
//............................................
|
||||
// Pre-initialize Distance
|
||||
for (n=0;n<N;n++){
|
||||
Distance[n]=100.0;
|
||||
}
|
||||
//............................................
|
||||
|
||||
//............................................
|
||||
// .........Loop over the spheres.............
|
||||
for (p=0;p<nspheres;p++){
|
||||
// Get the sphere from the list, map to local min
|
||||
cx = List_cx[p] - min_x;
|
||||
cy = List_cy[p] - min_y;
|
||||
cz = List_cz[p] - min_z;
|
||||
r = List_rad[p];
|
||||
// Check if
|
||||
// Range for this sphere in global indexing
|
||||
imin = int ((cx-r)/hx);
|
||||
imax = int ((cx+r)/hx)+2;
|
||||
jmin = int ((cy-r)/hy);
|
||||
jmax = int ((cy+r)/hy)+2;
|
||||
kmin = int ((cz-r)/hz);
|
||||
kmax = int ((cz+r)/hz)+2;
|
||||
// Obviously we have to do something at the edges
|
||||
if (imin<0) imin = 0;
|
||||
if (imin>Nx) imin = Nx;
|
||||
if (imax<0) imax = 0;
|
||||
if (imax>Nx) imax = Nx;
|
||||
if (jmin<0) jmin = 0;
|
||||
if (jmin>Ny) jmin = Ny;
|
||||
if (jmax<0) jmax = 0;
|
||||
if (jmax>Ny) jmax = Ny;
|
||||
if (kmin<0) kmin = 0;
|
||||
if (kmin>Nz) kmin = Nz;
|
||||
if (kmax<0) kmax = 0;
|
||||
if (kmax>Nz) kmax = Nz;
|
||||
// Loop over the domain for this sphere (may be null)
|
||||
for (i=imin;i<imax;i++){
|
||||
for (j=jmin;j<jmax;j++){
|
||||
for (k=kmin;k<kmax;k++){
|
||||
// x,y,z is distance in physical units
|
||||
x = i*hx;
|
||||
y = j*hy;
|
||||
z = k*hz;
|
||||
// if inside sphere, set to zero
|
||||
// get the position in the list
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
// Compute the distance
|
||||
distance = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z)) - r;
|
||||
// Assign the minimum distance
|
||||
if (distance < Distance[n]) Distance[n] = distance;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Map the distance to lattice units
|
||||
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
|
||||
}
|
||||
|
||||
inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
|
||||
{
|
||||
//.......................................................................
|
||||
|
|
@ -242,3 +333,13 @@ inline void WriteLocalSolidID(char *FILENAME, char *ID, int N)
|
|||
File.close();
|
||||
}
|
||||
|
||||
inline void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N)
|
||||
{
|
||||
double value;
|
||||
ofstream File(FILENAME,ios::binary);
|
||||
for (int n=0; n<N; n++){
|
||||
value = Distance[n];
|
||||
File.write((char*) &value, sizeof(value));
|
||||
}
|
||||
File.close();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
136315
|
||||
8.0 8.0 8.0
|
||||
192 192 192
|
||||
10 10 10
|
||||
2 2 2
|
||||
82 82 82
|
||||
399
|
||||
1.0 1.0 1.0
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
136315
|
||||
8.0 8.0 8.0
|
||||
192 192 192
|
||||
20 20 20
|
||||
2 2 2
|
||||
80 80 80
|
||||
399
|
||||
1.0 1.0 1.0
|
||||
|
||||
|
|
|
|||
BIN
domain/ReadDomain
Executable file
BIN
domain/ReadDomain
Executable file
Binary file not shown.
106
domain/ReadDomain.cpp
Normal file
106
domain/ReadDomain.cpp
Normal file
|
|
@ -0,0 +1,106 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include "Domain.h"
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
//.......................................................................
|
||||
// Variable declaration
|
||||
//.......................................................................
|
||||
int i,j,k,n,N;
|
||||
int Nx,Ny,Nz;
|
||||
int nspheres;
|
||||
double Lx,Ly,Lz;
|
||||
// parallel domain size (# of sub-domains)
|
||||
int nprocx,nprocy,nprocz;
|
||||
// int iproc,jproc,kproc;
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
// Reading the input file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> Nx;
|
||||
domain >> Ny;
|
||||
domain >> Nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
printf("********************************************************\n");
|
||||
printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
|
||||
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
|
||||
printf("Number of spheres = %i \n", nspheres);
|
||||
printf("Domain Length (x) = %f \n", Lx);
|
||||
printf("Domain Length (y) = %f \n", Ly);
|
||||
printf("Domain Length (z) = %f \n", Lz);
|
||||
printf("********************************************************\n");
|
||||
//.......................................................................
|
||||
N = Nx*Ny*Nz;
|
||||
printf("Read input media... \n");
|
||||
|
||||
char *id;
|
||||
id = new char[N];
|
||||
|
||||
double *SignDist;
|
||||
SignDist = new double[N];
|
||||
//.......................................................................
|
||||
// Read from a sphere packing
|
||||
//.......................................................................
|
||||
double *cx,*cy,*cz,*r;
|
||||
cx = new double[nspheres];
|
||||
cy = new double[nspheres];
|
||||
cz = new double[nspheres];
|
||||
r = new double[nspheres];
|
||||
//.......................................................................
|
||||
printf("Reading the sphere packing \n");
|
||||
ReadSpherePacking(nspheres,cx,cy,cz,r);
|
||||
//.......................................................................
|
||||
|
||||
//.......................................................................
|
||||
// Write out the files
|
||||
//.......................................................................
|
||||
int rank;
|
||||
char LocalRankString[8];
|
||||
char LocalRankFilename[40];
|
||||
for (k=0; k<nprocz; k++){
|
||||
for (j=0; j<nprocy; j++){
|
||||
for (i=0; i<nprocx; i++){
|
||||
//.......................................................................
|
||||
rank = k*nprocx*nprocy + j*nprocx + i;
|
||||
//.......................................................................
|
||||
sprintf(LocalRankString,"%05d",rank);
|
||||
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
|
||||
//.......................................................................
|
||||
// printf("Assigning the local phase ID \n");
|
||||
SignedDistance(SignDist,nspheres,cx,cy,cz,r,Lx,Ly,Lz,Nx,Ny,Nz,
|
||||
i,j,k,nprocx,nprocy,nprocz);
|
||||
|
||||
for (n=0; n<N; n++){
|
||||
if (SignDist[n] < 0.0) id[n] = 0;
|
||||
else id[n] = 1;
|
||||
}
|
||||
//.......................................................................
|
||||
// printf("Generating residual NWP \n");
|
||||
GenerateResidual(id,Nx,Ny,Nz,0.3);
|
||||
//.......................................................................
|
||||
WriteLocalSolidID(LocalRankFilename, id, N);
|
||||
//.......................................................................
|
||||
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
|
||||
WriteLocalSolidDistance(LocalRankFilename, SignDist, N);
|
||||
//.......................................................................
|
||||
printf("Finished rank = %i \n",rank);
|
||||
}
|
||||
}
|
||||
}
|
||||
//.......................................................................
|
||||
|
||||
}
|
||||
|
|
@ -22,16 +22,17 @@ int main(int argc, char **argv)
|
|||
// Reading the input file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> Nx;
|
||||
domain >> Ny;
|
||||
domain >> Nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
domain >> Nx;
|
||||
domain >> Ny;
|
||||
domain >> Nz;
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
|
||||
//.......................................................................
|
||||
printf("********************************************************\n");
|
||||
printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
|
||||
|
|
|
|||
8
example/ConstrainedBubble/Color.in
Normal file
8
example/ConstrainedBubble/Color.in
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
ID
|
||||
160 32 128
|
||||
1.0
|
||||
1.0e-2 0.95 0.1 0.9
|
||||
0.7
|
||||
0.0 0.0 0.0
|
||||
0 1.0 1.0
|
||||
1000 1000 1e-5
|
||||
4
example/ConstrainedBubble/Domain.in
Normal file
4
example/ConstrainedBubble/Domain.in
Normal file
|
|
@ -0,0 +1,4 @@
|
|||
1 1 1
|
||||
100 100 100
|
||||
229
|
||||
1.0 1.0 1.0
|
||||
18
example/ConstrainedBubble/WIA.pbs
Normal file
18
example/ConstrainedBubble/WIA.pbs
Normal file
|
|
@ -0,0 +1,18 @@
|
|||
#!/bin/bash
|
||||
#PBS -A GEO019
|
||||
#PBS -N ConstrainedBubble
|
||||
#PBS -j oe
|
||||
#PBS -l walltime=2:00:00,nodes=1
|
||||
#PBS -l gres=widow2%widow3
|
||||
|
||||
#cd /tmp/work/$USER
|
||||
date
|
||||
|
||||
cd $PBS_O_WORKDIR
|
||||
|
||||
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
|
||||
source $MODULESHOME/init/bash
|
||||
module swap cray-mpich2 cray-mpich2/5.6.3
|
||||
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
|
||||
export MPICH_RDMA_ENABLED_CUDA=1
|
||||
aprun -n 1 -N 1 ~/LBPM-WIA/bin/Color-WIA-CBUB
|
||||
|
|
@ -644,8 +644,10 @@ __global__ void DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *P
|
|||
//....Instantiate the density distributions
|
||||
// Generate Equilibrium Distributions and stream
|
||||
// Stationary value - distribution 0
|
||||
Den[2*n] += 0.3333333333333333*na;
|
||||
Den[2*n+1] += 0.3333333333333333*nb;
|
||||
// Den[2*n] += 0.3333333333333333*na;
|
||||
// Den[2*n+1] += 0.3333333333333333*nb;
|
||||
atomicAdd(&Den[2*n], 0.3333333333333333*na);
|
||||
atomicAdd(&Den[2*n+1], 0.3333333333333333*nb);
|
||||
// Non-Stationary equilibrium distributions
|
||||
feq[0] = 0.1111111111111111*(1+3*ux);
|
||||
feq[1] = 0.1111111111111111*(1-3*ux);
|
||||
|
|
|
|||
2014
gpu/lb2_Color_pBC_wia_mpi.cpp
Normal file
2014
gpu/lb2_Color_pBC_wia_mpi.cpp
Normal file
File diff suppressed because it is too large
Load Diff
2113
gpu/lb2_Color_wia_mpi.cpp
Normal file
2113
gpu/lb2_Color_wia_mpi.cpp
Normal file
File diff suppressed because it is too large
Load Diff
256
include/Array.cpp
Normal file
256
include/Array.cpp
Normal file
|
|
@ -0,0 +1,256 @@
|
|||
#include "Array.h"
|
||||
/*
|
||||
* Array.cpp
|
||||
*
|
||||
* Created by James Mcclure on 3/31/09.
|
||||
* Copyright 2009 __MyCompanyName__. All rights reserved.
|
||||
*
|
||||
*/
|
||||
|
||||
// *****************************************
|
||||
// ******** class IntArray ****************
|
||||
// *****************************************
|
||||
IntArray::IntArray()
|
||||
{
|
||||
m=n=o=Length=0;
|
||||
}
|
||||
|
||||
IntArray::IntArray(int size)
|
||||
{
|
||||
data = new int [size];
|
||||
Length = size;
|
||||
m = size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
}
|
||||
|
||||
IntArray::IntArray(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
IntArray::IntArray(int nx, int ny, int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
IntArray::~IntArray()
|
||||
{
|
||||
delete data;
|
||||
}
|
||||
|
||||
|
||||
void IntArray::New(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new int [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
void IntArray::New(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
void IntArray::New(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
int IntArray::e(int i)
|
||||
{
|
||||
return data[i];
|
||||
}
|
||||
|
||||
int IntArray::e(int i, int j)
|
||||
{
|
||||
return data[m*j+i];
|
||||
}
|
||||
|
||||
int IntArray::e(int i, int j, int k)
|
||||
{
|
||||
return data[m*n*k+m*j+i];
|
||||
}
|
||||
|
||||
// *****************************************
|
||||
// ******** class DoubleArray **************
|
||||
// *****************************************
|
||||
DoubleArray::DoubleArray()
|
||||
{
|
||||
m=n=o=Length=0;
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new double [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
void DoubleArray::New(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new double [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
void DoubleArray::New(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
void DoubleArray::New(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
DoubleArray::~DoubleArray()
|
||||
{
|
||||
delete data;
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i)
|
||||
{
|
||||
return data[i];
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i, int j)
|
||||
{
|
||||
return data[j*m+i];
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i, int j, int k)
|
||||
{
|
||||
return data[k*m*n+j*m+i];
|
||||
}
|
||||
|
||||
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
|
||||
{
|
||||
if (addLength<0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
|
||||
return DoubleArray();
|
||||
}
|
||||
|
||||
int newM,newN,newO;
|
||||
if (A.o>1) {
|
||||
if (addLength%(A.m*A.n)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
|
||||
return DoubleArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n;
|
||||
newO = A.o + addLength/(A.m*A.n);
|
||||
}
|
||||
else if (A.n>1) {
|
||||
if (addLength%(A.m)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
|
||||
return DoubleArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n + addLength/A.m;
|
||||
newO = 1;
|
||||
}
|
||||
else {
|
||||
newM = A.m + addLength;
|
||||
newN = 1;
|
||||
newO = 1;
|
||||
}
|
||||
|
||||
DoubleArray toReturn(newM,newN,newO);
|
||||
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(double));
|
||||
return toReturn;
|
||||
}
|
||||
extern IntArray IncreaseSize(IntArray &A, int addLength)
|
||||
{
|
||||
if (addLength<0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
|
||||
return IntArray();
|
||||
}
|
||||
|
||||
int newM,newN,newO;
|
||||
if (A.o>1) {
|
||||
if (addLength%(A.m*A.n)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
|
||||
return IntArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n;
|
||||
newO = A.o + addLength/(A.m*A.n);
|
||||
}
|
||||
else if (A.n>1) {
|
||||
if (addLength%(A.m)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
|
||||
return IntArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n + addLength/A.m;
|
||||
newO = 1;
|
||||
}
|
||||
else {
|
||||
newM = A.m + addLength;
|
||||
newN = 1;
|
||||
newO = 1;
|
||||
}
|
||||
|
||||
IntArray toReturn(newM,newN,newO);
|
||||
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(int));
|
||||
return toReturn;
|
||||
}
|
||||
|
||||
DoubleArray DoubleArray::Copy()
|
||||
{
|
||||
DoubleArray CopyInto(m,n,o);
|
||||
// Check that the allocation worked.
|
||||
if (CopyInto.Length!=Length) return CopyInto; // Failed. Already printed an error message.
|
||||
memcpy(CopyInto.Pointer(),Pointer(),Length*sizeof(double));
|
||||
return CopyInto;
|
||||
}
|
||||
|
||||
393
include/Array.h
Normal file
393
include/Array.h
Normal file
|
|
@ -0,0 +1,393 @@
|
|||
#ifndef ARRAY_H_INC
|
||||
#define ARRAY_H_INC
|
||||
|
||||
#include <iostream>
|
||||
|
||||
// ********** ARRAY CLASS INFO **************************************
|
||||
/*
|
||||
//..............................................................
|
||||
// Overview of array classes in Array.h
|
||||
//......... Declaration of array objects........................
|
||||
// supports integer data (IntArray) or double data (DoubleArray)
|
||||
// supports one, two or three dimensional arrays
|
||||
IntArray A(m); // array of integers, Length n
|
||||
IntArray A(m,n); // array of integers, dimensions: m x n
|
||||
DoubleArray A(m,n,o); // array of doubles, dimensions m x n x o
|
||||
//............ Access the size of the array.....................
|
||||
A.m; // size of first dimension
|
||||
A.n; // size of second dimension
|
||||
A.o; // size of third dimension
|
||||
A.Length; // total number of values stored
|
||||
//........... Access the array entries .........................
|
||||
A(i); // return data[i]
|
||||
A(i,j); // return data[j*m+i]
|
||||
A(i,j,k); // return data[k*m*n+j*m+i]
|
||||
//..............................................................
|
||||
*/
|
||||
|
||||
using namespace std;
|
||||
|
||||
class IntArray{
|
||||
public:
|
||||
IntArray Copy();
|
||||
|
||||
int Length;
|
||||
int m,n,o;
|
||||
int *data;
|
||||
IntArray();
|
||||
IntArray(int size);
|
||||
IntArray(int nx,int ny);
|
||||
IntArray(int nx,int ny,int nz);
|
||||
~IntArray();
|
||||
|
||||
void New(int size);
|
||||
void New(int nx, int ny);
|
||||
void New(int nx, int ny, int nz);
|
||||
|
||||
int & operator()(int index)
|
||||
{return data[index];}
|
||||
int & operator()(int i, int j)
|
||||
{ return data[j*m+i];}
|
||||
int & operator()(int i, int j, int k)
|
||||
{ return data[k*m*n+j*m+i];}
|
||||
|
||||
int e(int i);
|
||||
int e(int i, int j);
|
||||
int e(int i, int j, int k);
|
||||
|
||||
int *Pointer() {return data;}
|
||||
};
|
||||
|
||||
class DoubleArray{
|
||||
public:
|
||||
DoubleArray Copy();
|
||||
|
||||
int Length;
|
||||
int m;
|
||||
int n;
|
||||
int o;
|
||||
double *data;
|
||||
DoubleArray();
|
||||
DoubleArray(int size);
|
||||
DoubleArray(int nx, int ny);
|
||||
DoubleArray(int nx, int ny, int nz);
|
||||
~DoubleArray();
|
||||
|
||||
void New(int size);
|
||||
void New(int nx, int ny);
|
||||
void New(int nx, int ny, int nz);
|
||||
|
||||
double & operator()(int index)
|
||||
{return data[index];}
|
||||
double & operator()(int i, int j)
|
||||
{return data[j*m+i];}
|
||||
double & operator()(int i, int j, int k)
|
||||
{return data[k*m*n+j*m+i];}
|
||||
double e(int i);
|
||||
double e(int i, int j);
|
||||
double e(int i, int j, int k);
|
||||
|
||||
double *Pointer() {return data;}
|
||||
};
|
||||
|
||||
/*
|
||||
* Array.cpp
|
||||
*
|
||||
* Created by James Mcclure on 3/31/09.
|
||||
* Copyright 2009 __MyCompanyName__. All rights reserved.
|
||||
*
|
||||
*/
|
||||
|
||||
// *****************************************
|
||||
// ******** class IntArray ****************
|
||||
// *****************************************
|
||||
/*IntArray::IntArray();
|
||||
IntArray::IntArray(int size);
|
||||
IntArray::IntArray(int nx, int ny);
|
||||
IntArray::IntArray(int nx, int ny, int nz);
|
||||
IntArray::~IntArray();
|
||||
void IntArray::New(int size);
|
||||
void IntArray::New(int nx, int ny);
|
||||
void IntArray::New(int nx, int ny,int nz);
|
||||
int IntArray::e(int i);
|
||||
int IntArray::e(int i, int j);
|
||||
int IntArray::e(int i, int j, int k);
|
||||
|
||||
// *****************************************
|
||||
// ******** class DoubleArray **************
|
||||
// *****************************************
|
||||
DoubleArray::DoubleArray();
|
||||
DoubleArray::DoubleArray(int size);
|
||||
DoubleArray::DoubleArray(int nx, int ny);
|
||||
DoubleArray::DoubleArray(int nx, int ny,int nz);
|
||||
DoubleArray::~DoubleArray();
|
||||
void DoubleArray::New(int size);
|
||||
void DoubleArray::New(int nx, int ny);
|
||||
void DoubleArray::New(int nx, int ny,int nz);
|
||||
double DoubleArray::e(int i);
|
||||
double DoubleArray::e(int i, int j);
|
||||
double DoubleArray::e(int i, int j, int k);
|
||||
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength);
|
||||
extern IntArray IncreaseSize(IntArray &A, int addLength);
|
||||
*/
|
||||
|
||||
/*
|
||||
* Array.cpp
|
||||
*
|
||||
* Created by James Mcclure on 3/31/09.
|
||||
* Copyright 2009 __MyCompanyName__. All rights reserved.
|
||||
*
|
||||
*/
|
||||
/*
|
||||
* Array.cpp
|
||||
*
|
||||
* Created by James Mcclure on 3/31/09.
|
||||
* Copyright 2009 __MyCompanyName__. All rights reserved.
|
||||
*
|
||||
*/
|
||||
|
||||
IntArray::IntArray()
|
||||
{
|
||||
m=n=o=Length=0;
|
||||
}
|
||||
|
||||
IntArray::IntArray(int size)
|
||||
{
|
||||
data = new int [size];
|
||||
Length = size;
|
||||
m = size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
}
|
||||
|
||||
IntArray::IntArray(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
IntArray::IntArray(int nx, int ny, int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
IntArray::~IntArray()
|
||||
{
|
||||
// delete [] data;
|
||||
}
|
||||
|
||||
|
||||
void IntArray::New(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new int [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
void IntArray::New(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
void IntArray::New(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new int [Length];
|
||||
}
|
||||
|
||||
int IntArray::e(int i)
|
||||
{
|
||||
return data[i];
|
||||
}
|
||||
|
||||
int IntArray::e(int i, int j)
|
||||
{
|
||||
return data[m*j+i];
|
||||
}
|
||||
|
||||
int IntArray::e(int i, int j, int k)
|
||||
{
|
||||
return data[m*n*k+m*j+i];
|
||||
}
|
||||
|
||||
// *****************************************
|
||||
// ******** class DoubleArray **************
|
||||
// *****************************************
|
||||
DoubleArray::DoubleArray()
|
||||
{
|
||||
m=n=o=Length=0;
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new double [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
DoubleArray::DoubleArray(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
void DoubleArray::New(int size)
|
||||
{
|
||||
m=size;
|
||||
n = 1;
|
||||
o = 1;
|
||||
data = new double [size];
|
||||
Length = size;
|
||||
}
|
||||
|
||||
void DoubleArray::New(int nx, int ny)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = 1;
|
||||
Length = m*n;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
void DoubleArray::New(int nx, int ny,int nz)
|
||||
{
|
||||
m = nx;
|
||||
n = ny;
|
||||
o = nz;
|
||||
Length = m*n*o;
|
||||
data = new double [Length];
|
||||
}
|
||||
|
||||
DoubleArray::~DoubleArray()
|
||||
{
|
||||
// delete [] data;
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i)
|
||||
{
|
||||
return data[i];
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i, int j)
|
||||
{
|
||||
return data[j*m+i];
|
||||
}
|
||||
|
||||
double DoubleArray::e(int i, int j, int k)
|
||||
{
|
||||
return data[k*m*n+j*m+i];
|
||||
}
|
||||
|
||||
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
|
||||
{
|
||||
if (addLength<0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
|
||||
return DoubleArray();
|
||||
}
|
||||
|
||||
int newM,newN,newO;
|
||||
if (A.o>1) {
|
||||
if (addLength%(A.m*A.n)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
|
||||
return DoubleArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n;
|
||||
newO = A.o + addLength/(A.m*A.n);
|
||||
}
|
||||
else if (A.n>1) {
|
||||
if (addLength%(A.m)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
|
||||
return DoubleArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n + addLength/A.m;
|
||||
newO = 1;
|
||||
}
|
||||
else {
|
||||
newM = A.m + addLength;
|
||||
newN = 1;
|
||||
newO = 1;
|
||||
}
|
||||
|
||||
DoubleArray toReturn(newM,newN,newO);
|
||||
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(double));
|
||||
return toReturn;
|
||||
}
|
||||
extern IntArray IncreaseSize(IntArray &A, int addLength)
|
||||
{
|
||||
if (addLength<0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
|
||||
return IntArray();
|
||||
}
|
||||
|
||||
int newM,newN,newO;
|
||||
if (A.o>1) {
|
||||
if (addLength%(A.m*A.n)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
|
||||
return IntArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n;
|
||||
newO = A.o + addLength/(A.m*A.n);
|
||||
}
|
||||
else if (A.n>1) {
|
||||
if (addLength%(A.m)!=0) {
|
||||
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
|
||||
return IntArray();
|
||||
}
|
||||
newM = A.m;
|
||||
newN = A.n + addLength/A.m;
|
||||
newO = 1;
|
||||
}
|
||||
else {
|
||||
newM = A.m + addLength;
|
||||
newN = 1;
|
||||
newO = 1;
|
||||
}
|
||||
|
||||
IntArray toReturn(newM,newN,newO);
|
||||
memcpy(toReturn.Pointer(),A.Pointer(),A.Length*sizeof(int));
|
||||
return toReturn;
|
||||
}
|
||||
|
||||
DoubleArray DoubleArray::Copy()
|
||||
{
|
||||
DoubleArray CopyInto(m,n,o);
|
||||
// Check that the allocation worked.
|
||||
if (CopyInto.Length!=Length) return CopyInto; // Failed. Already printed an error message.
|
||||
memcpy(CopyInto.Pointer(),Pointer(),Length*sizeof(double));
|
||||
return CopyInto;
|
||||
}
|
||||
|
||||
#endif
|
||||
181
include/PointList.h
Normal file
181
include/PointList.h
Normal file
|
|
@ -0,0 +1,181 @@
|
|||
|
||||
#include <math.h>
|
||||
|
||||
struct Point {
|
||||
Point() : x(0.0), y(0.0), z(0.0) {}
|
||||
Point(double xv,double yv,double zv) : x(xv), y(yv), z(zv) {}
|
||||
|
||||
double x,y,z;
|
||||
};
|
||||
|
||||
inline Point operator+(const Point &A,const Point &B) {return Point(A.x+B.x,A.y+B.y,A.z+B.z);}
|
||||
inline Point operator-(const Point &A,const Point &B) {return Point(A.x-B.x,A.y-B.y,A.z-B.z);}
|
||||
inline Point operator*(const Point &A,double v) {return Point(A.x*v,A.y*v,A.z*v);}
|
||||
inline Point operator*(double v,const Point &A) {return Point(A.x*v,A.y*v,A.z*v);}
|
||||
inline Point operator/(const Point &A,double v) {return Point(A.x/v,A.y/v,A.z/v);}
|
||||
inline Point operator-(const Point &A) {return Point(-A.x,-A.y,-A.z);}
|
||||
|
||||
inline bool operator==(const Point &A,const Point &B) {return (A.x==B.x && A.y==B.y && A.z==B.z);}
|
||||
inline bool operator!=(const Point &A,const Point &B) {return (A.x!=B.x || A.y!=B.y || A.z!=B.z);}
|
||||
|
||||
inline double Norm(const Point &A) {return sqrt(A.x*A.x+A.y*A.y+A.z*A.z);}
|
||||
inline Point Cross(const Point &A,const Point &B) {return Point(A.y*B.z-A.z*B.y,B.x*A.z-A.x*B.z,A.x*B.y-A.y*B.x);}
|
||||
inline double Dot(const Point &A,const Point &B) {return (A.x*B.x+A.y*B.y+A.z*B.z);}
|
||||
inline double Distance(const Point &A,const Point &B) {return 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));}
|
||||
|
||||
/*
|
||||
class PointList{
|
||||
public:
|
||||
int length;
|
||||
int m;
|
||||
int index;
|
||||
|
||||
Point pt;
|
||||
|
||||
double *data;
|
||||
PointList();
|
||||
PointList(int size);
|
||||
~PointList();
|
||||
|
||||
void New(int size);
|
||||
|
||||
Point & operator()(int idx)
|
||||
{
|
||||
index = idx;
|
||||
pt.x = data[3*index];
|
||||
pt.y = data[3*index+1];
|
||||
pt.z = data[3*index+2];
|
||||
return pt;
|
||||
}
|
||||
|
||||
Point &operator=(Point &P){
|
||||
pt.x = P.x;
|
||||
pt.y = P.y;
|
||||
pt.z = P.z;
|
||||
data[3*index]=pt.x;
|
||||
data[3*index+1]=pt.y;
|
||||
data[3*index+2]=pt.z;
|
||||
}
|
||||
};
|
||||
|
||||
// *****************************************
|
||||
// ******** class PointList **************
|
||||
// *****************************************
|
||||
PointList::PointList()
|
||||
{
|
||||
m=length=0;
|
||||
}
|
||||
|
||||
PointList::PointList(int size)
|
||||
{
|
||||
m=size;
|
||||
data = new double [3*size];
|
||||
length = size;
|
||||
}
|
||||
|
||||
void PointList::New(int size)
|
||||
{
|
||||
m=size;
|
||||
data = new double [3*size];
|
||||
length = size;
|
||||
}
|
||||
|
||||
|
||||
PointList::~PointList()
|
||||
{
|
||||
delete data;
|
||||
}
|
||||
*/
|
||||
template <class T>
|
||||
class DTList {
|
||||
public:
|
||||
DTList() : Data(0), length(0), refCount(new size_t(1)), outOfRange() {}
|
||||
DTList(const DTList<T> &A) : Data(A.Data), length(A.length), refCount(A.refCount), outOfRange() {++(*refCount);}
|
||||
protected:
|
||||
DTList(size_t len) : Data(len<=0 ? 0 : new T[len]), length(len<=0 ? 0 : len), refCount(new size_t(1)), outOfRange() {}
|
||||
public:
|
||||
|
||||
virtual ~DTList() {
|
||||
// --(*refCount);
|
||||
// if (*refCount==0) {delete [] Data; delete refCount;}
|
||||
// Data = 0; refCount = 0; length=0;
|
||||
}
|
||||
|
||||
DTList<T> &operator=(const DTList<T> &A) {
|
||||
if (A.refCount!=refCount) { // Otherwise doing A=A.
|
||||
--(*refCount);
|
||||
if (*refCount==0) {delete [] Data; delete refCount;}
|
||||
refCount = A.refCount;
|
||||
++(*refCount);
|
||||
length = A.length;
|
||||
Data = A.Data;
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
size_t MemoryUsed(void) const {return length*sizeof(T);}
|
||||
|
||||
const T *Pointer(void) const {return Data;}
|
||||
|
||||
size_t IsEmpty(void) const {return (Data==0);}
|
||||
size_t Length(void) const {return length;}
|
||||
|
||||
#if DTRangeCheck
|
||||
const T operator()(size_t i) const {if (i>= length) {printf("DTList: out of range"); return outOfRange;} return Data[i];}
|
||||
#else
|
||||
const T operator()(size_t i) const {return Data[i];}
|
||||
#endif
|
||||
|
||||
protected:
|
||||
T *Data;
|
||||
size_t length;
|
||||
size_t *refCount;
|
||||
|
||||
// Should be static.
|
||||
T outOfRange;
|
||||
};
|
||||
|
||||
template <class T>
|
||||
class DTMutableList : public DTList<T> {
|
||||
public:
|
||||
DTMutableList() : DTList<T>() {}
|
||||
DTMutableList(size_t len) : DTList<T>(len) {}
|
||||
DTMutableList(const DTMutableList<T> &A) : DTList<T>(A) {}
|
||||
|
||||
DTMutableList<T> &operator=(const DTMutableList<T> &A) {DTList<T>::operator=(A); return *this;}
|
||||
|
||||
T *Pointer(void) {return DTList<T>::Data;}
|
||||
const T *Pointer(void) const {return DTList<T>::Data;}
|
||||
|
||||
#if DTRangeCheck
|
||||
T &operator()(size_t i) {if (i>=DTList<T>::length) {DTErrorOutOfRange("DTList<T>",i,DTList<T>::length); return DTList<T>::outOfRange;} return DTList<T>::Data[i];}
|
||||
T operator()(size_t i) const {if (i<0 || i>= DTList<T>::length) {DTErrorOutOfRange("DTList<T>",i,DTList<T>::length); return DTList<T>::outOfRange;} return DTList<T>::Data[i];}
|
||||
#else
|
||||
T &operator()(size_t i) {return DTList<T>::Data[i];}
|
||||
T operator()(size_t i) const {return DTList<T>::Data[i];}
|
||||
#endif
|
||||
|
||||
DTMutableList<T> &operator=(T v) {for (size_t i=0;i<DTList<T>::length;i++) DTList<T>::Data[i] = v; return *this;}
|
||||
};
|
||||
|
||||
template <class T> DTMutableList<T> TruncateSize(const DTList<T> &A,size_t length)
|
||||
{
|
||||
if (length>A.Length()) length = A.Length();
|
||||
DTMutableList<T> toReturn(length);
|
||||
const T *fromP = A.Pointer();
|
||||
T *toP = toReturn.Pointer();
|
||||
for (size_t i=0;i<length;i++) toP[i] = fromP[i];
|
||||
return toReturn;
|
||||
}
|
||||
|
||||
template <class T> DTMutableList<T> IncreaseSize(const DTList<T> &A,size_t addLength)
|
||||
{
|
||||
DTMutableList<T> toReturn(A.Length()+(addLength>=0 ? addLength : 0));
|
||||
size_t len = A.Length();
|
||||
const T *fromP = A.Pointer();
|
||||
T *toP = toReturn.Pointer();
|
||||
for (size_t i=0;i<len;i++) toP[i] = fromP[i];
|
||||
return toReturn;
|
||||
}
|
||||
|
||||
|
||||
3225
include/pmmc.h
Normal file
3225
include/pmmc.h
Normal file
File diff suppressed because it is too large
Load Diff
242
pmmc/AssignBlobID.cpp
Normal file
242
pmmc/AssignBlobID.cpp
Normal file
|
|
@ -0,0 +1,242 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <math.h>
|
||||
#include <Array.h>
|
||||
|
||||
using namespace std;
|
||||
//--------------------------------------------------------------------------------------------------------
|
||||
|
||||
inline double PhaseAverageScalar(int *PhaseID, double *Scalar, int N, int Np)
|
||||
{
|
||||
int i;
|
||||
// Store the averaged values for each phase
|
||||
double *PhaseAveragedValues;
|
||||
double *PhaseVolumes;
|
||||
PhaseAveragedValues = new double[Np];
|
||||
for (i=0; i<N; i++){
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
inline int ComputeBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray indicator,
|
||||
DoubleArray F, DoubleArray S, double vf, double vs, int startx, int starty,
|
||||
int startz, IntArray temp)
|
||||
{
|
||||
// Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
// F>vf => oil phase S>vs => in porespace
|
||||
// update the list of blobs, indicator mesh
|
||||
int m = F.m; // maxima for the meshes
|
||||
int n = F.n;
|
||||
int o = F.o;
|
||||
|
||||
|
||||
int cubes_in_blob=0;
|
||||
int nrecent = 1; // number of nodes added at most recent sweep
|
||||
temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps
|
||||
temp(1,0) = starty;
|
||||
temp(2,0) = startz;
|
||||
int ntotal = 1; // total number of nodes in blob
|
||||
indicator(startx,starty,startz) = nblobs;
|
||||
|
||||
int p,s,x,y,z,start,finish,nodx,nody,nodz;
|
||||
int imin=startx,imax=startx,jmin=starty,jmax=starty; // initialize maxima / minima
|
||||
int kmin=startz,kmax=startz;
|
||||
int d[26][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
|
||||
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1},
|
||||
{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1},
|
||||
{1,1,1},{1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1},{-1,1,-1},
|
||||
{-1,-1,1},{-1,-1,-1}}; // directions to neighbors
|
||||
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
|
||||
bool status = 1; // status == true => continue to look for points
|
||||
while (status == 1){
|
||||
start = ntotal - nrecent;
|
||||
finish = ntotal;
|
||||
nrecent = 0; // set recent points back to zero for next sweep through
|
||||
for (s=start;s<finish;s++){
|
||||
// Loop over recent points; look for new points
|
||||
x = temp(0,s);
|
||||
y = temp(1,s);
|
||||
z = temp(2,s);
|
||||
// Looop over the directions
|
||||
for (p=0;p<26;p++){
|
||||
nodx=x+d[p][0];
|
||||
if (nodx < 0 ){ nodx = m-1; } // Periodic BC for x
|
||||
if (nodx > m-1 ){ nodx = 0; }
|
||||
nody=y+d[p][1];
|
||||
if (nody < 0 ){ nody = n-1; } // Periodic BC for y
|
||||
if (nody > n-1 ){ nody = 0; }
|
||||
nodz=z+d[p][2];
|
||||
if (nodz < 0 ){ nodz = 0; } // No periodic BC for z
|
||||
if (nodz > o-1 ){ nodz = o-1; }
|
||||
if ( F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) == -1 ){
|
||||
// Node is a part of the blob - add it to the list
|
||||
temp(0,ntotal) = nodx;
|
||||
temp(1,ntotal) = nody;
|
||||
temp(2,ntotal) = nodz;
|
||||
ntotal++;
|
||||
nrecent++;
|
||||
// Update the indicator map
|
||||
indicator(nodx,nody,nodz) = nblobs;
|
||||
// Update the min / max for the cube loop
|
||||
if ( nodx < imin ){ imin = nodx; }
|
||||
if ( nodx > imax ){ imax = nodx; }
|
||||
if ( nody < jmin ){ jmin = nody; }
|
||||
if ( nody > jmax ){ jmax = nody; }
|
||||
if ( nodz < kmin ){ kmin = nodz; }
|
||||
if ( nodz > kmax ){ kmax = nodz; }
|
||||
}
|
||||
else if (F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
// Some kind of error in algorithm
|
||||
printf("Error in blob search algorithm!");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
if ( nrecent == 0){
|
||||
status = 0;
|
||||
}
|
||||
}
|
||||
// Use points in temporary storage array to add cubes to the list of blobs
|
||||
if ( imin > 0) { imin = imin-1; }
|
||||
// if ( imax < m-1) { imax = imax+1; }
|
||||
if ( jmin > 0) { jmin = jmin-1; }
|
||||
// if ( jmax < n-1) { jmax = jmax+1; }
|
||||
if ( kmin > 0) { kmin = kmin-1; }
|
||||
// if ( kmax < o-1) { kmax = kmax+1; }
|
||||
int i,j,k;
|
||||
bool add;
|
||||
for (k=kmin;k<kmax;k++){
|
||||
for (j=jmin;j<jmax;j++){
|
||||
for (i=imin;i<imax;i++){
|
||||
// If cube(i,j,k) has any nodes in blob, add it to the list
|
||||
// Loop over cube edges
|
||||
add = 0;
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) == nblobs ){
|
||||
// Cube corner is in this blob
|
||||
add = 1;
|
||||
}
|
||||
}
|
||||
if (add == 1){
|
||||
// Add cube to the list
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
cubes_in_blob++;
|
||||
// Loop again to check for overlap
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if (indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
printf("Overlapping cube!");
|
||||
cout << i << ", " << j << ", " << k << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return cubes_in_blob;
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
int m,n,o;
|
||||
int p,i,j,k;
|
||||
|
||||
double vF,vS;
|
||||
/* ****************************************************************
|
||||
IDENTIFY ALL BLOBS: F > vF, S > vS
|
||||
****************************************************************** */
|
||||
// Find blob domains, number of blobs
|
||||
int nblobs = 0; // number of blobs
|
||||
int ncubes = 0; // total number of nodes in any blob
|
||||
int N = (m-1)*(n-1)*(o-1); // total number of cubes
|
||||
IntArray blobs(3,N); // store indices for blobs (cubes)
|
||||
IntArray temp(3,N); // temporary storage array
|
||||
IntArray temp2; // temporary storage array
|
||||
IntArray b(50); // number of nodes in each blob
|
||||
|
||||
IntArray indicator(m,n,o);
|
||||
DoubleArray F(m,n,o);
|
||||
DoubleArray S(m,n,o);
|
||||
|
||||
// Loop over z=0 first -> blobs attached to this end considered "connected" for LB simulation
|
||||
i=0;
|
||||
int number=0;
|
||||
for (k=0;k<1;k++){
|
||||
for (j=0;j<n;j++){
|
||||
if ( F(i,j,k) > vF ){
|
||||
if ( S(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
number = number+ComputeBlob(blobs,nblobs,ncubes,indicator,F,S,vF,vS,i,j,k,temp);
|
||||
}
|
||||
}
|
||||
}
|
||||
// Specify the blob on the z axis
|
||||
b(nblobs) = number;
|
||||
nblobs++;
|
||||
for (k=0;k<o;k++){
|
||||
for (j=0;j<n;j++){
|
||||
for (i=1;i<m;i++){
|
||||
if ( indicator(i,j,k) == -1 ){
|
||||
if ( F(i,j,k) > vF ){
|
||||
if ( S(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
b(nblobs) = ComputeBlob(blobs,nblobs,ncubes,indicator,F,S,vF,vS,i,j,k,temp);
|
||||
nblobs++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Otherwise, this point has already been assigned - ignore
|
||||
|
||||
// Make sure list blob_nodes is large enough
|
||||
if ( nblobs > b.length-1){
|
||||
b = IncreaseSize(b,b.length);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Go over all cubes again -> add any that do not contain nw phase
|
||||
bool add=1; // Set to false if any corners contain nw-phase ( F > vF)
|
||||
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;
|
||||
for (k=0;k<o-1;k++){
|
||||
for (j=0;j<n-1;j++){
|
||||
for (i=0;i<m-1;i++){
|
||||
// Loop over cube corners
|
||||
add=1; // initialize to true - add unless corner occupied by nw-phase
|
||||
for (p=0;p<8;p++){
|
||||
nodx=i+cube[p][0];
|
||||
nody=j+cube[p][1];
|
||||
nodz=k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) > -1 ){
|
||||
// corner occupied by nw-phase -> do not add
|
||||
add = 0;
|
||||
}
|
||||
}
|
||||
if ( add == 1 ){
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
count_in++;
|
||||
}
|
||||
else { count_out++; }
|
||||
}
|
||||
}
|
||||
}
|
||||
b(nblobs) = count_in;
|
||||
nblobs++;
|
||||
}
|
||||
84
pmmc/CommonCurve.cpp
Normal file
84
pmmc/CommonCurve.cpp
Normal file
|
|
@ -0,0 +1,84 @@
|
|||
//--------------------------------------------------------------------------------------------------------
|
||||
inline double ContactAngle(DTPoint3D &pt, DTVectorField3D &gradF, DTVectorField3D &gradS)
|
||||
{
|
||||
double theta;
|
||||
int m = gradS.Grid().m();
|
||||
int n = gradS.Grid().n();
|
||||
int o = gradS.Grid().o();
|
||||
|
||||
double dx = gradS.Grid().dx();
|
||||
double dy = gradS.Grid().dy();
|
||||
double dz = gradS.Grid().dz();
|
||||
|
||||
DTPoint3D origin = gradS.Grid().Origin();
|
||||
|
||||
// Get the gradient arrays
|
||||
DTFloatArray Fx = gradF.X();
|
||||
DTFloatArray Fy = gradF.Y();
|
||||
DTFloatArray Fz = gradF.Z();
|
||||
|
||||
DTFloatArray Sx = gradS.X();
|
||||
DTFloatArray Sy = gradS.Y();
|
||||
DTFloatArray Sz = gradS.Z();
|
||||
|
||||
int i,j,k;
|
||||
// Determine the cube containing this point
|
||||
// i = int(floor((pt.x-origin.x)/dx)); // real space
|
||||
// j = int(floor((pt.y-origin.y)/dy));
|
||||
// k = int(floor((pt.z-origin.z)/dz));
|
||||
i = int(floor(pt.x));
|
||||
j = int(floor(pt.y));
|
||||
k = int(floor(pt.z));
|
||||
|
||||
double x,y,z; // map to 0,1 cube
|
||||
// x = (pt.x-origin.x-i*dx)/dx; // real space
|
||||
// y = (pt.y-origin.y-j*dy)/dy;
|
||||
// z = (pt.z-origin.z-k*dz)/dz;
|
||||
x = pt.x - double(i);
|
||||
y = pt.y - double(j);
|
||||
z = pt.z - double(k);
|
||||
|
||||
// evaluate contact angle at corners of the cube
|
||||
DTMutableDoubleArray cube(2,2,2);
|
||||
// theta = acos ( -(gradF*gradS) / (|gradF| |gradS|) )
|
||||
cube(0,0,0) = -( Fx(i,j,k)*Sx(i,j,k)+Fy(i,j,k)*Sy(i,j,k)+Fz(i,j,k)*Sz(i,j,k) )
|
||||
/( sqrt(pow(Fx(i,j,k),2)+pow(Fy(i,j,k),2)+pow(Fz(i,j,k),2))
|
||||
*sqrt(pow(Sx(i,j,k),2)+pow(Sy(i,j,k),2)+pow(Sz(i,j,k),2)) );
|
||||
cube(1,0,0) = -( Fx(i+1,j,k)*Sx(i+1,j,k)+Fy(i+1,j,k)*Sy(i+1,j,k)+Fz(i+1,j,k)*Sz(i+1,j,k) )
|
||||
/( sqrt(pow(Fx(i+1,j,k),2)+pow(Fy(i+1,j,k),2)+pow(Fz(i+1,j,k),2))
|
||||
*sqrt(pow(Sx(i+1,j,k),2)+pow(Sy(i+1,j,k),2)+pow(Sz(i+1,j,k),2)) );
|
||||
cube(0,1,0) = -( Fx(i,j+1,k)*Sx(i,j+1,k)+Fy(i,j+1,k)*Sy(i,j+1,k)+Fz(i,j+1,k)*Sz(i,j+1,k) )
|
||||
/( sqrt(pow(Fx(i,j+1,k),2)+pow(Fy(i,j+1,k),2)+pow(Fz(i,j+1,k),2))
|
||||
*sqrt(pow(Sx(i,j+1,k),2)+pow(Sy(i,j+1,k),2)+pow(Sz(i,j+1,k),2)) );
|
||||
cube(0,0,1) = -( Fx(i,j,k+1)*Sx(i,j,k+1)+Fy(i,j,k+1)*Sy(i,j,k+1)+Fz(i,j,k+1)*Sz(i,j,k+1) )
|
||||
/( sqrt(pow(Fx(i,j,k+1),2)+pow(Fy(i,j,k+1),2)+pow(Fz(i,j,k+1),2))
|
||||
*sqrt(pow(Sx(i,j,k+1),2)+pow(Sy(i,j,k+1),2)+pow(Sz(i,j,k+1),2)) );
|
||||
cube(1,1,0) = -( Fx(i+1,j+1,k)*Sx(i+1,j+1,k)+Fy(i+1,j+1,k)*Sy(i+1,j+1,k)+Fz(i+1,j+1,k)*Sz(i+1,j+1,k) )
|
||||
/( sqrt(pow(Fx(i+1,j+1,k),2)+pow(Fy(i+1,j+1,k),2)+pow(Fz(i+1,j+1,k),2))
|
||||
*sqrt(pow(Sx(i+1,j+1,k),2)+pow(Sy(i+1,j+1,k),2)+pow(Sz(i+1,j+1,k),2)) );
|
||||
cube(1,0,1) = -( Fx(i+1,j,k+1)*Sx(i+1,j,k+1)+Fy(i+1,j,k+1)*Sy(i+1,j,k+1)+Fz(i+1,j,k+1)*Sz(i+1,j,k+1) )
|
||||
/( sqrt(pow(Fx(i+1,j,k+1),2)+pow(Fy(i+1,j,k+1),2)+pow(Fz(i+1,j,k+1),2))
|
||||
*sqrt(pow(Sx(i+1,j,k+1),2)+pow(Sy(i+1,j,k+1),2)+pow(Sz(i+1,j,k+1),2)) );
|
||||
cube(0,1,1) = -( Fx(i,j+1,k+1)*Sx(i,j+1,k+1)+Fy(i,j+1,k+1)*Sy(i,j+1,k+1)+Fz(i,j+1,k+1)*Sz(i,j+1,k+1) )
|
||||
/( sqrt(pow(Fx(i,j+1,k+1),2)+pow(Fy(i,j+1,k+1),2)+pow(Fz(i,j+1,k+1),2))
|
||||
*sqrt(pow(Sx(i,j+1,k+1),2)+pow(Sy(i,j+1,k+1),2)+pow(Sz(i,j+1,k+1),2)) );
|
||||
cube(1,1,1) = -( Fx(i+1,j+1,k+1)*Sx(i+1,j+1,k+1)+Fy(i+1,j+1,k+1)*Sy(i+1,j+1,k+1)+Fz(i+1,j+1,k+1)*Sz(i+1,j+1,k+1) )
|
||||
/( sqrt(pow(Fx(i+1,j+1,k+1),2)+pow(Fy(i+1,j+1,k+1),2)+pow(Fz(i+1,j+1,k+1),2))
|
||||
*sqrt(pow(Sx(i+1,j+1,k+1),2)+pow(Sy(i+1,j+1,k+1),2)+pow(Sz(i+1,j+1,k+1),2)) );
|
||||
|
||||
// construct polynomial approximation for contact angle
|
||||
double a,b,c,d,e,f,g,h; // trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
||||
a = cube(0,0,0);
|
||||
b = cube(1,0,0)-a;
|
||||
c = cube(0,1,0)-a;
|
||||
d = cube(0,0,1)-a;
|
||||
e = cube(1,1,0)-a-b-c;
|
||||
f = cube(1,0,1)-a-b-d;
|
||||
g = cube(0,1,1)-a-c-d;
|
||||
h = cube(1,1,1)-a-b-c-d-e-f-g;
|
||||
|
||||
// evaluate at x,y,z
|
||||
theta = acos (a+b*x+c*y+d*z+e*x*y+f*x*z+g*y*z+h*x*y*z );
|
||||
|
||||
return theta;
|
||||
}
|
||||
231
pmmc/blob-pmmc.cpp
Normal file
231
pmmc/blob-pmmc.cpp
Normal file
|
|
@ -0,0 +1,231 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <math.h>
|
||||
|
||||
#ifndef ARRAY_H_INC
|
||||
#include "Array.h"
|
||||
#define ARRAY_H_INC
|
||||
#endif
|
||||
|
||||
#include "PointList.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
|
||||
double d;
|
||||
// double *phase;
|
||||
// double *distance;
|
||||
|
||||
int i,j,k,p,q,r, n, Nx,Ny,Nz;
|
||||
double vS,vF;
|
||||
double readval;
|
||||
|
||||
vS = 0.f;
|
||||
vF = 0.f;
|
||||
printf("Solid surface: S(x) = %f \n",vS);
|
||||
printf("Fluid surface: F(x) = %f \n",vF);
|
||||
|
||||
Nx = Ny = Nz = 64;
|
||||
printf("Domain size is: %i x %i x %i \n",Nx,Ny,Nz);
|
||||
|
||||
|
||||
DoubleArray F(Nx,Ny,Nz);
|
||||
DoubleArray S(Nx,Ny,Nz);
|
||||
|
||||
printf("Reading distance from a file... \n");
|
||||
ifstream DISTANCE("Distance.in",ios::binary);
|
||||
for (int k=0;k<Nz;k++){
|
||||
for (int j=0;j<Ny;j++){
|
||||
for (int i=0;i<Nx;i++){
|
||||
DISTANCE.read((char *) (&readval), sizeof(readval));
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
S(i,j,k) = readval;
|
||||
}
|
||||
}
|
||||
}
|
||||
DISTANCE.close();
|
||||
|
||||
printf("Reading phase from a file... \n");
|
||||
ifstream PHASE("Phase.in",ios::binary);
|
||||
for (int k=0;k<Nz;k++){
|
||||
for (int j=0;j<Ny;j++){
|
||||
for (int i=0;i<Nx;i++){
|
||||
PHASE.read((char *) (&readval), sizeof(readval));
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
F(i,j,k) = -readval;
|
||||
//phase[n] = d;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
PHASE.close();
|
||||
|
||||
IntArray indicator(Nx,Ny,Nz); // indicates if (i,j,k) has been assigned yet
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
// Initialize indicator fucntion to -1
|
||||
indicator(i,j,k) = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
printf("Execute blob identification algorithm... \n");
|
||||
|
||||
/* ****************************************************************
|
||||
IDENTIFY ALL BLOBS: F > vF, S > vS
|
||||
****************************************************************** */
|
||||
// Find blob domains, number of blobs
|
||||
int nblobs = 0; // number of blobs
|
||||
int ncubes = 0; // total number of nodes in any blob
|
||||
int N = (Nx-1)*(Ny-1)*(Nz-1); // total number of nodes
|
||||
IntArray blobs(3,N); // store indices for blobs (cubes)
|
||||
IntArray temp(3,N); // temporary storage array
|
||||
IntArray b(50); // number of nodes in each blob
|
||||
|
||||
// Loop over z=0 first -> blobs attached to this end considered "connected" for LB simulation
|
||||
i=0;
|
||||
int number=0;
|
||||
for (k=0;k<1;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
if ( F(i,j,k) > vF ){
|
||||
if ( S(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
number = number+ComputeBlob(blobs,nblobs,ncubes,indicator,F,S,vF,vS,i,j,k,temp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Specify the blob on the z axis
|
||||
if (ncubes > 0){
|
||||
b(nblobs) = number;
|
||||
printf("Number of blobs is: %i \n",nblobs);
|
||||
nblobs++;
|
||||
}
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=1;i<Nx;i++){
|
||||
if ( indicator(i,j,k) == -1 ){
|
||||
if ( F(i,j,k) > vF ){
|
||||
if ( S(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
b(nblobs) = ComputeBlob(blobs,nblobs,ncubes,indicator,F,S,vF,vS,i,j,k,temp);
|
||||
nblobs++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Otherwise, this point has already been assigned - ignore
|
||||
|
||||
// Make sure list blob_nodes is large enough
|
||||
if ( nblobs > b.Length-1){
|
||||
printf("Increasing size of blob list \n");
|
||||
b = IncreaseSize(b,b.Length);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Go over all cubes again -> add any that do not contain nw phase
|
||||
bool add=1; // Set to false if any corners contain nw-phase ( F > vF)
|
||||
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;
|
||||
for (k=0;k<Nz-1;k++){
|
||||
for (j=0;j<Ny-1;j++){
|
||||
for (i=0;i<Nx-1;i++){
|
||||
// Loop over cube corners
|
||||
add=1; // initialize to true - add unless corner occupied by nw-phase
|
||||
for (p=0;p<8;p++){
|
||||
nodx=i+cube[p][0];
|
||||
nody=j+cube[p][1];
|
||||
nodz=k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) > -1 ){
|
||||
// corner occupied by nw-phase -> do not add
|
||||
add = 0;
|
||||
}
|
||||
}
|
||||
if ( add == 1 ){
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
count_in++;
|
||||
}
|
||||
else { count_out++; }
|
||||
}
|
||||
}
|
||||
}
|
||||
b(nblobs) = count_in;
|
||||
nblobs++;
|
||||
/* ****************************************************************
|
||||
COMPUTE GRADIENT OF F,S CURVATURE OF F
|
||||
****************************************************************** */
|
||||
|
||||
// DTVectorField3D gradF = Gradient(Fluid);
|
||||
// DTVectorField3D gradS = Gradient(Solid);
|
||||
|
||||
/* ****************************************************************
|
||||
VARIABLES FOR THE PMMC ALGORITHM
|
||||
****************************************************************** */
|
||||
double area,awn,aws,ans,lwns;
|
||||
|
||||
// Initialize arrays for return variables (areas, volumes, etc. for all blobs)
|
||||
DoubleArray nw_areas(nblobs);
|
||||
DoubleArray ns_areas(nblobs);
|
||||
DoubleArray ws_areas(nblobs);
|
||||
DoubleArray nws_length(nblobs);
|
||||
DoubleArray volume(nblobs); // Volumes of every blob
|
||||
|
||||
for (int a=0; a<nblobs; a++){
|
||||
nw_areas(a) = 0.f;
|
||||
ns_areas(a) = 0.f;
|
||||
ws_areas(a) = 0.f;
|
||||
nws_length(a) = 0.f;
|
||||
volume(a) = 0.f;
|
||||
}
|
||||
|
||||
/* ****************************************************************
|
||||
RUN PMMC ON EACH BLOB
|
||||
****************************************************************** */
|
||||
printf("Running the PMMC Algorithm \n");
|
||||
printf("The number of blobs is %i \n",nblobs);
|
||||
int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg, n_nws_seg_beg;
|
||||
int start=0,finish;
|
||||
int a,c;
|
||||
int newton_steps = 0;
|
||||
double blob_volume;
|
||||
for (a=0;a<nblobs;a++){
|
||||
|
||||
finish = start+b(a);
|
||||
|
||||
ComputeAreasPMMC(blobs, start, finish, F, S, vF, vS,
|
||||
blob_volume, ans, aws, awn, lwns, Nx, Ny, Nz);
|
||||
|
||||
start = finish;
|
||||
|
||||
volume(a) = blob_volume;
|
||||
ws_areas(a) = aws;
|
||||
nw_areas(a) = awn;
|
||||
ns_areas(a) = ans;
|
||||
|
||||
// Last "blob" is just the ws interface
|
||||
if (a+1 < nblobs){
|
||||
printf("-------------------------------- \n");
|
||||
printf("Blob id = %i \n", a);
|
||||
printf("Blob volume = %f \n", blob_volume);
|
||||
printf("Area wn = %f \n", nw_areas(a));
|
||||
printf("Area ns = %f \n", ns_areas(a));
|
||||
printf("-------------------------------- \n");
|
||||
}
|
||||
|
||||
} // End of the blob loop
|
||||
|
||||
area = 0.0;
|
||||
for (a=0;a<nblobs;a++) area+= ws_areas(a);
|
||||
printf("Area ws = %f \n", area);
|
||||
|
||||
printf("Done. \n");
|
||||
}
|
||||
439
pmmc/pmmc-areas.cpp
Normal file
439
pmmc/pmmc-areas.cpp
Normal file
|
|
@ -0,0 +1,439 @@
|
|||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "pmmc.h"
|
||||
|
||||
// Computes the interfacial areas according to the Porous Medium marching cubes (pmmc)
|
||||
// algorithm developed by McClure, Adalsteinnsson et al. (2007) for a two fluid phase
|
||||
// system. The common curve length is also computed.
|
||||
// The wn, ws, ns and solid surfaces are written to disk in the form of triangle lists
|
||||
|
||||
int main(int argc, char * argv[])
|
||||
{
|
||||
|
||||
//...........................................................................
|
||||
int Nx,Ny,Nz,N;
|
||||
Nx = Ny = Nz = N;
|
||||
int i,j,k,p,q,r,n;
|
||||
//...........................................................................
|
||||
// Interfaces are defined by the F(x,y,z) = 0.0
|
||||
// input data must be mapped such that this corresponds to the
|
||||
// desired interface
|
||||
double fluid_isovalue = 0.0;
|
||||
double solid_isovalue = 0.0;
|
||||
//...........................................................................
|
||||
|
||||
cout << "Enter the domain length (x): " << endl;
|
||||
cin >> Nx;
|
||||
cout << "Enter the domain length (y): " << endl;
|
||||
cin >> Ny;
|
||||
cout << "Enter the domain length (z): " << endl;
|
||||
cin >> Nz;
|
||||
//...........................................................................
|
||||
|
||||
//...........................................................................
|
||||
DoubleArray SignDist(Nx,Ny,Nz);
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
//...........................................................................
|
||||
|
||||
printf("Reading distance from a file... \n");
|
||||
double readval;
|
||||
ifstream DISTANCE("Distance.in",ios::binary);
|
||||
for (int k=0;k<Nz;k++){
|
||||
for (int j=0;j<Ny;j++){
|
||||
for (int i=0;i<Nx;i++){
|
||||
DISTANCE.read((char *) (&readval), sizeof(readval));
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
SignDist(i,j,k) = readval;
|
||||
}
|
||||
}
|
||||
}
|
||||
DISTANCE.close();
|
||||
|
||||
printf("Reading phase from a file... \n");
|
||||
ifstream PHASE("Phase.in",ios::binary);
|
||||
for (int k=0;k<Nz;k++){
|
||||
for (int j=0;j<Ny;j++){
|
||||
for (int i=0;i<Nx;i++){
|
||||
PHASE.read((char *) (&readval), sizeof(readval));
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Phase(i,j,k) = readval;
|
||||
//phase[n] = d;
|
||||
}
|
||||
}
|
||||
}
|
||||
PHASE.close();
|
||||
|
||||
/* ****************************************************************
|
||||
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 n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg;
|
||||
int c;
|
||||
int newton_steps = 0;
|
||||
//...........................................................................
|
||||
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||
IntArray cubeList(3,ncubes);
|
||||
int nc=0;
|
||||
//...........................................................................
|
||||
// Set up the cube list (very regular in this case due to lack of blob-ID)
|
||||
for (k=0; k<Nz-2; k++){
|
||||
for (j=0; j<Ny-2; j++){
|
||||
for (i=0; i<Nx-2; i++){
|
||||
cubeList(0,nc) = i;
|
||||
cubeList(1,nc) = j;
|
||||
cubeList(2,nc) = k;
|
||||
nc++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
n_nw_tris_beg = 0;// n_nw_tris;
|
||||
n_ns_tris_beg = 0;//n_ns_tris;
|
||||
n_ws_tris_beg = 0;//n_ws_tris;
|
||||
|
||||
// if there is a solid phase interface in the grid cell
|
||||
if (Interface(SignDist,solid_isovalue,i,j,k) == 1){
|
||||
|
||||
/////////////////////////////////////////
|
||||
/// CONSTRUCT THE LOCAL SOLID SURFACE ///
|
||||
/////////////////////////////////////////
|
||||
|
||||
// find the local solid surface
|
||||
// SOL_SURF(SignDist,0.0,Phase,fluid_isovalue,i,j,k, Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
|
||||
// local_sol_tris,n_local_sol_tris,values);
|
||||
|
||||
// find the local solid surface using the regular Marching Cubes algorithm
|
||||
SolidMarchingCubes(SignDist,0.0,Phase,fluid_isovalue,i,j,k,Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
|
||||
local_sol_tris,n_local_sol_tris,values);
|
||||
|
||||
/////////////////////////////////////////
|
||||
//////// TRIM THE SOLID SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
/* TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
|
||||
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
|
||||
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts,
|
||||
Phase, SignDist, i, j, k, newton_steps);
|
||||
*/
|
||||
TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
|
||||
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
|
||||
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts);
|
||||
|
||||
/////////////////////////////////////////
|
||||
//////// WRITE COMMON LINE POINTS ///////
|
||||
//////// TO MAIN ARRAYS ///////
|
||||
/////////////////////////////////////////
|
||||
// SORT THE LOCAL COMMON LINE POINTS
|
||||
/////////////////////////////////////////
|
||||
// Make sure the first common line point is on a face
|
||||
// Common curve points are located pairwise and must
|
||||
// be searched and rearranged accordingly
|
||||
for (p=0; p<n_local_nws_pts-1; p++){
|
||||
P = local_nws_pts(p);
|
||||
if ( P.x == 1.0*i || P.x ==1.0*(i+1)||
|
||||
P.y == 1.0*j || P.y == 1.0*(j+1) ||
|
||||
P.z == 1.0*k || P.z == 1.0*(k+1) ){
|
||||
if (p%2 == 0){
|
||||
// even points
|
||||
// Swap the pair of points
|
||||
local_nws_pts(p) = local_nws_pts(0);
|
||||
local_nws_pts(0) = P;
|
||||
P = local_nws_pts(p+1);
|
||||
local_nws_pts(p+1) = local_nws_pts(1);
|
||||
local_nws_pts(1) = P;
|
||||
p = n_local_nws_pts;
|
||||
|
||||
}
|
||||
else{
|
||||
// odd points - flip the order
|
||||
local_nws_pts(p) = local_nws_pts(p-1);
|
||||
local_nws_pts(p-1) = P;
|
||||
p-=2;
|
||||
|
||||
}
|
||||
// guarantee exit from the loop
|
||||
}
|
||||
}
|
||||
// Two common curve points per triangle
|
||||
// 0-(1=2)-(3=4)-...
|
||||
for (p=1; p<n_local_nws_pts-1; p+=2){
|
||||
A = local_nws_pts(p);
|
||||
for (q=p+1; q<n_local_nws_pts; q++){
|
||||
B = local_nws_pts(q);
|
||||
if ( A.x == B.x && A.y == B.y && A.z == B.z){
|
||||
if (q%2 == 0){
|
||||
// even points
|
||||
// Swap the pair of points
|
||||
local_nws_pts(q) = local_nws_pts(p+1);
|
||||
local_nws_pts(p+1) = B;
|
||||
B = local_nws_pts(q+1);
|
||||
local_nws_pts(q+1) = local_nws_pts(p+2);
|
||||
local_nws_pts(p+2) = B;
|
||||
q = n_local_nws_pts;
|
||||
|
||||
}
|
||||
else{
|
||||
// odd points - flip the order
|
||||
local_nws_pts(q) = local_nws_pts(q-1);
|
||||
local_nws_pts(q-1) = B;
|
||||
q-=2;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
map = n_nws_pts = 0;
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(0);
|
||||
for (p=2; p < n_local_nws_pts; p+=2){
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(p);
|
||||
|
||||
}
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(n_local_nws_pts-1);
|
||||
|
||||
for (q=0; q < n_nws_pts-1; q++){
|
||||
nws_seg(0,n_nws_seg) = map+q;
|
||||
nws_seg(1,n_nws_seg) = map+q+1;
|
||||
n_nws_seg++;
|
||||
}
|
||||
// End of the common line sorting algorithm
|
||||
/////////////////////////////////////////
|
||||
|
||||
|
||||
/////////////////////////////////////////
|
||||
////// CONSTRUCT THE nw SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
if ( n_local_nws_pts > 0){
|
||||
n_nw_tris =0;
|
||||
EDGE(Phase, fluid_isovalue, SignDist, i,j,k, Nx, Ny, Nz, nw_pts, n_nw_pts, nw_tris, n_nw_tris,
|
||||
nws_pts, n_nws_pts);
|
||||
}
|
||||
else {
|
||||
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
|
||||
}
|
||||
}
|
||||
|
||||
/////////////////////////////////////////
|
||||
////// CONSTRUCT THE nw SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
|
||||
else if (Fluid_Interface(Phase,SignDist,fluid_isovalue,i,j,k) == 1){
|
||||
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
|
||||
}
|
||||
//******END OF BLOB PMMC*********************************************
|
||||
|
||||
//*******************************************************************
|
||||
// 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_nw_tris_beg = 0;// n_nw_tris;
|
||||
// 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 \n", awn);
|
||||
printf("Area ns = %f \n", ans);
|
||||
printf("Area ws = %f \n", aws);
|
||||
printf("Area s = %f \n", As);
|
||||
printf("Length wns = %f \n", lwns);
|
||||
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);
|
||||
*/
|
||||
}
|
||||
2690
pmmc/pmmc.cpp
Normal file
2690
pmmc/pmmc.cpp
Normal file
File diff suppressed because it is too large
Load Diff
BIN
tests/cylindertest
Executable file
BIN
tests/cylindertest
Executable file
Binary file not shown.
422
tests/pmmc_cylinder.cpp
Normal file
422
tests/pmmc_cylinder.cpp
Normal file
|
|
@ -0,0 +1,422 @@
|
|||
#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 (char *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 n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg;
|
||||
int c;
|
||||
int newton_steps = 0;
|
||||
//...........................................................................
|
||||
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||
IntArray cubeList(3,ncubes);
|
||||
int nc=0;
|
||||
//...........................................................................
|
||||
// Set up the cube list (very regular in this case due to lack of blob-ID)
|
||||
for (k=0; k<Nz-2; k++){
|
||||
for (j=0; j<Ny-2; j++){
|
||||
for (i=0; i<Nx-2; i++){
|
||||
cubeList(0,nc) = i;
|
||||
cubeList(1,nc) = j;
|
||||
cubeList(2,nc) = k;
|
||||
nc++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//...........................................................................
|
||||
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;
|
||||
|
||||
n_nw_tris_beg = 0;// n_nw_tris;
|
||||
n_ns_tris_beg = 0;//n_ns_tris;
|
||||
n_ws_tris_beg = 0;//n_ws_tris;
|
||||
|
||||
// if there is a solid phase interface in the grid cell
|
||||
if (Interface(SignDist,solid_isovalue,i,j,k) == 1){
|
||||
|
||||
/////////////////////////////////////////
|
||||
/// CONSTRUCT THE LOCAL SOLID SURFACE ///
|
||||
/////////////////////////////////////////
|
||||
|
||||
// find the local solid surface
|
||||
// SOL_SURF(SignDist,0.0,Phase,fluid_isovalue,i,j,k, Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
|
||||
// local_sol_tris,n_local_sol_tris,values);
|
||||
|
||||
// find the local solid surface using the regular Marching Cubes algorithm
|
||||
SolidMarchingCubes(SignDist,0.0,Phase,fluid_isovalue,i,j,k,Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
|
||||
local_sol_tris,n_local_sol_tris,values);
|
||||
|
||||
/////////////////////////////////////////
|
||||
//////// TRIM THE SOLID SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
/* TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
|
||||
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
|
||||
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts,
|
||||
Phase, SignDist, i, j, k, newton_steps);
|
||||
*/
|
||||
TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
|
||||
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
|
||||
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts);
|
||||
|
||||
/////////////////////////////////////////
|
||||
//////// WRITE COMMON LINE POINTS ///////
|
||||
//////// TO MAIN ARRAYS ///////
|
||||
/////////////////////////////////////////
|
||||
// SORT THE LOCAL COMMON LINE POINTS
|
||||
/////////////////////////////////////////
|
||||
// Make sure the first common line point is on a face
|
||||
// Common curve points are located pairwise and must
|
||||
// be searched and rearranged accordingly
|
||||
for (p=0; p<n_local_nws_pts-1; p++){
|
||||
P = local_nws_pts(p);
|
||||
if ( P.x == 1.0*i || P.x ==1.0*(i+1)||
|
||||
P.y == 1.0*j || P.y == 1.0*(j+1) ||
|
||||
P.z == 1.0*k || P.z == 1.0*(k+1) ){
|
||||
if (p%2 == 0){
|
||||
// even points
|
||||
// Swap the pair of points
|
||||
local_nws_pts(p) = local_nws_pts(0);
|
||||
local_nws_pts(0) = P;
|
||||
P = local_nws_pts(p+1);
|
||||
local_nws_pts(p+1) = local_nws_pts(1);
|
||||
local_nws_pts(1) = P;
|
||||
p = n_local_nws_pts;
|
||||
|
||||
}
|
||||
else{
|
||||
// odd points - flip the order
|
||||
local_nws_pts(p) = local_nws_pts(p-1);
|
||||
local_nws_pts(p-1) = P;
|
||||
p-=2;
|
||||
|
||||
}
|
||||
// guarantee exit from the loop
|
||||
}
|
||||
}
|
||||
// Two common curve points per triangle
|
||||
// 0-(1=2)-(3=4)-...
|
||||
for (p=1; p<n_local_nws_pts-1; p+=2){
|
||||
A = local_nws_pts(p);
|
||||
for (q=p+1; q<n_local_nws_pts; q++){
|
||||
B = local_nws_pts(q);
|
||||
if ( A.x == B.x && A.y == B.y && A.z == B.z){
|
||||
if (q%2 == 0){
|
||||
// even points
|
||||
// Swap the pair of points
|
||||
local_nws_pts(q) = local_nws_pts(p+1);
|
||||
local_nws_pts(p+1) = B;
|
||||
B = local_nws_pts(q+1);
|
||||
local_nws_pts(q+1) = local_nws_pts(p+2);
|
||||
local_nws_pts(p+2) = B;
|
||||
q = n_local_nws_pts;
|
||||
|
||||
}
|
||||
else{
|
||||
// odd points - flip the order
|
||||
local_nws_pts(q) = local_nws_pts(q-1);
|
||||
local_nws_pts(q-1) = B;
|
||||
q-=2;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
map = n_nws_pts = 0;
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(0);
|
||||
for (p=2; p < n_local_nws_pts; p+=2){
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(p);
|
||||
|
||||
}
|
||||
nws_pts(n_nws_pts++) = local_nws_pts(n_local_nws_pts-1);
|
||||
|
||||
for (q=0; q < n_nws_pts-1; q++){
|
||||
nws_seg(0,n_nws_seg) = map+q;
|
||||
nws_seg(1,n_nws_seg) = map+q+1;
|
||||
n_nws_seg++;
|
||||
}
|
||||
// End of the common line sorting algorithm
|
||||
/////////////////////////////////////////
|
||||
|
||||
|
||||
/////////////////////////////////////////
|
||||
////// CONSTRUCT THE nw SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
if ( n_local_nws_pts > 0){
|
||||
n_nw_tris =0;
|
||||
EDGE(Phase, fluid_isovalue, SignDist, i,j,k, Nx, Ny, Nz, nw_pts, n_nw_pts, nw_tris, n_nw_tris,
|
||||
nws_pts, n_nws_pts);
|
||||
}
|
||||
else {
|
||||
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
|
||||
}
|
||||
}
|
||||
|
||||
/////////////////////////////////////////
|
||||
////// CONSTRUCT THE nw SURFACE /////////
|
||||
/////////////////////////////////////////
|
||||
|
||||
else if (Fluid_Interface(Phase,SignDist,fluid_isovalue,i,j,k) == 1){
|
||||
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
|
||||
}
|
||||
//******END OF BLOB PMMC*********************************************
|
||||
|
||||
//*******************************************************************
|
||||
// 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_nw_tris_beg = 0;// n_nw_tris;
|
||||
// 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);
|
||||
*/
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user