Форум

The code is shown below:
It's working perfect,

What I'm trying to do is code what I attached in the excel sheet into this one below:

// Upscaling-FMM.cpp : Defines the entry point for the console application.
// 28jan-2012.cpp : Defines the entry point for the console application.
// 27jan2012.cpp : Defines the entry point for the console application.
//
// 11-dec-noon.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <math.h>
#include <time.h>
#include <ctime>
#include <fstream>
#include <iomanip>
//#include <vector>
//#include <algorithm>
#include <stdio.h>
#include <stdlib.h>
using namespace std;

//=========================================================================================================
double *dxc;
double *dyc;
double *dzc;
double *dx;
double *dy;
double *dz;
//=========================================================================================================

//#define eps 2.2204460492503131e-16
#define eps 3.e-16
#define doublemax 1e50
#define INF 2e50
//#define INF 111
#define listINF 2.345e50
//#define listINF 111
#ifndef min
#define min(a,b) ((a) < (b) ? (a): (b))
#endif
#ifndef max
#define max(a,b) ((a) > (b) ? (a): (b))
#endif

/* Find minimum value of an array and return its index */
int minarray(double *A, int l) {
int i;
int minind=0;
for (i=0; i<l; i++) { if(A<A[minind]){ minind=i; } }
return minind;
}

double pow2(double val) { return val*val; }

int iszero(double a){ return a*a<eps; }
int isnotzero(double a){ return a*a>eps; }

void roots(double* Coeff, double* ans) {
double a=Coeff[0];
double b=Coeff[1];
double c=Coeff[2];
double r1, r2;
double d;

d=max(pow2(b)-4.0*a*c,0);
if(isnotzero(a)) {
ans[0]= (-b - sqrt(d)) / (2.0*a);
ans[1]= (-b + sqrt(d)) / (2.0*a);
}
else {
r1=(-b - sqrt(d));
r2=(-b + sqrt(d));
if(isnotzero(r1))
{
if(isnotzero(r2))
{
ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r2;
}
else
{
ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r1;
}
}
else if(isnotzero(r2))
{
ans[0]= (2.0*c)/r2; ans[1]= (2.0*c)/r2;
}
else
{
ans[0]=0; ans[1]=0;
}
}
}

int maxarray(double *A, int l) {
int i;
int maxind=0;
for (i=0; i<l; i++) { if(A>A[maxind]){ maxind=i; } }
return maxind;
}

__inline int mindex3(int x, int y, int z, int sizx, int sizy) { return x+y*sizx+z*sizx*sizy; }

__inline bool IsFinite(double x) { return (x <= doublemax && x >= -doublemax ); }
__inline bool IsInf(double x) { return (x >= doublemax ); }

__inline bool IsListInf(double x){ return (x == listINF ); }

__inline bool isntfrozen3d(int i, int j, int k, int *dims, bool *Frozen) {

return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==0);
}
__inline bool isfrozen3d(int i, int j, int k, int *dims, bool *Frozen) {

return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==1);
}

int p2x(int x) /* 2^x */
{
/* return pow(2,x); */
int y=1;
int p2x[16] ={1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768};
while(x>15) { x=x-15; y=y*32768; }
return y*p2x[x];
}

void show_list(double **listval, int *listprop) {
int z, k;
double mm;

for(z=0;z<listprop[1]; z++) {
for(k=0;k<p2x(z+1); k++) {
mm=k;
if((z>0)&&(listval[z-1][(int)floor(mm/2)]>listval[z][k])) {
printf("*%6.2f", listval[z][k]);
}
else {
printf(" %6.2f", listval[z][k]);
}
}
printf("\n");
}
}

void initialize_list(double ** listval, int *listprop) {
/* Loop variables */
int i;
/* Current Length, Orde and Current Max Length */
listprop[0]=0; listprop[1]=1; listprop[2]=2;
/* Make first orde storage of 2 values */
listval[0]=new double[2];
/* Initialize on infinite */
for(i=0;i<2;i++) { listval[0]=listINF; }
}

void destroy_list(double ** listval, int *listprop) {
/* Loop variables */
int i, list_orde;
/* Get list orde */
list_orde=listprop[1];
/* Free memory */
for(i=0;i<list_orde;i++) { free(listval); }
free(listval);
free(listprop);
}

void list_add(double ** listval, int *listprop, double val) {
/* List parameters */
int list_length, list_orde, list_lengthmax;
/* Loop variable */
int i, j;
/* Temporary list location */
int listp;
/* Get list parameters */
list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2];
/* If list is full expand list */
if(list_length==list_lengthmax) {
list_lengthmax=list_lengthmax*2;
for (i=list_orde; i>0; i--) {
listval=listval[i-1];
listval = (double *)realloc(listval, p2x(i+1)*sizeof(double));
for(j=p2x(i); j<p2x(i+1); j++) { listval[j]=listINF; }
}
listval[0]=new double[2];
listval[0][0]=min(listval[1][0], listval[1][1]);
listval[0][1]=listINF;
list_orde++;
}
/* Add a value to the list */
listp=list_length;
list_length++;
listval[list_orde-1][listp]=val;
/* Update the links minimum */
for (i=list_orde-1; i>0; i--) {
listp=(int)floor(((double)listp)/2);
if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; }
}
/* Set list parameters */
listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax;
}

int list_minimum(double ** listval, int *listprop) {
/* List parameters */
int list_length, list_orde, list_lengthmax;
/* Loop variable */
int i;
/* Temporary list location */
int listp;
/* Index of Minimum */
int minindex;
/* Get list parameters */
list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2];
/* Follow the minimum through the binary tree */
listp=0;
for(i=0;i<(list_orde-1);i++) {

/* <= ????? */
if(listval[listp]<=listval[listp+1]) { listp=listp*2; } else { listp=(listp+1)*2; }
}
i=list_orde-1;
if(listval[listp]<=listval[listp+1]){minindex=listp; } else { minindex=listp+1; }
return minindex;
}
void list_remove(double ** listval, int *listprop, int index) {
/* List parameters */
int list_length, list_orde, list_lengthmax;
/* Loop variable */
int i;
/* Temp index */
int index2;
double val;
double valmin;
/* Get list parameters */
list_length=listprop[0];
list_orde=listprop[1];
list_lengthmax=listprop[2];
/* Temporary store current value */
val=listval[list_orde-1][index];
valmin=listINF;
/* Replace value by infinite */
listval[list_orde-1][index]=listINF;
/* Follow the binary tree to replace value by minimum values from */
/** the other values in the binary tree. */

i=list_orde-1;
while(true) {
if((index%2)==0) { index2=index+1; } else { index2=index-1; }
if(val<listval[index2]) {
index=(int)floor(((double)index2)/2.0);
if(listval[index2]<valmin) { valmin=listval[index2]; }
if(i==0) { break; }
listval[i-1][index]=valmin;
i--; if(i==0) { break; }
}
else { break; }
}
}

void list_remove_replace(double ** listval, int *listprop, int index) {
/* List parameters */
int list_length, list_orde, list_lengthmax;
/* Loop variable */
int i, listp;
/* Temporary store value */
double val;
int templ;
/* Get list parameters */
list_length=listprop[0];
list_orde=listprop[1];
list_lengthmax=listprop[2];
/* Remove the value the minimum number and goup in tree by comparision its adjacent number to the parents */
list_remove(listval, listprop, index);
/* Replace the removed value by the last in the list. (if it was */
/* not already the last value) */
if(index<(list_length-1)) {
/* Temporary store last value in the list */
val=listval[list_orde-1][list_length-1];
/* Remove last value in the list */
list_remove(listval, listprop, list_length-1);
/* Add a value to the list */
listp=index;
listval[list_orde-1][index]=val;
/* Update the links minimum */
for (i=(list_orde-1); i>0; i--) {
listp=(int)floor(((double)listp)/2);
if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; }
}
}
/* List is now shorter */
list_length--;
/* Remove trunk of binary tree / Free memory if list becomes shorter */
if(list_orde>2&&IsListInf(listval[0][1])) {
//cout<<"truncate-start"<<endl;
//show_list(listval, listprop);
list_orde--;
list_lengthmax=list_lengthmax/2;
/* Remove trunk array */
free(listval[0]);
/* Move the other arrays one up */
templ=2;
for (i=0; i<list_orde; i++) {
listval=listval[i+1];
/* Resize arrays to their new shorter size */
listval = (double *)realloc(listval, templ*sizeof(double));
templ*=2;
}

}
/* Set list parameters */
listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax;
//show_list(listval, listprop);
//cout<<"truncate-end"<<endl;
}

void listupdate(double **listval, int *listprop, int index, double val) {
/* List parameters */
int list_length, list_orde, list_lengthmax;
/* loop variable */
int i, listp;
/* Get list parameters */
list_length=listprop[0];
list_orde=listprop[1];
list_lengthmax=listprop[2];
/* Add a value to the list */
listp=index;
listval[list_orde-1][index]=val;
/* Update the links minimum */
for (i=(list_orde-1); i>0; i--) {
listp=(int)floor(((double)listp)/2);
if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; }
}
/* Set list parameters */
listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax;
}

__inline int mindex2(int x, int y, int sizx) { return x+y*sizx; }

__inline bool isntfrozen2d(int i, int j, int *dims, bool *Frozen)
{
return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==0);
}
__inline bool isfrozen2d(int i, int j, int *dims, bool *Frozen)
{
return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==1);
}
//=========================================================================================================
double second_derivative(double Txm1, double Txm2, double Txp1, double Txp2) {
bool ch1, ch2;
double Tm;
Tm=INF;
ch1=(Txm2<Txm1)&&IsFinite(Txm1); ch2=(Txp2<Txp1)&&IsFinite(Txp1);
if(ch1&&ch2) { Tm =min( (4.0*Txm1-Txm2)/3.0 , (4.0*Txp1-Txp2)/3.0);}
else if(ch1) { Tm =(4.0*Txm1-Txm2)/3.0; }
else if(ch2) { Tm =(4.0*Txp1-Txp2)/3.0; }
return Tm;
}
double CalculateDistance(double *T, double Fijk, int *dims, int i, int j, int k, bool *Frozen) {

/* Loop variables */
int q, t;

/* Current location */
int in, jn, kn;

/* Derivatives */
double Tm[3]={0, 0, 0};
double Coeff[3];

/* local derivatives in distance image */
double Txm1, Txm2, Txp1, Txp2;
double Tym1, Tym2, Typ1, Typ2;
double Tzm1, Tzm2, Tzp1, Tzp2;
double Tt, Tt2,xm1,xp1,ym1,yp1,zm1,zp1,x1,x2,y1,y2,z1,z2;

/* Return values root of polynomial */
double ansroot[2]={0, 0};

/* Order derivatives in a certain direction */
int Order[3]={0, 0, 0};

/* Neighbours 4x2 */
int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1};

/* Stencil constants */
double G1[3]={1, 1, 1};

x1 = dx[mindex3(i, j, k, dims[0], dims[1])];
y1 = dy[mindex3(i, j, k, dims[0], dims[1])];
z1 = dz[mindex3(i, j, k, dims[0], dims[1])];

/*Get First order derivatives (only use frozen pixel) */
in=i-1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txm1=T[mindex3(in, jn, kn, dims[0], dims[1])];xm1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txm1=INF; xm1=0; }
in=i+1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txp1=T[mindex3(in, jn, kn, dims[0], dims[1])];xp1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txp1=INF; xp1=0; }
in=i+0; jn=j-1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tym1=T[mindex3(in, jn, kn, dims[0], dims[1])];ym1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tym1=INF; ym1=0;}
in=i+0; jn=j+1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Typ1=T[mindex3(in, jn, kn, dims[0], dims[1])];yp1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Typ1=INF; yp1=0;}
in=i+0; jn=j+0; kn=k-1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzm1=T[mindex3(in, jn, kn, dims[0], dims[1])];zm1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzm1=INF; zm1=0;}
in=i+0; jn=j+0; kn=k+1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzp1=T[mindex3(in, jn, kn, dims[0], dims[1])];zp1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzp1=INF; zp1=0;}

/*Make 1e order derivatives in x and y direction */
Tm[0] = min( Txm1 , Txp1); if(IsFinite(Tm[0])){ Order[0]=1; } else { Order[0]=0; }
Tm[1] = min( Tym1 , Typ1); if(IsFinite(Tm[1])){ Order[1]=1; } else { Order[1]=0; }
Tm[2] = min( Tzm1 , Tzp1); if(IsFinite(Tm[2])){ Order[2]=1; } else { Order[2]=0; }
// X2=XM1 IF CONDITIONAL EXPRESSION IS TRUE AND XP1 IF IT IS FALSE
x2 = Tm[0] == Txm1 ? xm1 : xp1;
y2 = Tm[1] == Tym1 ? ym1 : yp1;
z2 = Tm[2] == Tzm1 ? zm1 : zp1;

G1[0] = abs(x2 + x1)/2;
G1[1] = abs(y2 + y1)/2;
G1[2] = abs(z2 + z1)/2;

/*Calculate the distance using x and y direction */
Coeff[0]=0; Coeff[1]=0; Coeff[2]=-1/(max(pow2(Fijk),eps));

for (t=0; t<3; t++) {
switch(Order[t]) {
case 1:
Coeff[0]+=G1[t]; Coeff[1]+=-2.0*Tm[t]*G1[t]; Coeff[2]+=pow2(Tm[t])*G1[t];
break;
}
}

roots(Coeff, ansroot);
Tt=max(ansroot[0], ansroot[1]);
return Tt;
}
//starting your modificatiob

//==============================================================================
int main()
{
clock_t startClock,finishClock,currentClock;
vector<int>::iterator pos;
double timeCount,temp_time;
startClock = clock();

//-----your sort function goes here here--------
// Declare size of two dimensional array and initialize.
int nx,ny,nz,nt,i,j,k,initial_zero,index,inf;
cout << endl << "enter nx == " ;
cin>> nx;

cout << endl << "enter ny == " ;
cin>> ny;

cout << endl << "enter nz == " ;
cin>> nz;

nt=nx*ny*nz;

startClock = clock();
//ofstream timefile;
cout<<"t1= ";
cout<<startClock/1000<<endl;

int s,x,y,z,XYZ_index,w,IJK_index,itt,npixels,q;
double Tt;
double *F, *SourcePoints;
int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1};
/* Neighbour list */
int neg_free;
int neg_pos;
double *neg_listv;
double *neg_listx;
double *neg_listy;
double *neg_listz;
double *neg_listo;
/* Matrix containing the Frozen Pixels" */
bool *Frozen;
/* Size of input image */
int *dims_c;
int dims[3];
/* Size of SourcePoints array */
int num_source;

cout << endl << "enter num_source == " ;
cin>> num_source;

// assign the dimensions
dims[0]=nx;
dims[1]=ny;
dims[2]=nz;

npixels=dims[0]*dims[1]*dims[2];

//double* dx = NULL;
dx = new double[nt];

//double* dy = NULL;
dy = new double[nt];

//double* dz = NULL;
dz = new double[nt];

for (int ii=0;ii<npixels;ii++){
dx[ii]=1;
dy[ii]=1;
dz[ii]=1;
}

/* The output distance image */
double *T;
int *listprop;
double **listval;
double *Final_T;
double *perm;
double *permx;
double *permx_orig;
double *time_orig;
double *dtc;
double *porosity;
double *porosity_orig;
double *source;

permx_orig=new double[npixels];
time_orig=new double[npixels];
dtc=new double[npixels];
porosity_orig=new double[npixels];
Final_T=new double[npixels];
perm=new double[npixels];
permx=new double[npixels];
porosity=new double[npixels];

source=new double[num_source*num_source*3];

int kk;
kk=-1;

for(int i=0;i<num_source;i++){
for(int j=0;j<num_source;j++){
k=0;
kk=kk+1;

source[3*kk]=nx/2+1;
source[3*kk+1]=ny/2+1;
source[3*kk+2]=nz/2+1;
cout<<source[3*(kk)]<<"-"<<source[3*kk+1]<<"-"<<source[3*kk+2]<<endl;

}
}

for (int ii=0;ii<npixels;ii++){
permx[ii]=1;
permx_orig[ii]=1;
porosity[ii]=1;
porosity_orig[ii]=1;
}
int ii;

T= Final_T;
F=permx;
SourcePoints=source;

/* Pixels which are processed and have a final distance are frozen */
Frozen = new bool[npixels];
for(q=0;q<npixels;q++){Frozen[q]=0; T[q]=-1;}

/*Free memory to store neighbours of the (segmented) region */
neg_free = 1000000;
neg_pos=0;

neg_listx = new double [neg_free];
neg_listy = new double [neg_free];
neg_listz = new double [neg_free];

/* List parameters array */
listprop=new int [3];
/* Make jagged list to store a maximum of 2^64 values */
//listval= (double **)malloc( 64* sizeof(double *) );/// num of rows=64?
listval= new double* [64];
/* Initialize parameter list */
initialize_list(listval, listprop);
neg_listv=listval[listprop[1]-1];

for (s=0; s<num_source*num_source; s++) {
/*starting point */
x= (int)SourcePoints[0+s*3]-1;
y= (int)SourcePoints[1+s*3]-1;
z= (int)SourcePoints[2+s*3]-1;
XYZ_index=mindex3(x, y, z, dims[0], dims[1]);

Frozen[XYZ_index]=1;
T[XYZ_index]=0;
}

for (s=0; s<num_source*num_source; s++) {
/*starting point */
x= (int)SourcePoints[0+s*3]-1;
y= (int)SourcePoints[1+s*3]-1;

z= (int)SourcePoints[2+s*3]-1;

XYZ_index=mindex3(x, y, z, dims[0], dims[1]);
for (w=0; w<6; w++) {
/*Location of neighbour */
i=x+ne[w];
j=y+ne[w+6];
k=z+ne[w+12];

IJK_index=mindex3(i, j, k, dims[0], dims[1]);

/*Check if current neighbour is not yet frozen and inside the */

/*picture */
if(isntfrozen3d(i, j, k, dims, Frozen)) {
Tt=(1/(max(F[IJK_index],eps)));
/*Update distance in neigbour list or add to neigbour list */
if(T[IJK_index]>0) {
if(neg_listv[(int)T[IJK_index]]>Tt) {
listupdate(listval, listprop, (int)T[IJK_index], Tt);
}
}
else {
/*If running out of memory at a new block */
//if(neg_pos>=neg_free) {
// neg_free+=100000;
// neg_listx = (double *)realloc(neg_listx, neg_free*sizeof(double) );
// neg_listy = (double *)realloc(neg_listy, neg_free*sizeof(double) );
// neg_listz = (double *)realloc(neg_listz, neg_free*sizeof(double) );
//}
list_add(listval, listprop, Tt);
neg_listv=listval[listprop[1]-1];
neg_listx[neg_pos]=i;
neg_listy[neg_pos]=j;
neg_listz[neg_pos]=k;
T[IJK_index]=neg_pos;
neg_pos++;
}
}
}
}

/*Loop through all pixels of the image */
for (itt=0; itt<(npixels); itt++) /* */ {
/*Get the pixel from narrow list (boundary list) with smallest */
/*distance value and set it to current pixel location */
index=list_minimum(listval, listprop);
neg_listv=listval[listprop[1]-1];
/* Stop if pixel distance is infinite (all pixels are processed) */
if(IsInf(neg_listv[index])) { break; }

/*index=minarray(neg_listv, neg_pos); */
x=(int)neg_listx[index]; y=(int)neg_listy[index]; z=(int)neg_listz[index];
XYZ_index=mindex3(x, y, z, dims[0], dims[1]);
Frozen[XYZ_index]=1;
T[XYZ_index]=neg_listv[index];
//cout<< T[XYZ_index]<<endl;
//show_list(listval, listprop);
/*Remove min value by replacing it with the last value in the array */
//1- remove the min value and go up in tree by comparing it with its adjacent bothers and parent
//2- remove the last value in the list and again go up in tree
//3- replace the minimum by the last value and again go up in the tree
//- if we have extra tree delete them
list_remove_replace(listval, listprop, index) ;
//show_list(listval, listprop);
neg_listv=listval[listprop[1]-1];
if(index<(neg_pos-1)) {
neg_listx[index]=neg_listx[neg_pos-1];
neg_listy[index]=neg_listy[neg_pos-1];
neg_listz[index]=neg_listz[neg_pos-1];
T[(int)mindex3((int)neg_listx[index], (int)neg_listy[index], (int)neg_listz[index], dims[0], dims[1])]=index;
}
neg_pos =neg_pos-1;

/*Loop through all 6 neighbours of current pixel */
for (w=0;w<6;w++) {
/*Location of neighbour */
i=x+ne[w]; j=y+ne[w+6]; k=z+ne[w+12];
IJK_index=mindex3(i, j, k, dims[0], dims[1]);

/*Check if current neighbour is not yet frozen and inside the */
/*picture */
if(isntfrozen3d(i, j, k, dims, Frozen)) {
//Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, usesecond, usecross, Frozen);
Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, Frozen);

/*Update distance in neigbour list or add to neigbour list */
IJK_index=mindex3(i, j, k, dims[0], dims[1]);
if((T[IJK_index]>-1)&&T[IJK_index]<=listprop[0]) {
if(neg_listv[(int)T[IJK_index]]>Tt) {
listupdate(listval, listprop, (int)T[IJK_index], Tt);
}
}
else {
list_add(listval, listprop, Tt);
//show_list(listval, listprop);
neg_listv=listval[listprop[1]-1];
neg_listx[neg_pos]=i; neg_listy[neg_pos]=j; neg_listz[neg_pos]=k;

T[IJK_index]=neg_pos;
neg_pos++;
}
}
}
}
/* Free memory */
/* Destroy parameter list */
destroy_list(listval, listprop);
free(neg_listx);
free(neg_listy);
free(neg_listz);
free(Frozen);
finishClock=clock();
timeCount = finishClock - startClock ;
cout<<"run time = ";
cout<<timeCount/1000/60<<endl;

//==================================================================================
ofstream myfile1;
myfile1.open ("fine-time.txt");
//myfile << "Writing this to a file.\n";

double temp=0;
myfile1 << setw(16) ;
for (int k=1; k<nz+1;k++){
for (int j=1; j<ny+1;j++){
for (int i=1; i<nx+1;i++){
index= (k-1) * nx * ny +(j-1) * (nx) + (i-1) ;
/*if(Final_T[index]==inf){
Final_T[index]=1.;
}*/
myfile1 <<Final_T[index]<< setw(16) ;
// save fine scale time
time_orig[index]= Final_T[index];
//myfile1 <<F[index]<< setw(16) ;
temp=temp+Final_T[index];

}

myfile1 << endl ;
}
}
myfile1.close();
cout<<"Summation-time = ";
cout<<temp/1000;

//=====================================================================================
}