Allink  v0.1
ElPolyOutput.cpp
1 #include "ElPoly.h"
2 
3 void ElPoly::Conv2Tecplot(int NBin,int How){
4  int NType = 3;
5  double Round = 0.001;
6  int Nx = MIN(NBin,NEdge);
7  int Ny = MIN(NBin,NEdge);
8  double **Plot = (double **)calloc(NType,sizeof(double));
9  double **Count = (double **)calloc(NType,sizeof(double));
10  for(int t=0;t<NType;t++){
11  Plot[t] = (double *)calloc(Nx*Ny,sizeof(double));
12  Count[t] = (double *)calloc(Nx*Ny,sizeof(double));
13  }
14  FILE *TecPlot = fopen("TecPlot.dat","w");
15  if(How == 0){//Dens
16  fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"oil density\",\"density phob\",\"density phil\"\n");
17  fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
18  }
19  else if(How == 1){
20  fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"pressure\",\"density phob\",\"density phil\"\n");
21  fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
22  }
23  else if(How == 2){
24  fprintf(TecPlot,"VARIABLES = \"R\", \"Z\", \"height\",\"density phob\",\"density phil\"\n");
25  fprintf(TecPlot,"ZONE J=%d, K=%d, F=POINT\n",Nx,Ny);
26  }
27  for(int p=0;p<pNPart();p++){
28  int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*Nx);
29  if(vx < 0 || vx >= Nx) continue;
30  int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*Ny);
31  if(vy < 0 || vy >= Ny) continue;
32  if(How == 0){//Dens
33  int t = Pm[p].Typ;
34  t = (t+1)%3;
35  // if(t==2){t=0;Plot[1][vx*Nx+vy] += -1.;}
36  // else
37  Plot[t][vx*Nx+vy] += pPos(p,2);
38  Count[t][vx*Nx+vy] += 1.;
39  }
40  else if(How == 1){
41  for(int t=0;t<3;t++){
42  Plot[t][vx*Nx+vy] += pVel(p,t);
43  Count[t][vx*Nx+vy] += 1.;
44  }
45  }
46  else if(How == 2){
47  int t = 0;//Pm[p].Typ;
48  Plot[t][vx*Nx+vy] += pPos(p,2);
49  Count[t][vx*Nx+vy] += 1.;
50  }
51  }
52  Matrice Mask(5,5);
53  Mask.FillGaussian(.5,3.);
54  int NDim = 2;
55  int IfMinImConv = 0;
56  for(int t=0;t<3;t++){
57  //Mask.ConvoluteMatrix(Plot[t],Nx,2);
58  Mask.ConvoluteMatrix(Plot[t],Nx,NDim,IfMinImConv);
59  }
60  for(int vx=0;vx<Nx;vx++){
61  for(int vy=0;vy<Ny;vy++){
62  double r = vx*pEdge(0)/(double)Nx;
63  double z = vy*pEdge(1)/(double)Ny - pEdge(1)*.5;
64  double Norm0 = Count[0][vx*Nx+vy] > 0. ? 1./Count[0][vx*Nx+vy] : 1.;
65  double Norm1 = Count[1][vx*Nx+vy] > 0. ? 1./Count[1][vx*Nx+vy] : 1.;
66  double Norm2 = Count[2][vx*Nx+vy] > 0. ? 1./Count[2][vx*Nx+vy] : 1.;
67  fprintf(TecPlot,"%lf %lf %lf %lf %lf\n",r,z,
68  Plot[0][vx*Nx+vy]*Norm0,Plot[1][vx*Nx+vy]*Norm1,Plot[2][vx*Nx+vy]*Norm2);
69  }
70  }
71  // if(Pm[p].Vel[0]+Pm[p].Vel[1]+Pm[p].Vel[2]<.1){
72  // Dens1 = -5.;
73  // Dens2 = -5.;
74  // Dens3 = -5.;
75  // }
76  // if(pPos(p,0) > 0 && pPos(p,0) < 5.)
77  // if(pPos(p,1)>13. && pPos(p,1)<19.)
78  // if(Pm[p].Vel[0] < 0.01){
79  // Dens1 = -3.;
80  // Dens2 = -3.;
81  // Dens3 = -3.;
82  // }
83  fclose(TecPlot);
84  for(int t=0;t<NType;t++){
85  free(Plot[t]);
86  free(Count[t]);
87  }
88  free(Plot);
89  free(Count);
90 }
92  FILE *OutVmd = fopen("Sim.vtf","w");
93  int cOff = 0;
94  int pOff = 0;
95  char Shape[30];
96  int NanoPoint = 10;
97  for(int b=0;b<pNBlock();b++){
98  int bType = 0;
99  if(!strncmp(Block[b].Name,"PEP",3)){
100  bType = 2;
101  }
102  for(int c=cOff;c<cOff+pNChain(b);c++){
103  for(int p=pOff;p<pOff+pNPCh(b);p++){
104  int Type = pType(p) + bType;
105  if(Type == 0)
106  fprintf(OutVmd,"a %d r 0.8 n A resid %d res %s\n",p,pChain(p),Block[b].Name);
107  else if(Type == 1)
108  fprintf(OutVmd,"a %d r 0.8 n B resid %d res %s\n",p,pChain(p),Block[b].Name);
109  else if(Type == 2)
110  fprintf(OutVmd,"a %d r 0.8 n D resid %d res %s\n",p,pChain(p),Block[b].Name);
111  else if(Type == 3)
112  fprintf(OutVmd,"a %d r 0.8 n E resid %d res %s\n",p,pChain(p),Block[b].Name);
113  }
114  pOff += pNPCh(b);
115  fprintf(OutVmd,"b %d::%d\n",c*pNPCh(b),(c+1)*pNPCh(b)-1);
116  }
117  cOff += pNChain(b);
118  }
119  for(int n=0;n<pNNano();n++){
120  ShapeId(Nano[n].Shape,Shape);
121  for(int i=0;i<NanoPoint;i++){
122  fprintf(OutVmd,"a %d r %lf n Nano resid %d res %s\n",n*NanoPoint+pNPart()+i,Nano[n].Rad,cOff+n,Shape);
123  }
124  fprintf(OutVmd,"b %d::%d\n",n*NanoPoint+pNPart(),n*NanoPoint+pNPart()+NanoPoint-1);
125  }
126  for(int f=NFile[0];f<NFile[1];f++){
127  Processing(f);
128  if(OpenRisk(cFile[f],BF_CHAIN)) return ;
129  fprintf(OutVmd,"\ntimestep\npbc %lf %lf %lf\n",pEdge(0),pEdge(1),pEdge(2));
130  for(int p=0;p<pNPart();p++){
131  fprintf(OutVmd,"%lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2));
132  }
133  for(int n=0;n<pNNano();n++){
134  double Pos[3];
135  for(int d=0;d<3;d++){
136  Pos[d] = pNanoPos(n,d) - .5*Nano[n].Height*Nano[n].Axis[d];
137  }
138  for(int i=0;i<NanoPoint;i++){
139  for(int d=0;d<3;d++){
140  Pos[d] += Nano[n].Height/(double)NanoPoint*Nano[n].Axis[d];
141  }
142  fprintf(OutVmd,"%lf %lf %lf\n",Pos[0],Pos[1],Pos[2]);
143  }
144  }
145  }
146 }
147 //-----------------------POVRAY--------------------------
149  double Pos[3];
150  for(int d=0;d<3;d++) Pos[d] = pPos(p,d)*InvScaleUn;
151  int Typ = pType(p) < 6 ? pType(p) : 5;
152  int Chc = Ch[Pm[p].CId].Type;
153  if(pType(p) == 0 && CHAIN_IF_TYPE(Chc,CHAIN_ADDED) )Typ = 3;
154  fprintf(DrawOutFile,"ElSphere(<%.4f, %.4f, %.4f>,",Pos[0],Pos[1],Pos[2]);
155  fprintf(DrawOutFile,"SphRad,rgbt<%.4f, %.4f, %.4f, %.4f>)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
156 }
158  if(VAR_IF_TYPE(Nano[n].Shape,SHAPE_CYL)){
159  double PosP[3];
160  double PosN[3];
161  int Typ = 2;
162  for(int d=0;d<3;d++){
163  PosP[d] = pNanoPos(n,d) + .5*Nano[n].Height*Nano[n].Axis[d];
164  PosN[d] = pNanoPos(n,d) - .5*Nano[n].Height*Nano[n].Axis[d];
165  PosP[d] *= InvScaleUn;
166  PosN[d] *= InvScaleUn;
167  }
168  fprintf(DrawOutFile,"ElCylinder(<%.2f, %.2f, %.2f>,",PosP[0],PosP[1],PosP[2]);
169  fprintf(DrawOutFile,"<%.2f, %.2f, %.2f>%.2f,",PosN[0],PosN[1],PosN[2],Nano[n].Rad*InvScaleUn);
170  fprintf(DrawOutFile,"rgbt<%.4f, %.4f, %.4f, %.4f>,1)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
171  }
172  else{
173  // Point2Shape(Nano[n].Shape);
174  // DrField(128,SQR(Nano[n].Rad),n,DrawOutFile);
175  }
176 }
177 void ElPoly::DrBondPovRay(double *Pos1,double *Pos2,float *Color){
178  fprintf(DrawOutFile,"ElCylinder(<%lf, %lf, %lf>,",Pos1[0],Pos1[1],Pos1[2]);
179  fprintf(DrawOutFile,"<%lf, %lf, %lf>",Pos2[0],Pos2[1],Pos2[2]);
180  fprintf(DrawOutFile,"CylRad,rgbt<%.4f, %.4f, %.4f, %.4f>,1)\n",Color[0],Color[1],Color[2],Color[3]-1.);
181 }
183  //particles, lines
184  HeaderPovRay();
185  for(int f=NFile[0];f<NFile[1];f++){
186  SigErr(DrawOutFile != NULL,"DrawOutFile already allocated, can't use the file");
187  char FName[60];
188  sprintf(FName,"PovSnap%05d.pov",pStep());
189  DrawOutFile = fopen(FName,"w");
190  int ImSize[2] = {1000,1000};
191  fprintf(DrawOutFile,"// POV 3.x input script : plot.pov\n// command: povray +W%d +H%d -I%s -O%s.tga +P +X +A +FT +C\n",ImSize[0],ImSize[1],FName,FName);
192  fprintf(DrawOutFile,"#if (version < 3.5)\n#error \"POV3DisplayDevice has been compiled for POV-Ray 3.5 or above.\\nPlease upgrade POV-Ray.\"\n#end\n");
193  fprintf(DrawOutFile,"#include \"PovHeader.inc\"\n");
194  fprintf(DrawOutFile,"// System \n");
195  for(int b=0,NPep=0;b<pNBlock();b++){
196  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
197  DrPartPovRay(p);
198  int Typ = Pm[p].Typ;
199  double Pos1[3];
200  double Pos2[3];
201  for(int l=0;l<Ln[p].NLink;l++){
202  int link = Ln[p].Link[l];
203  if(p == link) continue;
204  for(int d=0;d<3;d++){
205  Pos1[d] = pPos(p,d)*InvScaleUn;
206  Pos2[d] = pPos(p,link)*InvScaleUn;
207  if(Pos1[d] - Pos2[d] > .5*InvScaleUn)
208  Pos2[d] += pEdge(d)*InvScaleUn;
209  else if(Pos1[d] - Pos2[d] < -.5*InvScaleUn)
210  Pos2[d] -= pEdge(d)*InvScaleUn;
211  }
212  DrBondPovRay(Pos1,Pos2,ColorType[Typ]);
213  }
214  }
215  }
216  for(int n=0;n<pNNano();n++){
217  DrNanoPovRay(n);
218  }
219  fclose(DrawOutFile);
220  }
221 }
222 #include "ElPolyDrawSurf.h"
223 void ElPoly::DrField(int NGrid,double IsoLevel,int nNano,FILE *FWrite){
224  int Typ = 2;
225  double CubeDist[8];
226  double EdgeVertex[12][3];
227  double EdgeNormal[12][3];
228  double InvNGrid = 1./(double)NGrid;
229  double Pos[3];
230  double Pos1[3];
231  double Cm[3];
232  for(int d=0;d<3;d++){
233  Cm[d] = .5*pEdge(d)*InvScaleUn;
234  }
235  for(int gx=0;gx<NGrid;gx++){
236  Pos[0] = gx*InvNGrid*pEdge(0);
237  for(int gy=0;gy<NGrid;gy++){
238  Pos[1] = gy*InvNGrid*pEdge(1);
239  for(int gz=0;gz<NGrid;gz++){
240  Pos[2] = gz*InvNGrid*pEdge(2);
241  for(int v=0;v<8;v++){
242  for(int d=0;d<3;d++){
243  Pos1[d] = Pos[d] + VertCube[v][d]*InvNGrid*pEdge(d);
244  }
245  CubeDist[v] = NanoDist2(Pos1,nNano);
246  }
247  int Flag = 0;
248  for(int v=0;v<8;v++){
249  if(CubeDist[v] <= IsoLevel)
250  Flag |= 1<<v;
251  }
252  int CubeShape = CubeTop[Flag];
253  if(CubeShape==0) continue;
254  for(int e=0;e<12;e++){
255  if(CubeShape & (1<<e)){
256  double Delta = CubeDist[EdgeConn[e][1]] - CubeDist[EdgeConn[e][0]];
257  double OffSet = (IsoLevel-CubeDist[EdgeConn[e][0]])/Delta;
258  if(Delta == 0.0){
259  OffSet = .5;
260  }
261  EdgeNormal[e][0] = NanoDist2(Pos[0]-0.01,Pos[1],Pos[2],nNano)
262  - NanoDist2(Pos[0]+0.01,Pos[1],Pos[2],nNano);
263  EdgeNormal[e][1] = NanoDist2(Pos[0],Pos[1]-0.01,Pos[2],nNano)
264  - NanoDist2(Pos[0],Pos[1]+0.01,Pos[2],nNano);
265  EdgeNormal[e][2] = NanoDist2(Pos[0],Pos[1],Pos[2]-0.01,nNano)
266  - NanoDist2(Pos[0],Pos[1],Pos[2]+0.01,nNano);
267  double Norm = sqrt(SQR(EdgeNormal[e][0])+SQR(EdgeNormal[e][1])+SQR(EdgeNormal[e][2]));
268  for(int d=0;d<3;d++){
269  EdgeVertex[e][d] = Pos[d] + (VertCube[EdgeConn[e][0]][d]+OffSet*EdgeDir[e][d])*InvNGrid*pEdge(d);
270  EdgeVertex[e][d] *= InvScaleUn;
271  EdgeNormal[e][d] *= 1./Norm;
272  }
273  }
274  }
275  for(int t=0;t<5;t++){
276  if(TrConnTable[Flag][3*t] < 0.) break;
277  fprintf(FWrite,"ElTriangle (");
278  for(int d=0;d<3;d++){
279  int v = TrConnTable[Flag][3*t+d];
280  fprintf(FWrite,"<%.2f, %.2f, %.2f>,",EdgeVertex[v][0]-Cm[0],EdgeVertex[v][1]-Cm[1],EdgeVertex[v][2]-Cm[2]);
281  fprintf(FWrite,"<%.2f, %.2f, %.2f>,",EdgeNormal[v][0],EdgeNormal[v][1],EdgeNormal[v][2]);
282  }
283  fprintf(FWrite,"rgbt<%.3f, %.3f, %.3f, %.3f>)\n",ColorType[Typ][0],ColorType[Typ][1],ColorType[Typ][2],ColorType[Typ][3]-1.);
284  }
285  }
286  }
287  }
288 }
289 void ElPoly::Conv2rzd(int NBin){
290  // FILE *FOut = fopen("data.dat","w");
291  // for(int p=0;p<pNPart();p++){
292  // fprintf(FOut,"%lf %lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0));
293  // }
294  // fclose(FOut);
295  // return;
296  int NType = 3;
297  double Round = 0.001;
298  double **Plot = (double **)calloc(NType,sizeof(double));
299  for(int t=0;t<NType;t++)
300  Plot[t] = (double *)calloc(NBin*NBin,sizeof(double));
301  double *Count = (double *)calloc(NBin*NBin*NType,sizeof(double));
302  FILE *TecPlot = fopen("data.dat","w");
303  for(int p=0;p<pNPart();p++){
304  int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*NBin);
305  if(vx < 0 || vx >= NBin) continue;
306  int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*NBin);
307  if(vy < 0 || vy >= NBin) continue;
308  for(int t=0;t<3;t++){
309  Plot[t][vx*NBin+vy] += pVel(p,t);
310  Count[(vx*NBin+vy)*NType+t] += 1.;
311  }
312  }
313  //smooth and write dens
314  Matrice Mask(5,5);
315  Mask.FillGaussian(.5,3.);
316  Mask.Print();
317  int NDim = 2;
318  int IfMinImConv = 0;
319  for(int t=0;t<3;t++){
320  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
321  Mask.ConvoluteMatrix(Plot[t],NBin,NDim,IfMinImConv);
322  }
323  for(int vx=0;vx<NBin;vx++){
324  for(int vy=0;vy<NBin;vy++){
325  double r = vx*pEdge(0)/(double)NBin;
326  double z = vy*pEdge(1)/(double)NBin - pEdge(1)*.5;
327  double Norm0 = Count[(vx*NBin+vy)*NType+0] > 0. ? 1./Count[(vx*NBin+vy)*NType+0] : 1.;
328  double Norm1 = Count[(vx*NBin+vy)*NType+1] > 0. ? 1./Count[(vx*NBin+vy)*NType+1] : 1.;
329  double Norm2 = Count[(vx*NBin+vy)*NType+2] > 0. ? 1./Count[(vx*NBin+vy)*NType+2] : 1.;
330  fprintf(TecPlot,"%lf %lf %lf %lf %lf\n",r,z,
331  Plot[0][vx*NBin+vy]*Norm0,Plot[1][vx*NBin+vy]*Norm1,Plot[2][vx*NBin+vy]*Norm2);
332  // fprintf(TecPlot,"%lf %lf %lf \n",r,z,Plot[1][vx*NBin+vy]*Norm0);
333  }
334  }
335  fclose(TecPlot);
336  for(int t=0;t<NType;t++)
337  free(Plot[t]);
338  free(Plot);
339  free(Count);
340 }
341 void ElPoly::Conv2xyzd(int NBin){
342  // FILE *FOut = fopen("data.dat","w");
343  // for(int p=0;p<pNPart();p++){
344  // fprintf(FOut,"%lf %lf %lf %lf\n",pPos(p,0),pPos(p,1),pPos(p,2),pVel(p,0));
345  // }
346  // fclose(FOut);
347  // return;
348  int NType = 3;
349  double Round = 0.001;
350  double **Plot = (double **)calloc(NType,sizeof(double));
351  for(int t=0;t<NType;t++)
352  Plot[t] = (double *)calloc(NBin*NBin*NBin,sizeof(double));
353  double *Count = (double *)calloc(NBin*NBin*NBin*NType,sizeof(double));
354  FILE *TecPlot = fopen("data.dat","w");
355  for(int p=0;p<pNPart();p++){
356  int vx = (int)((pPos(p,0)+Round)*pInvEdge(0)*NBin);
357  if(vx < 0 || vx >= NBin) continue;
358  int vy = (int)((pPos(p,1)+Round)*pInvEdge(1)*NBin);
359  if(vy < 0 || vy >= NBin) continue;
360  int vz = (int)((pPos(p,2)+Round)*pInvEdge(2)*NBin);
361  if(vz < 0 || vz >= NBin) continue;
362  for(int t=0;t<3;t++){
363  Plot[t][(vx*NBin+vy)*NBin+vz] += pVel(p,t);
364  Count[((vx*NBin+vy)*NBin+vz)*NType+t] += 1.;
365  }
366  }
367  //smooth and write dens
368  for(int vx=0;vx<NBin;vx++){
369  for(int vy=0;vy<NBin;vy++){
370  for(int vz=0;vz<NBin;vz++){
371  double x = vx*pEdge(0)/(double)NBin;
372  double y = vy*pEdge(1)/(double)NBin;
373  double z = vz*pEdge(1)/(double)NBin;
374  double Norm0 = Count[((vx*NBin+vy)*NBin+vz)*NType+0] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+0] : 1.;
375  double Norm1 = Count[((vx*NBin+vy)*NBin+vz)*NType+1] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+1] : 1.;
376  double Norm2 = Count[((vx*NBin+vy)*NBin+vz)*NType+2] > 0. ? 1./Count[((vx*NBin+vy)*NBin+vz)*NType+2] : 1.;
377  fprintf(TecPlot,"%lf %lf %lf %lf %lf %lf\n",x,y,z,
378  Plot[0][(vx*NBin+vy)*NBin+vz]*Norm0,Plot[1][(vx*NBin+vy)*NBin+vz]*Norm1,Plot[2][(vx*NBin+vy)*NBin+vz]*Norm2);
379  // fprintf(TecPlot,"%lf %lf %lf %lf\n",x,y,z,Plot[0][vx*NBin+vy]*Norm0);
380  }
381  }
382  }
383  fclose(TecPlot);
384  for(int t=0;t<NType;t++)
385  free(Plot[t]);
386  free(Plot);
387  free(Count);
388 }
390  FILE *HeaderPov = fopen("PovHeader.inc","w");
391  double SphRad = 0.005;
392  double CylRad = 0.002;
393  double PerspAngle = 15;
394  int ImSize[2] = {1000,1000};
395  fprintf(HeaderPov,"// POV 3.x input script : plot.pov\n// command: povray +W%d +H%d -ISnap.pov -OSnap.pov.tga +P +X +A +FT +C\n",ImSize[0],ImSize[1]);
396  fprintf(HeaderPov,"#if (version < 3.5)\n#error \"POV3DisplayDevice has been compiled for POV-Ray 3.5 or above.\\nPlease upgrade POV-Ray.\"\n#end\n");
397  //include
398  fprintf(HeaderPov,"// #include \"colors.inc\"\n \
399 // #include \"stones.inc\"\n\
400 // #include \"textures.inc\"\n\
401 // #include \"shapes.inc\"\n\
402 // #include \"glass.inc\"\n\
403 // #include \"metals.inc\"\n\
404 // #include \"woods.inc\"\n");
405  fprintf(HeaderPov,"\
406 #declare clip_on=array[3] {0, 0, 0};\n\
407 #declare clip=array[3];\n\
408 #declare scaledclip=array[3];\n\
409 #declare SphRad=%lf;\n\
410 #declare CylRad=%lf;\n\
411 #declare line_width=0.0020;\n",SphRad,CylRad);
412  // //materials
413  // fprintf(HeaderPov,"\
414  // #declare Fotoni = photons {\n\
415  // target\n\
416  // refraction on\n\
417  // reflection on\n\
418  // }\n");
419  //define a point
420  fprintf(HeaderPov,"\
421 #macro ElPoint (P1, R1, C1)\n\
422  #local T = texture { finish { ambient 1.0 diffuse 0.0 phong 0.0 specular 0.0 } pigment { C1 } }\n \
423  #if(clip_on[2])\n\
424  intersection {\n\
425  sphere {P1, R1 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
426  clip[2]\n\
427  }\n\
428  #else\n\
429  sphere {P1, R1 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
430  #end\n\
431 #end\n");
432  //define a line
433  fprintf(HeaderPov,"\
434 #macro ElLine (P1, P2, C1)\n\
435  #local T = texture { finish { ambient 1.0 diffuse 0.0 phong 0.0 specular 0.0 } pigment { C1 } }\n\
436  #if(clip_on[2])\n\
437  intersection {\n\
438  cylinder {P1, P2, line_width texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
439  clip[2]\n\
440  }\n\
441  #else\n\
442  cylinder {P1, P2, line_width texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
443  #end\n\
444 #end\n");
445  //define a sphere
446  fprintf(HeaderPov,"#macro ElSphere(P1, R1, C1)\n\
447  #local T = texture { pigment { C1 } }\n\
448  #local M = material{\n \
449  texture{\n\
450  pigment{gradient y color_map{[0.4 C1][0.4 C1]}}\n\
451  finish{ambient 0 diffuse 0.4 specular 1 roughness 0.0001 reflection 0.25}\n\
452  }\n\
453  interior{ior 1.33}\n \
454  }\n\
455  #if(clip_on[2])\n\
456  intersection {\n\
457  sphere {P1, R1 material {M} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
458  clip[2]\n\
459  }\n\
460  #else\n\
461  sphere {P1, R1 material {M} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
462  #end\n\
463 #end\n");
464  //define a cylinder
465  fprintf(HeaderPov,"\
466 #macro ElCylinder (P1, P2, R1, C1, O1)\n\
467  #local T = texture { pigment { C1 } }\n\
468  #if(clip_on[2])\n\
469  intersection {\n\
470  cylinder {P1, P2, R1 #if(O1) open #end texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow\n\
471 }\n\
472  clip[2]\n\
473  }\n\
474  #else\n\
475  cylinder {P1, P2, R1 #if(O1) open #end texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
476  #end\n\
477 #end\n");
478  //define a triangle
479  fprintf(HeaderPov,"\
480 #macro ElTriangle (P1, N1, P2, N2, P3, N3, C1)\n\
481  #local T = texture { pigment { C1 } }\n\
482  smooth_triangle {P1, N1, P2, N2, P3, N3 texture {T} #if(clip_on[1]) clipped_by {clip[1]} #end no_shadow}\n\
483 #end\n");
484  //set the camera
485  double Loc[3] = {0.5,4.0,4.};
486  double LAt[3] = {0.5,0.5,0.5};
487  double Up[3] = {0.,6.,0.};
488  double Right[3] = {4.8,0.,0.};
489  double Dir[3] = {0.,0.,4.};
490  double Light0[3] = {0.5,4.0,4.0};
491  double Light1[3] = {0.5,4.0,2.0};
492  double CBack[3] = {1.,1.,1.};
493  #ifdef __glut_h__
494  // Draw *Dr;
495  // Loc[0] = Dr->xp; Loc[1] = Dr->yp; Loc[2] = Dr->zp+Dr->zw;
496  // Light0[0] = Dr->xl0; Light0[1] = Dr->yl0; Light0[2] = Dr->zl0;
497  // Light1[0] = Dr->xl1; Light1[1] = Dr->yl1; Light1[2] = Dr->zl1;
498  // CBack[0] = Dr->Rback; CBack[1] = Dr->Gback; CBack[2] = Dr->Bback;
499  #endif //__glut_h__
500  fprintf(HeaderPov,"camera {\n\
501 //perspective orthographic fisheye ultra_wide_angle...\n\
502  location <%.3f, %.3f, %.3f>\n\
503  look_at <%.3f, %.3f, %.3f>\n\
504  sky <0. 0. 1.>\n\
505 // up y\n\
506 // right 1.33*x\n\
507 // direction <%.3f, %.3f, %.3f>\n\
508 // translate <0.0, 1.0, 1.0>\n\
509 // rotate <0.0, 1.0, 1.0>\n\
510  angle %f\n \
511 }\n",
512  Loc[0],Loc[1],Loc[2],
513  LAt[0],LAt[1],LAt[2],
514  Dir[0],Dir[1],Dir[2],
515  PerspAngle);
516  //set the light0
517  fprintf(HeaderPov,"\
518 light_source {\n\
519  <%lf, %lf, %lf>\n\
520  color rgb<1.000, 1.000, 1.000>\n\
521  parallel\n\
522  point_at <.5, .5, .5>\n\
523 }\n",Light0[0],Light0[1],Light0[2]);
524  //set the light1
525  fprintf(HeaderPov,"\
526 light_source {\n\
527  <%lf, %lf, %lf>\n\
528  color rgb<1.000, 1.000, 1.000>\n\
529  parallel\n\
530  point_at <.5, .5, .5>\n\
531 }\n",Light1[0],Light1[1],Light1[2]);
532  //set the fog
533  fprintf(HeaderPov,"\
534  fog { fog_type 2\n\
535  distance 50\n\
536  color rgb<1.000, 1.000, 1.000>\n\
537  fog_offset 0.1\n\
538  fog_alt 1.5\n\
539  turbulence 1.8\n\
540  }\n");
541  //background
542  fprintf(HeaderPov,"\
543 background {\n\
544  color rgb<%lf, %lf, %lf>\n\
545 }\n",CBack[0],CBack[1],CBack[2]);
546  //texture
547  fprintf(HeaderPov,"\
548 #default { texture {\n\
549  finish { ambient 0.000 diffuse 0.650 phong 0.1 phong_size 40.000 specular 0.500 }\n\
550 } }\n");
551 }
552 void ElPoly::ConvLattice(int NSample,char *FName){
553  int NType = 3;
554  double **Plot = (double **) calloc(3,sizeof(double));
555  for(int t=0;t<NType;t++){
556  Plot[t] = (double *)calloc(NSample*NSample,sizeof(double));
557  }
558  LoadDensFile(Plot,NSample);
559  FILE *FWrite = fopen(FName,"w");
560  fprintf(FWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(0),pEdge(1),pEdge(2),NSample,ChooseDraw(EL_PART));
561  double InvNSample = 1./(double)NSample;
562  for(int t=0;t<3;t++){
563  for(int sx=0;sx<NSample;sx++){
564  for(int sy=0;sy<NSample;sy++){
565  //if(fabs(Plot[t][sx*NSample+sy]) < 0.01) continue;
566  double x = sx*InvNSample*pEdge(CLat1);
567  double y = sy*InvNSample*pEdge(CLat2);
568  fprintf(FWrite,"{t[%d %d %d] x(%lf %lf %lf)}\n",t*pNPart()+sx*NSample+sy,t,t,x,y,Plot[t][sx*NSample+sy]);
569  }
570  }
571  }
572  for(int t=0;t<3;t++) free(Plot[t]);
573  free(Plot);
574 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int CId
Chain Identifier.
Definition: VarData.h:224
void Conv2Povray()
esport the data in pov file format for rendering
void Conv2xyzd(int NSample)
esport the data in radius depth density file format
double InvScaleUn
Define the shrink factor between the box edges.
Definition: ElPoly.h:457
double NanoDist2(double *Pos, int n)
Pointer to a generic function.
Definition: VarData.h:764
int Type
Type of the chain (see list CHAIN_)
Definition: VarData.h:246
bool OpenRisk(char *InFile, int BF)
Opens a file without reallocationg.
Definition: VarData.cpp:126
double pVel(int p, int d)
Return the velocity.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
void ShapeId(int iShape, char *Shape)
Identifier of the shape.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
int pChain(int p)
Return the chain.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
double Height
Height of the cylinder.
Definition: VarData.h:449
int pType(int p)
Return the type.
int CLat2
lateral coordinate
Definition: VarData.h:1078
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
void DrPartPovRay(int p)
PovRay draw function.
void DrBondPovRay(double *Pos1, double *Pos2, float *Color)
Draw the bonds.
double pNanoPos(int n, int d)
Return back folded nano position.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
void HeaderPovRay()
Print the header for povray.
int pNPCh()
Number of particle per chain.
int pNBlock()
Number of blocks.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
FILE * DrawOutFile
Output file for drawing.
Definition: ElPoly.h:463
int pNNano()
Number of nanoparticles.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
Matrice computes the algebric operations on matrices.
void DrNanoPovRay(int n)
PovRay draw function.
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
int EndIdx
End particle position.
Definition: VarData.h:267
double pPos(int p, int d)
Return back folded position.
int pNChain()
Number of chain.
void Conv2rzd(int NSample)
esport the data in radius depth density file format
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
int CLat1
lateral coordinate
Definition: VarData.h:1076
void Conv2Vmd()
esport the data in vmd file format
double Axis[3]
Rotation axis.
Definition: VarData.h:435
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
void Conv2Tecplot(int NBin, int How)
Prepare a countor plot for tecplot.
Definition: ElPolyOutput.cpp:3
void DrField(int NGrid, double IsoLevel, int nNano, FILE *FWrite)
Draw a scalar field.
void ConvLattice(int NSample, char *FName)
Convert into a square lattice.
int pStep()
Number of steps.
int pNPart()
Number of particle.
void LoadDensFile(double **Plot, int NBin)
Load in the array Plot the density of the system.
Definition: VarDataEl.cpp:150