/* This file is to convert binary .prn file from Taras' I3ELVIS code into VTK-files in ASCII format */
/* and interpolate marker's compostion field, temperature and viscosity into rectilinear Grid */
/* with specified resolution(xresol,yresol,zresol) */
/**/
/* Original: Guizhi Zhu, September 8,2010 */
/* Modified and integrated into i3mg.c by Ria Fischer, April 2012 */
/* Added interpolation from marker's c1,c2 field, Dec. 2013 */

void prn2vtr(int f0)
{
	/* Load Information from data file ------------------------------- */
	/* void loader() */
	/* bondv[] - bondary value */
	/* bondm[] - bondary mode 0=Not, -1=Value, 1,2...=LinNum+1 */
	/* m1,m2,m3 - node X,Y,Z number */
	/* Counter */
	int n1;
	char nn1, sf0[15], vtk_line[200];
	char *pos;
	long int m1, m2, m3, m4;
	long int mm1, mcmax;
	char szint, szlong, szfloat, szdouble, szcur;
	char intitule[] = "# vtk DataFile Version 2.0";
	char data_title[] = "\ntest ecrit. coord \n";
	char type_ascii[] = "ASCII \n";
	char dataset_type[] = "DATASET RECTILINEAR_GRID \n";
	char dimensions[] = "DIMENSIONS ";
	float ival0;
	double ival1;
	char file_title[] = "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
	/* interpolate compositon field with specified resolution, you can change it to meet your need */
	/* XXX switch for higher resolution */
	xresol = (xnumx - 1) * 1 + 1;
	yresol = (ynumy - 1) * 1 + 1;
	zresol = (znumz - 1) * 1 + 1;
	nresol = xresol * yresol * zresol;
	/**/
	if (maxmat < nresol)
	{
		printf("\nError: Space out in val0, val1 or lin0!\n");
		exit(0);
	}
	/* Define input and output file name */
	strcpy(fl1in, fl1out);
	pos = strstr(fl1out, "_");
	//if(pos<1)
	//{
	//	//ERROR: wrong format of name: Please choose name: name_nr.prn
	// 	printf("Wrong format of file name: Please choose: name_nr.prn"; exit(0);}
	//}
	strncpy(pos, "_", 2);
	//*(pos+1) = (int)0;
	sprintf(sf0, "%i.vtr", f0);
	strcat(fl1out, sf0);
	
	/*write composition,temperature and viscosity in binary format according to what resolution you specified */
	fl2 = fopen(fl1out, "w");
	if (fl2 == 0)
	{
		printf("\nError: Cannot open file %s!\n", fl1out);
		exit(0);
	}
	
	/* for markers' composition , please note the resolution*/
	/* Check coordinates, calc step */
	xminx = 0;
	xmaxx = xsize;
	xresx = (xmaxx - xminx) / (double)(xresol - 1);
	yminy = 0;
	ymaxy = ysize;
	yresy = (ymaxy - yminy) / (double)(yresol - 1);
	zminz = 0;
	zmaxz = zsize;
	zresz = (zmaxz - zminz) / (double)(zresol - 1);
	// Interpolate markers onto nodes
	compo_water_recalc();
	
	/**** it is the following part to write xml file ****/
	szint = sizeof(int);
	szlong = sizeof(long);
	szfloat = sizeof(float);
	szdouble = sizeof(double);
	
	//printf("szint is  %d; szfloat is %d",szint,szfloat );
	
	fprintf(fl2, "%s", file_title);
	fprintf(fl2, "\t<RectilinearGrid WholeExtent=\" %d %ld %d %ld %d %ld \">\n", 0, xresol - 1, 0, yresol - 1, 0,
			zresol - 1);
	fprintf(fl2, "\t\t<Piece Extent=\"%d %ld %d %ld %d %ld \">\n", 0, xresol - 1, 0, yresol - 1, 0, zresol - 1);
	fprintf(fl2, "\t\t\t<PointData Scalars=\"Composition\" Vectors=\"Velocity\">\n");
	
	m4 = 0;
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Granite\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Dacite\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Temperature\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Viscosity\"   format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Density\"     format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Pressure\"     format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Velocity\" NumberOfComponents=\"3\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld\"/>\n",
			m4);
	m4 += nresol * szfloat * 3 + szint;
	
	fprintf(fl2, "\t\t\t</PointData>\n");
	fprintf(fl2, "\t\t\t<CellData>\n");
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Int32\"   Name=\"Composition\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\" %ld \"/>\n",
			m4);
	m4 += (xresol - 1) * (yresol - 1) * (zresol - 1) * szint + szint;
	
	fprintf(fl2, "\t\t\t</CellData>\n");
	fprintf(fl2, "\t\t\t<Coordinates>\n");
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Xcoordinate\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\"%ld\"/>\n",
			m4);
	m4 += xresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Ycoordinate\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\"%ld\"/>\n",
			m4);
	m4 += yresol * szfloat + szint;
	
	fprintf(fl2,
			"\t\t\t\t<DataArray type=\"Float32\" Name=\"Zcoordinate\" format=\"appended\" RangeMin=\"\" RangeMax=\"\" offset=\"%ld\"/>\n",
			m4);
	m4 += zresol * szfloat + szint;
	
	fprintf(fl2, "\t\t\t</Coordinates>\n");
	fprintf(fl2, "\t\t</Piece>\n");
	fprintf(fl2, "\t</RectilinearGrid>\n");
	fprintf(fl2, "\t<AppendedData encoding=\"raw\">\n");
	fprintf(fl2, "_");
	
	/* firstly write the length of the following granite data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following granite data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(c1[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}

	/* firstly write the length of the following dacite data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following dacite data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(c2[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}

	/* firstly write the length of the following temperature data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following temperature data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(tk_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}

	/* firstly write the length of the following viscosity data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following viscosity data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(nu_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}
	
	/* firstly write the length of the following density data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following density data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(ro_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}
	
	/* firstly write the length of the following pressure data to read */
	m1 = nresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following pressure data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(pr_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}
	
	/* firstly write the length of the following velocity vector data to read */
	m1 = nresol * szfloat * 3;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following velocity vector data to read */
	for (m3 = 0; m3 < zresol; m3++)
		for (m2 = 0; m2 < yresol; m2++)
			for (m1 = 0; m1 < xresol; m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				ival0 = (float)(vx_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
				ival0 = (float)(vy_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
				ival0 = (float)(vz_vtk[mcmax]);
				fwrite(&ival0, szfloat, 1, fl2);
			}
	
	/* firstly write the length of the following composition data to read */
	m1 = (xresol - 1) * (yresol - 1) * (zresol - 1) * szint;
	fwrite(&m1, szint, 1, fl2);
	/* then write the following composition data to read , seperate new crust and others*/
	for (m3 = 0; m3 < (zresol - 1); m3++)
		for (m2 = 0; m2 < (yresol - 1); m2++)
			for (m1 = 0; m1 < (xresol - 1); m1++)
			{
				mcmax = m3 * (xresol * yresol) + yresol * m1 + m2;
				/* if (lin0[mcmax]==16)
				 {n1=2;}
				 else
				 {n1=1;}
				 */
				n1 = lin0[mcmax];
				fwrite(&n1, szint, 1, fl2);

			}
	
	/* firstly write the length of the xcoordinates to read */
	m1 = xresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the x-coordinates */
	for (m1 = 0; m1 < xresol; m1++)
	{
		ival0 = (float)(m1 * xresx / 1000.);
		fwrite(&ival0, szfloat, 1, fl2);
	}
	
	/* firstly write the length of the ycoordinates to read */
	m1 = yresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the y-coordinates */
	for (m2 = 0; m2 < yresol; m2++)
	{
		ival0 = (float)(m2 * yresy / 1000.);
		fwrite(&ival0, szfloat, 1, fl2);
	}
	
	/* firstly write the length of the zcoordinates to read */
	m1 = zresol * szfloat;
	fwrite(&m1, szint, 1, fl2);
	/* then write the z-coordinates */
	for (m3 = 0; m3 < zresol; m3++)
	{
		ival0 = (float)(m3 * zresz / 1000.);
		fwrite(&ival0, szfloat, 1, fl2);
	}
	
	fprintf(fl2, "\n\t</AppendedData>\n");
	fprintf(fl2, "</VTKFile>\n");
	/**** it is over now, to celebrate finishing it ****/

	fclose(fl2);
	
	/* write timestep and file name to time.pvd file */
	fl = fopen("time.pvd", "r+");
	if (fl == 0)
	{
		printf("\nError: Cannot open file 'time.pvd'!\n");
		exit(0);
	}
	fseek(fl, -25, SEEK_END);
	sprintf(vtk_line, "\t\t<DataSet timestep=\"%1.12lf \" file=\"./%s\"/>\n\t</Collection>\n</VTKFile>",
			timesum / 3.15576e+13, fl1out);
	fputs(vtk_line, fl);
	fclose(fl);
}

/* Divide mesh into Voronoi cells */
/* by assigning a closest marker to each nodal point of the mesh */
/* which allows to assign marker properties to the nodal points of the mesh */
void compo_water_recalc()
{
	/* Counters */
	long int m1, m2, m3, m4, m5, m6, m7, m8, m9;
	long int m80;
	long int m10, m20, m30, m10c, m20c, m30c, m10v, m20v, m30v;
	long int m40, m40pr, m40vx, m40vy, m40vz;
	long int mm1, mm10;

	double x, y, z;
	double dx, dy, dz, dxc, dyc, dzc, dxv, dyv, dzv;
	double rx, ry, rz, rxc, ryc, rzc, rxv, ryv, rzv;
	double ival, ival0, ival1, ival2, ival3, ival4, ival5, ival6, ival7;
	double ea, na, ua, ea1, na1, ua1;
	double dea, dna, dua, dea1, dna1, dua1;
	double swt0, swt1, swt2, swt3, swt4, swt5, swt6, swt7;

	/* compute composition */
	/* initial value */
	for (m1 = 0; m1 < nresol; m1++)
	{
		val0[m1] = 1e+30;		// distance of marker
		lin0[m1] = -1;			// number of marker
		sol2[m1] = 0.0;
		c1[m1] = 0.0;			// granite composition
		c2[m1] = 0.0;			// dacite composition
	}
	/* search for markers */
	/* Assign each marker to the closest node using L2-norm -> Voronoi tesselation */
	for (m4 = 0; m4 < marknum; m4++)
		if ((markt[m4] < 50 || markt[m4] >= 100) && (double)(markx[m4]) > xminx && (double)(markx[m4]) < xmaxx
				&& (double)(marky[m4]) > yminy && (double)(marky[m4]) < ymaxy && (double)(markz[m4]) > zminz
				&& (double)(markz[m4]) < zmaxz)
		{
			if (markt[m4] >= 100) markt[m4] -= 100;
			
			/* Define relative position of the marker */
			ea = ((double)(markx[m4])) / xresx;
			m1 = (long int)(ea);
			if (m1 < 0) m1 = 0;
			if (m1 > xresol - 2) m1 = xresol - 2;
			na = ((double)(marky[m4])) / yresy;
			m2 = (long int)(na);
			if (m2 < 0) m2 = 0;
			if (m2 > yresol - 2) m2 = yresol - 2;
			ua = ((double)(markz[m4])) / zresz;
			m3 = (long int)(ua);
			if (m3 < 0) m3 = 0;
			if (m3 > zresol - 2) m3 = zresol - 2;

			/* Calc distances and weights to surrounding nodes */
			ea = markx[m4] - ((double)(m1)) * xresx;
			ea1 = ABSV(ea - xresx);
			na = marky[m4] - ((double)(m2)) * yresy;
			na1 = ABSV(na - yresy);
			ua = markz[m4] - ((double)(m3)) * zresz;
			ua1 = ABSV(ua - zresz);

			dea = (1 - ea / xresx);
			dea1 = (1 - dea);
			dna = (1 - na / yresy);
			dna1 = (1 - dna);
			dua = (1 - ua / zresz);
			dua1 = (1 - dua);
			swt0 = dua * dea * dna;
			swt1 = dua * dea * dna1;
			swt2 = dua * dea1 * dna;
			swt3 = dua * dea1 * dna1;
			swt4 = dua1 * dea * dna;
			swt5 = dua1 * dea * dna1;
			swt6 = dua1 * dea1 * dna;
			swt7 = dua1 * dea1 * dna1;

			ea *= ea;
			ea1 *= ea1;
			na *= na;
			na1 *= na1;
			ua *= ua;
			ua1 *= ua1;

			/* Check distances to each of the eight nodes surrounding the marker */
			/* Find node which is closest to marker */
			//
			//     4---6
			//    /|  /|
			//   / 5-/-7
			//  / / / /
			// 0---2 /
			// |/  |/
			// 1---3
			//
			//	nodes are addressed in the following order: 0,1,3,2,6,7,5,4
			//
			m5 = m3 * xresol * yresol + m1 * yresol + m2;
			ival0 = sqrt(ua + ea + na);
			ival1 = sqrt(ua + ea + na1);
			ival2 = sqrt(ua + ea1 + na);
			ival3 = sqrt(ua + ea1 + na1);
			ival4 = sqrt(ua1 + ea + na);
			ival5 = sqrt(ua1 + ea + na1);
			ival6 = sqrt(ua1 + ea1 + na);
			ival7 = sqrt(ua1 + ea1 + na1);
			{
				m6 = m5;
				if (swt0 > 0)
				{
					ival = markc1[m4] * swt0;
					c1[m6] += ival;
					ival = markc2[m4] * swt0;
					c2[m6] += ival;
					sol2[m6] += swt0;
				}
				ival = val0[m6];
				if (ival > ival0)
				{
					val0[m6] = ival0;
					lin0[m6] = m4;
				}
				m6 += 1;
				if (swt1 > 0)
				{
					ival = markc1[m4] * swt1;
					c1[m6] += ival;
					ival = markc2[m4] * swt1;
					c2[m6] += ival;
					sol2[m6] += swt1;
				}
				ival = val0[m6];
				if (ival > ival1)
				{
					val0[m6] = ival1;
					lin0[m6] = m4;
				}
				m6 += yresol;
				if (swt3 > 0)
				{
					ival = markc1[m4] * swt3;
					c1[m6] += ival;
					ival = markc2[m4] * swt3;
					c2[m6] += ival;
					sol2[m6] += swt3;
				}
				ival = val0[m6];
				if (ival > ival3)
				{
					val0[m6] = ival3;
					lin0[m6] = m4;
				}
				m6 -= 1;
				if (swt2 > 0)
				{
					ival = markc1[m4] * swt2;
					c1[m6] += ival;
					ival = markc2[m4] * swt2;
					c2[m6] += ival;
					sol2[m6] += swt2;
				}
				ival = val0[m6];
				if (ival > ival2)
				{
					val0[m6] = ival2;
					lin0[m6] = m4;
				}
				m6 += (yresol * xresol);
				if (swt6 > 0)
				{
					ival = markc1[m4] * swt6;
					c1[m6] += ival;
					ival = markc2[m4] * swt6;
					c2[m6] += ival;
					sol2[m6] += swt6;
				}
				ival = val0[m6];
				if (ival > ival6)
				{
					val0[m6] = ival6;
					lin0[m6] = m4;
				}
				m6 += 1;
				if (swt7 > 0)
				{
					ival = markc1[m4] * swt7;
					c1[m6] += ival;
					ival = markc2[m4] * swt7;
					c2[m6] += ival;
					sol2[m6] += swt7;
				}
				ival = val0[m6];
				if (ival > ival7)
				{
					val0[m6] = ival7;
					lin0[m6] = m4;
				}
				m6 -= yresol;
				if (swt5 > 0)
				{
					ival = markc1[m4] * swt5;
					c1[m6] += ival;
					ival = markc2[m4] * swt5;
					c2[m6] += ival;
					sol2[m6] += swt5;
				}
				ival = val0[m6];
				if (ival > ival5)
				{
					val0[m6] = ival5;
					lin0[m6] = m4;
				}
				m6 -= 1;
				if (swt4 > 0)
				{
					ival = markc1[m4] * swt4;
					c1[m6] += ival;
					ival = markc2[m4] * swt4;
					c2[m6] += ival;
					sol2[m6] += swt4;
				}
				ival = val0[m6];
				if (ival > ival4)
				{
					val0[m6] = ival4;
					lin0[m6] = m4;
				}
			}
		}/*end if, end for loop*/
	/* Search for nodes which have not been assigned to a marker and search for the closest marker */
	for (m3 = 0; m3 < zresol; m3++)
		for (m1 = 0; m1 < xresol; m1++)
			for (m2 = 0; m2 < yresol; m2++)
			{
				m4 = m3 * xresol * yresol + m1 * yresol + m2;
				if (lin0[m4] == -1)
				{
					val0[m4] = 1e30;
					m5 = 1;
					
					ival1 = ((double)(m1)) * xresx;
					ival2 = ((double)(m2)) * yresy;
					ival3 = ((double)(m3)) * zresz;
					
					do
					{
						for (m6 = m3 - m5; m6 <= m3 + m5; m6++)
							for (m7 = m1 - m5; m7 <= m1 + m5; m7++)
								for (m8 = m2 - m5; m8 <= m2 + m5; m8++)
									if ((m7 == m1 - m5 || m7 == m1 + m5 || m8 == m2 - m5 || m8 == m2 + m5
											|| m6 == m3 - m5 || m6 == m3 + m5) && m7 >= 0 && m7 < xresol && m8 >= 0
											&& m8 < yresol && m6 >= 0 && m6 < zresol)
									{
										m9 = m6 * yresol * xresol + m7 * yresol + m8;
										
										m80 = lin0[m9];
										ival = val0[m9];
										if (m80 >= 0 && ival > 0)
										{
											ea = markx[m80] - ival1;
											ea *= ea;
											na = marky[m80] - ival2;
											na *= na;
											ua = markz[m80] - ival3;
											ua *= ua;
											ival0 = sqrt(ea + na + ua);
											if (val0[m4] > ival0)
											{
												val0[m4] = ival0;
												lin0[m4] = -(m80 + 2);
											}
										}
									}
						m5++;
					}
					while (m5 <= 3);
				}
				if (sol2[m4] == 0)
				{
					m5 = 1;
					swt0 = 0;

					do
					{
						for (m6 = m3 - m5; m6 <= m3 + m5; m6++)
							for (m7 = m1 - m5; m7 <= m1 + m5; m7++)
								for (m8 = m2 - m5; m8 <= m2 + m5; m8++)
									if ((m7 == m1 - m5 || m7 == m1 + m5 || m8 == m2 - m5 || m8 == m2 + m5
											|| m6 == m3 - m5 || m6 == m3 + m5) && m7 >= 0 && m7 < xresol && m8 >= 0
											&& m8 < yresol && m6 >= 0 && m6 < zresol)
									{
										m9 = m6 * yresol * xresol + m7 * yresol + m8;

										if (sol2[m9] >= 0)
										{
											// weights in each direction
											dea = (2 * m5 - (((double)(m7)) - ((double)(m1)))) / (2 * m5);
											dna = (2 * m5 - (((double)(m8)) - ((double)(m2)))) / (2 * m5);
											dua = (2 * m5 - (((double)(m6)) - ((double)(m3)))) / (2 * m5);
											// total weight
											swt0 = dua * dea * dna;

											ival = c1[m9] / sol2[m9] * swt0;
											c1[m4] += ival;
											ival = c2[m9] / sol2[m9] * swt0;
											c2[m4] += ival;
											sol2[m4] += swt0;
										}
									}
						m5++;
						if (m5 > 10) printf("\nError in prn2vtr: no marker found for node %ld!\n", m4);
					}
					while (sol2[m4] == 0);
				}
			}
	
	/*****/
	/* Assign compositional type */
	/* Assign values of c1 and c2 */
	for (m1 = 0; m1 < nresol; m1++)
	{
		if (lin0[m1] >= 0)
		{
			m2 = markt[lin0[m1]];
			lin0[m1] = m2;
		}
		if (lin0[m1] < -1)
		{
			m2 = markt[-lin0[m1] - 2];
			lin0[m1] = m2;
		}
		if (sol2[m1] > 0)
		{
			c1[m1] /= sol2[m1];
			c2[m1] /= sol2[m1];
		}
		else
		{
			printf("\nError in prn2vtr: Cannot assign composition. No marker found for node %ld!\n", m1);
			exit(0);
		}
	}
	/**/
	/* End CHEMICAL COMPONENT INTERPOLATION ----------------------------*/

	/* To start melt extraction */
	/* Clear wt */
	/*for (m4=0;m4<nresol;m4++)
	 {
	 
	 extract_cell0[m4]=0;
	 
	 }
	 */
	/*gridmod=10;
	 
	 
	 for (m4=0;m4<marknum;m4+=gridmod)
	 {
	 if ((markt[m4]<50 ) && (double)(markx[m4])>xminx && (double)(markx[m4])<xmaxx && (double)(marky[m4])>yminy && (double)(marky[m4])<ymaxy&& (double)(markz[m4])>zminz && (double)(markz[m4])<zmaxz)
	 {
	 ea=((double)(markx[m4]))/xresx;
	 m10=(long int)(ea);
	 if(m10<0) m10=0;
	 if(m10>xresol-2) m10=xresol-2;
	 na=((double)(marky[m4]))/yresy;
	 m20=(long int)(na);
	 if(m20<0) m20=0;
	 if(m20>yresol-2) m20=yresol-2;
	 ua=((double)(markz[m4]))/zresz;
	 m30=(long int)(ua);
	 if(m30<0) m30=0;
	 if(m30>zresol-2) m30=zresol-2;
	 for(m3=m30;m3<=m30+1;m3++)
	 {
	 ddz=(markz[m4]-((double)(m30))*zresz)/zresz;
	 ddz=1.0-ABSV(ddz);
	 for(m1=m10;m1<=m10+1;m1++)
	 {
	 ddx=(markx[m4]-((double)(m10))*xresx)/xresx;
	 ddx=1.0-ABSV(ddx);
	 for(m2=m20;m2<=m20+1;m2++)
	 {
	 ddy=(marky[m4]-((double)(m20))*yresy)/yresy;
	 ddy=1.0-ABSV(ddy);
	 swt=ddx*ddy*ddz;
	 m5=m3*xresol*yresol+m1*yresol+m2;
	 sol0[m5]+=swt;
	 extract_cell0[m5]+=markex[m4]*swt;
	 }
	 
	 }
	 
	 }
	 }
	 
	 
	 }
	 for (m1=0;m1<nresol;m1++)
	 {
	 if(sol0[m1])
	 {
	 extract_cell[m1]=extract_cell0[m1]/sol0[m1];
	 sol0[m1]=0;
	 
	 }
	 }
	 
	 */
	/**/
	/* Visualisation parameters output preparation */
	for (m3 = 0; m3 < zresol; m3++)
		for (m1 = 0; m1 < xresol; m1++)
			for (m2 = 0; m2 < yresol; m2++)
			{
				/* Define visualization node */
				m4 = m3 * xresol * yresol + m1 * yresol + m2;
				/* Coordinates */
				x = (double)(m1) * xresx;
				y = (double)(m2) * yresy;
				z = (double)(m3) * zresz;
				/* Out of grid */
				if (x < 0) x = 0;
				if (x > xsize) x = xsize;
				if (y < 0) y = 0;
				if (y > ysize) y = ysize;
				if (z < 0) z = 0;
				if (z > zsize) z = zsize;
				/* Find reference node: Left, upper, frontal corner of cell */
				m10 = m1serch(x);
				m20 = m2serch(y);
				m30 = m3serch(z);
				/* Define P reference node at cell center */
				m10c = m10;
				m20c = m20;
				m30c = m30;
				/* Define V reference node */
				m10v = m10;
				m20v = m20;
				m30v = m30;
				/* Shift pressure nodes by half a cell towards origin */
				if (x > (gx[m10c] + gx[m10c + 1]) / 2.0) m10c += 1;
				if (y > (gy[m20c] + gy[m20c + 1]) / 2.0) m20c += 1;
				if (z > (gz[m30c] + gz[m30c + 1]) / 2.0) m30c += 1;
				/* Special case: left, frontal, upper boundary conditions */
				if (m10c < 1) m10c = 1;
				if (m20c < 1) m20c = 1;
				if (m30c < 1) m30c = 1;
				/* Special case: right, back, lower boundary conditions */
				if (m10c > xnumx - 2) m10c = xnumx - 2;
				if (m20c > ynumy - 2) m20c = ynumy - 2;
				if (m30c > znumz - 2) m30c = znumz - 2;
				/* Shift velocity nodes by half a cell away from origin */
				if (x < (gx[m10v] + gx[m10v + 1]) / 2.0) m10v -= 1;
				if (y < (gy[m20v] + gy[m20v + 1]) / 2.0) m20v -= 1;
				if (z < (gz[m30v] + gz[m30v + 1]) / 2.0) m30v -= 1;
				/* Special case: left, frontal, upper boundary conditions */
				if (m10v < 0) m10v = 0;
				if (m20v < 0) m20v = 0;
				if (m30v < 0) m30v = 0;
				/* Special case: right, back, lower boundary condition */
				if (m10v > xnumx - 3) m10v = xnumx - 3;
				if (m20v > ynumy - 3) m20v = ynumy - 3;
				if (m30v > znumz - 3) m30v = znumz - 3;
				/* Cell dimensions */
				dx = gx[m10 + 1] - gx[m10];
				dy = gy[m20 + 1] - gy[m20];
				dz = gz[m30 + 1] - gz[m30];
				/* P cell dimensions */
				dxc = (gx[m10c + 1] - gx[m10c - 1]) / 2.0;
				dyc = (gy[m20c + 1] - gy[m20c - 1]) / 2.0;
				dzc = (gz[m30c + 1] - gz[m30c - 1]) / 2.0;
				/* V cell dimensions */
				dxv = (gx[m10v + 2] - gx[m10v]) / 2.0;
				dyv = (gy[m20v + 2] - gy[m20v]) / 2.0;
				dzv = (gz[m30v + 2] - gz[m30v]) / 2.0;
				/* Normalized distance of visualization node to reference node */
				rx = (x - gx[m10]) / dx;
				ry = (y - gy[m20]) / dy;
				rz = (z - gz[m30]) / dz;
				/* Normalized distance of visualization node to P reference node */
				rxc = (x - (gx[m10c - 1] + gx[m10c]) / 2.0) / dxc;
				ryc = (y - (gy[m20c - 1] + gy[m20c]) / 2.0) / dyc;
				rzc = (z - (gz[m30c - 1] + gz[m30c]) / 2.0) / dzc;
				/* Normalized distance of visualization node to Vx,Vy,Vz reference node */
				rxv = (x - (gx[m10v] + gx[m10v + 1]) / 2.0) / dxv;
				ryv = (y - (gy[m20v] + gy[m20v + 1]) / 2.0) / dyv;
				rzv = (z - (gz[m30v] + gz[m30v + 1]) / 2.0) / dzv;
				/* Define reference node */
				m40 = m30 * xnumx * ynumy + m10 * ynumy + m20;
				m40pr = m30c * xnumx * ynumy + m10c * ynumy + m20c;
				m40vx = m30v * xnumx * ynumy + m10 * ynumy + m20v;
				m40vy = m30v * xnumx * ynumy + m10v * ynumy + m20;
				m40vz = m30 * xnumx * ynumy + m10v * ynumy + m20v;
				/**/
				/* Temperature */
				tk_vtk[m4] = (1.0 - rx) * (1.0 - ry) * (1.0 - rz) * tk[m40];
				tk_vtk[m4] += (1.0 - rx) * (ry) * (1.0 - rz) * tk[m40 + 1];
				tk_vtk[m4] += (rx) * (1.0 - ry) * (1.0 - rz) * tk[m40 + ynumy];
				tk_vtk[m4] += (rx) * (ry) * (1.0 - rz) * tk[m40 + ynumy + 1];
				;
				tk_vtk[m4] += (1.0 - rx) * (1.0 - ry) * (rz) * tk[m40 + xnumx * ynumy];
				tk_vtk[m4] += (1.0 - rx) * (ry) * (rz) * tk[m40 + xnumx * ynumy + 1];
				tk_vtk[m4] += (rx) * (1.0 - ry) * (rz) * tk[m40 + xnumx * ynumy + ynumy];
				tk_vtk[m4] += (rx) * (ry) * (rz) * tk[m40 + xnumx * ynumy + 1 + ynumy];
				/* Viscosity */
				nu_vtk[m4] = (1.0 - rx) * (1.0 - ry) * (1.0 - rz) * nu[m40];
				nu_vtk[m4] += (1.0 - rx) * (ry) * (1.0 - rz) * nu[m40 + 1];
				nu_vtk[m4] += (rx) * (1.0 - ry) * (1.0 - rz) * nu[m40 + ynumy];
				nu_vtk[m4] += (rx) * (ry) * (1.0 - rz) * nu[m40 + ynumy + 1];
				;
				nu_vtk[m4] += (1.0 - rx) * (1.0 - ry) * (rz) * nu[m40 + xnumx * ynumy];
				nu_vtk[m4] += (1.0 - rx) * (ry) * (rz) * nu[m40 + xnumx * ynumy + 1];
				nu_vtk[m4] += (rx) * (1.0 - ry) * (rz) * nu[m40 + xnumx * ynumy + ynumy];
				nu_vtk[m4] += (rx) * (ry) * (rz) * nu[m40 + xnumx * ynumy + 1 + ynumy];
				/* gravity potential */
				//gp_vtk[m4]=(1.0-rx)*(1.0-ry)*(1.0-rz)*gp[m40];
				//gp_vtk[m4]+=(1.0-rx)*(ry)*(1.0-rz)*gp[m40+1];
				//gp_vtk[m4]+=(rx)*(1.0-ry)*(1.0-rz)*gp[m40+ynumy];
				//gp_vtk[m4]+=(rx)*(ry)*(1.0-rz)*gp[m40+ynumy+1];;
				//gp_vtk[m4]+=(1.0-rx)*(1.0-ry)*(rz)*gp[m40+xnumx*ynumy];
				//gp_vtk[m4]+=(1.0-rx)*(ry)*(rz)*gp[m40+xnumx*ynumy+1];
				//gp_vtk[m4]+=(rx)*(1.0-ry)*(rz)*gp[m40+xnumx*ynumy+ynumy];
				//gp_vtk[m4]+=(rx)*(ry)*(rz)*gp[m40+xnumx*ynumy+1+ynumy];
				/* Density */
				ro_vtk[m4] = (1.0 - rx) * (1.0 - ry) * (1.0 - rz) * ro[m40];
				ro_vtk[m4] += (1.0 - rx) * (ry) * (1.0 - rz) * ro[m40 + 1];
				ro_vtk[m4] += (rx) * (1.0 - ry) * (1.0 - rz) * ro[m40 + ynumy];
				ro_vtk[m4] += (rx) * (ry) * (1.0 - rz) * ro[m40 + ynumy + 1];
				;
				ro_vtk[m4] += (1.0 - rx) * (1.0 - ry) * (rz) * ro[m40 + xnumx * ynumy];
				ro_vtk[m4] += (1.0 - rx) * (ry) * (rz) * ro[m40 + xnumx * ynumy + 1];
				ro_vtk[m4] += (rx) * (1.0 - ry) * (rz) * ro[m40 + xnumx * ynumy + ynumy];
				ro_vtk[m4] += (rx) * (ry) * (rz) * ro[m40 + xnumx * ynumy + ynumy + 1];
				/* Pressure */
				/* Solid pressure is not defined at nodes, but at cell center */
				pr_vtk[m4] = (1.0 - rxc) * (1.0 - ryc) * (1.0 - rzc) * pr[m40pr];
				pr_vtk[m4] += (1.0 - rxc) * (ryc) * (1.0 - rzc) * pr[m40pr + 1];
				pr_vtk[m4] += (rxc) * (1.0 - ryc) * (1.0 - rzc) * pr[m40pr + ynumy];
				pr_vtk[m4] += (rxc) * (ryc) * (1.0 - rzc) * pr[m40pr + ynumy + 1];
				pr_vtk[m4] += (1.0 - rxc) * (1.0 - ryc) * (rzc) * pr[m40pr + ynumy * xnumx];
				pr_vtk[m4] += (1.0 - rxc) * (ryc) * (rzc) * pr[m40pr + ynumy * xnumx + 1];
				pr_vtk[m4] += (rxc) * (1.0 - ryc) * (rzc) * pr[m40pr + ynumy * xnumx + ynumy];
				pr_vtk[m4] += (rxc) * (ryc) * (rzc) * pr[m40pr + ynumy * xnumx + ynumy + 1];
				/*  Vx computation */
				/* Vx is not defined at nodes, but at left and right cell side center */
				vx_vtk[m4] = (1.0 - rx) * (1.0 - ryv) * (1.0 - rzv) * vx[m40vx];
				vx_vtk[m4] += (1.0 - rx) * (ryv) * (1.0 - rzv) * vx[m40vx + 1];
				vx_vtk[m4] += (rx) * (1.0 - ryv) * (1.0 - rzv) * vx[m40vx + ynumy];
				vx_vtk[m4] += (rx) * (ryv) * (1.0 - rzv) * vx[m40vx + ynumy + 1];
				vx_vtk[m4] += (1.0 - rx) * (1.0 - ryv) * (rzv) * vx[m40vx + ynumy * xnumx];
				vx_vtk[m4] += (1.0 - rx) * (ryv) * (rzv) * vx[m40vx + ynumy * xnumx + 1];
				vx_vtk[m4] += (rx) * (1.0 - ryv) * (rzv) * vx[m40vx + ynumy * xnumx + ynumy];
				vx_vtk[m4] += (rx) * (ryv) * (rzv) * vx[m40vx + ynumy * xnumx + ynumy + 1];
				/*  Vy computation */
				/* Vy is not defined at nodes, but at top and bottom cell side center */
				vy_vtk[m4] = (1.0 - rxv) * (1.0 - ry) * (1.0 - rzv) * vy[m40vy];
				vy_vtk[m4] += (1.0 - rxv) * (ry) * (1.0 - rzv) * vy[m40vy + 1];
				vy_vtk[m4] += (rxv) * (1.0 - ry) * (1.0 - rzv) * vy[m40vy + ynumy];
				vy_vtk[m4] += (rxv) * (ry) * (1.0 - rzv) * vy[m40vy + ynumy + 1];
				vy_vtk[m4] += (1.0 - rxv) * (1.0 - ry) * (rzv) * vy[m40vy + ynumy * xnumx];
				vy_vtk[m4] += (1.0 - rxv) * (ry) * (rzv) * vy[m40vy + ynumy * xnumx + 1];
				vy_vtk[m4] += (rxv) * (1.0 - ry) * (rzv) * vy[m40vy + ynumy * xnumx + ynumy];
				vy_vtk[m4] += (rxv) * (ry) * (rzv) * vy[m40vy + ynumy * xnumx + ynumy + 1];
				/*  Vz computation */
				/* Vz is not defined at nodes, but at front and back cell side center */
				vz_vtk[m4] = (1.0 - rxv) * (1.0 - ryv) * (1.0 - rz) * vz[m40vz];
				vz_vtk[m4] += (1.0 - rxv) * (ryv) * (1.0 - rz) * vz[m40vz + 1];
				vz_vtk[m4] += (rxv) * (1.0 - ryv) * (1.0 - rz) * vz[m40vz + ynumy];
				vz_vtk[m4] += (rxv) * (ryv) * (1.0 - rz) * vz[m40vz + ynumy + 1];
				vz_vtk[m4] += (1.0 - rxv) * (1.0 - ryv) * (rz) * vz[m40vz + ynumy * xnumx];
				vz_vtk[m4] += (1.0 - rxv) * (ryv) * (rz) * vz[m40vz + ynumy * xnumx + 1];
				vz_vtk[m4] += (rxv) * (1.0 - ryv) * (rz) * vz[m40vz + ynumy * xnumx + ynumy];
				vz_vtk[m4] += (rxv) * (ryv) * (rz) * vz[m40vz + ynumy * xnumx + ynumy + 1];
			}
}