Revision 80d469e116108f7e0dada3ef59563953b9676f8c authored by Eh Tan on 09 December 2006, 00:27:47 UTC, committed by Eh Tan on 09 December 2006, 00:27:47 UTC
1 parent 0b9c4d2
Tracer_advection.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
* Clint Conrad, Michael Gurnis, and Eun-seo Choi.
* Copyright (C) 1994-2005, California Institute of Technology.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/*
Tracer_advection.c
A program which initiates the distribution of tracers
and advects those tracers in a time evolving velocity field.
Called and used from the CitCOM finite element code.
Written 2/96 M. Gurnis for Citcom in cartesian geometry
Modified by Lijie in 1998 and by Vlad and Eh in 2005 for CitcomS
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <mpi.h>
#include <math.h>
#include <sys/types.h>
#ifdef HAVE_MALLOC_H
#include <malloc.h>
#endif
#include <stdlib.h> /* for "system" command */
#include "element_definitions.h"
#include "global_defs.h"
#include "parsing.h"
void tracer_input(struct All_variables *E)
{
int m=E->parallel.me;
input_int("tracer",&(E->control.tracer),"0",m);
input_string("tracer_file",E->control.tracer_file,"tracer.dat",m);
}
void tracer_initial_settings(E)
struct All_variables *E;
{
void tracer_setup();
void tracer_advection();
E->problem_tracer_setup=tracer_setup;
E->problem_tracer_advection=tracer_advection;
}
void tracer_setup(E)
struct All_variables *E;
{
int i,j,k,node;
int m,ntr;
float idummy,xdummy,ydummy,zdummy;
FILE *fp;
int nox,noy,noz,gnox,gnoy,gnoz;
char output_file[255];
MPI_Comm world;
MPI_Status status;
gnox=E->mesh.nox;
gnoy=E->mesh.noy;
gnoz=E->mesh.noz;
nox=E->lmesh.nox;
noy=E->lmesh.noy;
noz=E->lmesh.noz;
sprintf(output_file,"%s",E->control.tracer_file);
fp=fopen(output_file,"r");
if (fp == NULL) {
fprintf(E->fp,"(Tracer_advection #1) Cannot open %s\n", output_file);
exit(8);
}
fscanf(fp,"%d",&(E->Tracer.NUM_TRACERS));
E->Tracer.tracer_x=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
E->Tracer.tracer_y=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
E->Tracer.tracer_z=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
E->Tracer.itcolor=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
E->Tracer.x_space=(float*) malloc((gnox+1)*sizeof(float));
E->Tracer.y_space=(float*) malloc((gnoy+1)*sizeof(float));
E->Tracer.z_space=(float*) malloc((gnoz+1)*sizeof(float));
/***comment by Vlad 03/15/2005
each processor holds its own number of tracers
***/
ntr=1;
for(i=1;i<=E->Tracer.NUM_TRACERS;i++) {
fscanf(fp,"%f %f %f %f", &idummy, &xdummy, &ydummy, &zdummy);
if(xdummy >= E->sx[1][1][1] && xdummy <= E->sx[1][1][nox*noy*noz]) {
if(ydummy >= E->sx[1][2][1] && ydummy <= E->sx[1][2][nox*noy*noz]) {
if(zdummy >= E->sx[1][3][1] && zdummy <= E->sx[1][3][nox*noy*noz]) {
E->Tracer.itcolor[ntr]=idummy;
E->Tracer.tracer_x[ntr]=xdummy;
E->Tracer.tracer_y[ntr]=ydummy;
E->Tracer.tracer_z[ntr]=zdummy;
ntr++;
}
}
}
}
/***comment by Vlad 3/30/2005
E->Tracer.LOCAL_NUM_TRACERS is the initial number
of tracers in each processor
***/
E->Tracer.LOCAL_NUM_TRACERS=ntr-1;
/***comment by Vlad 1/26/2005
reading the local mesh coordinate
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(i=1;i<=nox;i++)
{
j=1;
k=1;
node=k+(i-1)*noz+(j-1)*nox*noz;
E->Tracer.x_space[i]=E->sx[m][1][node];
}
for(j=1;j<=noy;j++)
{
i=1;
k=1;
node=k+(i-1)*noz+(j-1)*nox*noz;
E->Tracer.y_space[j]=E->sx[m][2][node];
}
for(k=1;k<=noz;k++)
{
i=1;
j=1;
node=k+(i-1)*noz+(j-1)*nox*noz;
E->Tracer.z_space[k]=E->sx[m][3][node];
}
} /* end of m */
return;
}
void tracer_advection(E)
struct All_variables *E;
{
int i,j,k,l,m,n,o,p;
int n_x,n_y,n_z;
int node1,node2,node3,node4,node5,node6,node7,node8;
int nno,nox,noy,noz,gnox,gnoy,gnoz;
int iteration;
float x_tmp, y_tmp, z_tmp;
float volume, tr_dx, tr_dy, tr_dz, dx, dy, dz;
float w1,w2,w3,w4,w5,w6,w7,w8;
float tr_v[NCS][4];
MPI_Comm world;
MPI_Status status[4];
MPI_Status status_count;
MPI_Request request[4];
float xmin,xmax,ymin,ymax,zmin,zmax;
float *tr_color_1,*tr_x_1,*tr_y_1,*tr_z_1;
float *tr_color_new[13],*tr_x_new[13],*tr_y_new[13],*tr_z_new[13];
int *Left_proc_list,*Right_proc_list;
int *jump_new,*count_new;
int *jj;
int proc;
int Previous_num_tracers,Current_num_tracers;
int locx,locy,locz;
int left,right,up,down,back,front;
world=E->parallel.world;
gnox=E->mesh.nox;
gnoy=E->mesh.noy;
gnoz=E->mesh.noz;
nox=E->lmesh.nox;
noy=E->lmesh.noy;
noz=E->lmesh.noz;
nno=nox*noy*noz;
Left_proc_list=(int*) malloc(6*sizeof(int));
Right_proc_list=(int*) malloc(6*sizeof(int));
jump_new=(int*) malloc(6*sizeof(int));
count_new=(int*) malloc(6*sizeof(int));
jj=(int*) malloc(6*sizeof(int));
tr_x_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_y_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_z_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_color_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
for(i=0;i<=11;i++){
tr_color_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_x_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_y_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
tr_z_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
}
/*** comment by Vlad 3/25/2005
This part of code gets the bounding box of the local mesh.
***/
xmin=E->Tracer.x_space[1];
xmax=E->Tracer.x_space[nox];
ymin=E->Tracer.y_space[1];
ymax=E->Tracer.y_space[noy];
zmin=E->Tracer.z_space[1];
zmax=E->Tracer.z_space[noz];
/*fprintf(stderr,"%d %d\n", E->parallel.nprocx, E->parallel.loc2proc_map[0][0][1][0]);*/
/*** comment by Tan2 1/25/2005
Copy the velocity array.
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(i=1;i<=nno;i++)
for(j=1;j<=3;j++) {
E->GV[m][j][i]=E->sphere.cap[m].V[j][i];
}
}
/*** comment by vlad 03/17/2005
advecting tracers in each processor
***/
for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
n_x=0;
n_y=0;
n_z=0;
/*printf("%s %d %s %d %s %f\n","La pasul",E->monitor.solution_cycles, "Procesorul", E->parallel.me, "advecteaza tracerul",E->Tracer.tracer_y[n]);
//printf("%d %d %d %f %f %f %f\n", E->monitor.solution_cycles, E->parallel.me,E->Tracer.LOCAL_NUM_TRACERS,E->Tracer.itcolor[n],E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n]);
//printf("%d %d %d %f %f %f %f\n", E->monitor.solution_cycles, E->parallel.me,E->Tracer.LOCAL_NUM_TRACERS,tr_color[n],tr_x[n],tr_y[n],tr_z[n]);
*/
/* mid point method uses 2 iterations */
for(iteration=1;iteration<=2;iteration++)
{
if(iteration ==1) {
x_tmp=E->Tracer.tracer_x[n];
y_tmp=E->Tracer.tracer_y[n];
z_tmp=E->Tracer.tracer_z[n];
}
/*** comment by Tan2 1/25/2005
Find the element that contains the tracer.
nodex n_x n_x+1
| * |
<----------->
tr_dx
<-------------------------------->
dx
***/
for(i=1;i<nox;i++) {
if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
tr_dx=x_tmp-E->Tracer.x_space[i];
dx=E->Tracer.x_space[i+1]-E->Tracer.x_space[i];
n_x=i;
}
}
for(j=1;j<noy;j++) {
if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
tr_dy=y_tmp-E->Tracer.y_space[j];
dy=E->Tracer.y_space[j+1]-E->Tracer.y_space[j];
n_y=j;
}
}
for(k=1;k<noz;k++) {
if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
tr_dz=z_tmp-E->Tracer.z_space[k];
dz=E->Tracer.z_space[k+1]-E->Tracer.z_space[k];
n_z=k;
}
}
/*** comment by Tan2 1/25/2005
Calculate shape functions from tr_dx, tr_dy, tr_dz
This assumes linear element
***/
/* compute volumetic weighting functions */
w1=tr_dx*tr_dz*tr_dy;
w2=(dx-tr_dx)*tr_dz*tr_dy;
w3=tr_dx*(dz-tr_dz)*tr_dy;
w4=(dx-tr_dx)*(dz-tr_dz)*tr_dy;
w5=tr_dx*tr_dz*(dy-tr_dy);
w6=(dx-tr_dx)*tr_dz*(dy-tr_dy);
w7=tr_dx*(dz-tr_dz)*(dy-tr_dy);
w8=(dx-tr_dx)*(dz-tr_dz)*(dy-tr_dy);
volume=dx*dz*dy;
/*** comment by Tan2 1/25/2005
Calculate the 8 node numbers of current element
***/
node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
node6 = n_z + n_x*noz + n_y*noz*nox;
node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
node8 = n_z+1 + n_x*noz + n_y*noz*nox;
/* printf("%d %d %d %d %d %d %d %d %d\n",E->parallel.me, node1,node2,node3,node4,node5,node6,node7,node8);
//printf("%d %f %f %f %f %f %f %f %f\n", E->parallel.me, E->GV[1][2][node1], E->GV[1][2][node2], E->GV[1][2][node3], E->GV[1][2][node4], E->GV[1][2][node5], E->GV[1][2][node6], E->GV[1][2][node7], E->GV[1][2][node8]);
*/
/*** comment by Tan2 1/25/2005
Interpolate the velocity on the tracer position
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(j=1;j<=3;j++) {
tr_v[m][j]=w8*E->GV[m][j][node1]
+w7*E->GV[m][j][node2]
+w6*E->GV[m][j][node3]
+w5*E->GV[m][j][node4]
+w4*E->GV[m][j][node5]
+w3*E->GV[m][j][node6]
+w2*E->GV[m][j][node7]
+w1*E->GV[m][j][node8];
tr_v[m][j]=tr_v[m][j]/volume;
}
/*** comment by Tan2 1/25/2005
advect tracer using mid-point method (2nd order accuracy)
***/
/* mid point method */
if(iteration == 1) {
x_tmp = x_tmp + (E->advection.timestep/2.0)*tr_v[m][1]/E->Tracer.z_space[n_z];
y_tmp = y_tmp + (E->advection.timestep/2.0)*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
z_tmp = z_tmp + (E->advection.timestep/2.0)*tr_v[m][3];
}
if( iteration == 2) {
E->Tracer.tracer_x[n] += E->advection.timestep*tr_v[m][1]/E->Tracer.z_space[n_z];
E->Tracer.tracer_y[n] += E->advection.timestep*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
E->Tracer.tracer_z[n] += E->advection.timestep*tr_v[m][3];
}
} /* end of m */
} /* end of iteration loop */
}/* end of tracer loop */
/*** Comment by Vlad 3/25/2005
MPI for the tracer-advection code
***/
m = 0;
locx = E->parallel.me_loc[1];
locy = E->parallel.me_loc[2];
locz = E->parallel.me_loc[3];
/* Am I the left-most proc.? If not, who is on my left? */
if (locy == 0)
left = -1;
else
left = E->parallel.loc2proc_map[m][locx][locy-1][locz];
/* Am I the right-most proc.? If not, who is on my right? */
if (locy == E->parallel.nprocy-1)
right = -1;
else
right = E->parallel.loc2proc_map[m][locx][locy+1][locz];
/* Am I the lower-most proc.? If not, who is beneath me? */
if (locz == 0)
down = -1;
else
down = E->parallel.loc2proc_map[m][locx][locy][locz-1];
/* Am I the upper-most proc.? If not, who is above me? */
if (locz == E->parallel.nprocz-1)
up = -1;
else
up = E->parallel.loc2proc_map[m][locx][locy][locz+1];
/* Am I the back-most proc.? If not, who is behind me? */
if (locx == 0)
back = -1;
else
back = E->parallel.loc2proc_map[m][locx-1][locy][locz];
/* Am I the front-most proc.? If not, who is in front of me? */
if (locx == E->parallel.nprocx-1)
front = -1;
else
front = E->parallel.loc2proc_map[m][locx+1][locy][locz];
Left_proc_list[0]=left;
Left_proc_list[1]=right;
Left_proc_list[2]=down;
Left_proc_list[3]=up;
Left_proc_list[4]=back;
Left_proc_list[5]=front;
Right_proc_list[0]=right;
Right_proc_list[1]=left;
Right_proc_list[2]=up;
Right_proc_list[3]=down;
Right_proc_list[4]=front;
Right_proc_list[5]=back;
jump_new[0]=0;
jump_new[1]=0;
jump_new[2]=0;
jump_new[3]=0;
jump_new[4]=0;
jump_new[5]=0;
count_new[0]=0;
count_new[1]=0;
count_new[2]=0;
count_new[3]=0;
count_new[4]=0;
count_new[5]=0;
jj[0]=1;
jj[1]=0;
jj[2]=3;
jj[3]=2;
jj[4]=5;
jj[5]=4;
E->Tracer.TEMP_TRACERS=0;
Current_num_tracers=0;
for(i=0;i<=11;i++){
for(j=0;j<=E->Tracer.NUM_TRACERS;j++){
tr_color_new[i][j]=999;
tr_x_new[i][j]=999;
tr_y_new[i][j]=999;
tr_z_new[i][j]=999;
tr_color_1[j]=999;
tr_x_1[j]=999;
tr_y_1[j]=999;
tr_z_1[j]=999;
}
}
i=0;
j=0;
k=0;
l=0;
m=0;
o=0;
p=0;
for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++){
if(E->Tracer.tracer_y[n]>ymax) {
if(E->Tracer.tracer_y[n]+100 != 100) {
tr_color_new[0][i]=E->Tracer.itcolor[n];
tr_x_new[0][i]=E->Tracer.tracer_x[n];
tr_y_new[0][i]=E->Tracer.tracer_y[n];
tr_z_new[0][i]=E->Tracer.tracer_z[n];
i++;
jump_new[0]=i;
}
}
else if(E->Tracer.tracer_y[n]<ymin) {
if(E->Tracer.tracer_y[n]+100 != 100) {
tr_color_new[1][j]=E->Tracer.itcolor[n];
tr_x_new[1][j]=E->Tracer.tracer_x[n];
tr_y_new[1][j]=E->Tracer.tracer_y[n];
tr_z_new[1][j]=E->Tracer.tracer_z[n];
j++;
jump_new[1]=j;
}
}
else if(E->Tracer.tracer_z[n]>zmax) {
if(E->Tracer.tracer_z[n]+100 != 100) {
tr_color_new[2][k]=E->Tracer.itcolor[n];
tr_x_new[2][k]=E->Tracer.tracer_x[n];
tr_y_new[2][k]=E->Tracer.tracer_y[n];
tr_z_new[2][k]=E->Tracer.tracer_z[n];
k++;
jump_new[2]=k;
}
}
else if(E->Tracer.tracer_z[n]<zmin) {
if(E->Tracer.tracer_z[n]+100 != 100) {
tr_color_new[3][l]=E->Tracer.itcolor[n];
tr_x_new[3][l]=E->Tracer.tracer_x[n];
tr_y_new[3][l]=E->Tracer.tracer_y[n];
tr_z_new[3][l]=E->Tracer.tracer_z[n];
l++;
jump_new[3]=l;
}
}
else if(E->Tracer.tracer_x[n]>xmax) {
if(E->Tracer.tracer_x[n]+100 != 100) {
tr_color_new[4][m]=E->Tracer.itcolor[n];
tr_x_new[4][m]=E->Tracer.tracer_x[n];
tr_y_new[4][m]=E->Tracer.tracer_y[n];
tr_z_new[4][m]=E->Tracer.tracer_z[n];
m++;
jump_new[4]=m;
}
}
else if(E->Tracer.tracer_x[n]<xmin) {
if(E->Tracer.tracer_x[n]+100 != 100) {
tr_color_new[5][o]=E->Tracer.itcolor[n];
tr_x_new[5][o]=E->Tracer.tracer_x[n];
tr_y_new[5][o]=E->Tracer.tracer_y[n];
tr_z_new[5][o]=E->Tracer.tracer_z[n];
o++;
jump_new[5]=o;
}
}
else {
tr_color_1[p]=E->Tracer.itcolor[n];
tr_x_1[p]=E->Tracer.tracer_x[n];
tr_y_1[p]=E->Tracer.tracer_y[n];
tr_z_1[p]=E->Tracer.tracer_z[n];
p++;
}
}
Previous_num_tracers=E->Tracer.LOCAL_NUM_TRACERS;
Current_num_tracers=Previous_num_tracers-jump_new[0]-jump_new[1]-jump_new[2]-jump_new[3]-jump_new[4]-jump_new[5];
for(p=1;p<=Current_num_tracers;p++){
E->Tracer.itcolor[p]=tr_color_1[p-1];
E->Tracer.tracer_x[p]=tr_x_1[p-1];
E->Tracer.tracer_y[p]=tr_y_1[p-1];
E->Tracer.tracer_z[p]=tr_z_1[p-1];
}
for(i=0;i<=5;i++){
j=jj[i];
if (Left_proc_list[i] >= 0) {
proc=Left_proc_list[i];
MPI_Irecv(tr_color_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 11+i, world, &request[0]);
MPI_Irecv(tr_x_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 12+i, world, &request[1]);
MPI_Irecv(tr_y_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 13+i, world, &request[2]);
MPI_Irecv(tr_z_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 14+i, world, &request[3]);
}
if (Right_proc_list[i] >= 0) {
proc=Right_proc_list[i];
MPI_Send(tr_color_new[i], jump_new[i], MPI_FLOAT, proc, 11+i, world);
MPI_Send(tr_x_new[i], jump_new[i], MPI_FLOAT, proc, 12+i, world);
MPI_Send(tr_y_new[i], jump_new[i], MPI_FLOAT, proc, 13+i, world);
MPI_Send(tr_z_new[i], jump_new[i], MPI_FLOAT, proc, 14+i, world);
}
if (Left_proc_list[i] >= 0) {
MPI_Waitall(4, request, status);
status_count = status[0];
MPI_Get_count(&status_count, MPI_FLOAT, &count_new[i]);
}
E->Tracer.TEMP_TRACERS=E->Tracer.TEMP_TRACERS+count_new[i]-jump_new[i];
E->Tracer.LOCAL_NUM_TRACERS=Previous_num_tracers+E->Tracer.TEMP_TRACERS;
if(i <= 1){
for(n=Current_num_tracers+count_new[j];n<=Current_num_tracers+count_new[i]+count_new[j]-1;n++) {
m=Current_num_tracers+count_new[j];
E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
}
}
else if (i <= 4) {
for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];n<=Current_num_tracers+count_new[0]+count_new[1]+count_new[i]+count_new[j]-1;n++) {
m=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];
E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
}
}
else {
for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];n<=E->Tracer.LOCAL_NUM_TRACERS-1;n++) {
m=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];
E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
}
}
}
free (tr_color_1);
free (tr_x_1);
free (tr_y_1);
free (tr_z_1);
for(i=0;i<=11;i++){
free (tr_color_new[i]);
free (tr_x_new[i]);
free (tr_y_new[i]);
free (tr_z_new[i]);
}
return;
}
Computing file changes ...