MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
spherepack.c
1  /*"
2  License
3  This file is part of Modena.
4 
5  Modena is free software; you can redistribute it and/or modify it under
6  the terms of the GNU General Public License as published by the Free
7  Software Foundation, either version 3 of the License, or (at your option)
8  any later version.
9 
10  Modena is distributed in the hope that it will be useful, but WITHOUT ANY
11  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  FOR A PARTICULAR PURPOSE. See the GNU General Public License
13  for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with Modena. If not, see <http://www.gnu.org/licenses/>.
17 
18  Description
19  Tool for Packing Spheres
20 ****************************************************************
21 
22  Kugelpackprogramm
23 
24  packt Kugeln gleicher und ungleicher Durchmesser durch
25  eine Rearrangementmethode nach Jodrey&Tory
26  und Bargiel&Moscinski
27 
28  programmiert: Antje Elsner
29 
30  (Listenstruktur aus:
31  Moscinski, J., and Bargiel, M., "C-language program
32  for simulation of irregular close packing of hard spheres",
33  Comp Phys Comm, 64 (1991) 183-192
34  http://cpc.cs.qub.ac.uk/summaries/ABZF_v1_0.html)
35 
36  Version: 12.28
37 
38  Datum: 28.08.2008
39 
40 ****************************************************************
41 */
42 
43 #include <stdio.h>
44 #include <math.h>
45 #include <stdlib.h>
46 #include <time.h>
47 #include <string.h>
48 #include <limits.h>
49 
50 #include "spdeclarations.h"
51 
52 //Definitionen
53 #define FANZ 7
54 #define TRUE 1
55 #define FALSE 0
56 
57 #define PI 3.141592653589793238462643383279
58 
59 //Flags
60 int Flag_Random= 0;
61 int Flag_End = FALSE;
62 
63 short Time_out[6];
64 
65 //Dateierweiterungen
66 const char gen_par_psfx[] = ".prj"; //.prj - generelle Parameterdatei
67 const char beg_coord_psfx[]= ".sco"; //.sco - Anfangszustand bei zuf�lliger Erzeugung
68 const char res_coord_psfx[]= ".rco"; //.rco - Ergebniskoordinaten
69 const char inp_coord_psfx[]= ".pco"; //.pco - vorgegebene Koordinatendatei
70 const char res_data_psfx[] = ".rst"; //.rst - Ergebnisdaten
71 
72 //Dateinamen
73 char prj_name[80]; //Projektname
74 char gen_par_name[80]; //.prj - generelle Parameterdatei
75 char beg_coord_name[80]; //.sco - Anfangszustand bei zuf�lliger Erzeugung
76 char res_coord_name[80]; //.rco - Ergebniskoordinaten
77 char inp_coord_name[80]; //.pco - vorgegebene Koordinatendatei
78 char res_data_name[80]; //.rst - Ergebnisdaten
79 
80 //Dateipointer
81 FILE* gen_par_file= NULL; //.prj - generelle Parameterdatei
82 FILE* beg_coord_file= NULL; //.sco - Anfangszustand bei zuf�lliger Erzeugung
83 FILE* res_coord_file= NULL; //.rco - Ergebniskoordinaten
84 FILE* inp_coord_file= NULL; //.pco - vorgegebene Koordinatendatei
85 FILE* res_data_file= NULL; //.rst - Ergebnisdaten
86 
87 //Globale Variablen
88 int No_parts;
89 int No_cells_x, No_cells_y, No_cells_z;
90 int No_cells, Ncell_min;
91 long Max_steps= 100000, Ntau= 1638400;
92 double Epsilon= 0.6, Epsilon_scl;
93 double Din, Dout0, Dout, Dout2;
94 double Pactual, Pnom0= 0.6, Pnomin, Pactual_old;
95 double Diam_dens;
96 double Box_x, Box_y, Box_z, Half_x, Half_y, Half_z;
97 double Relax;
98 double Diam_min, Diam_max;
99 double Dmin, Dmax;
100 
101 //Speicherbereiche
102 int *Link_head, *Link_list;
103 int *Ncell_bound_x, *Ncell_bound_y, *Ncell_bound_z;
104 double *Pbc_x, *Pbc_y, *Pbc_z;
105 double *X, *Y, *Z;
106 double *Diam;
107 double *Force_x, *Force_y, *Force_z;
108 double *k_freq, *l_freq, *g_r;
109 
110 //Parameter f�r Verteilungsfunktionen
111 int Distrfunct= 0;
112 double Alpha = 2.0;
113 double Lambda = 0.04;
114 double Mu = 1.0;
115 double GSigma = 1.0;
116 double Vv = 0.5;
117 double Exp = -0.3;
118 double *Percentage;
119 double *Perc_values;
120 int Perc_number= 0;
121 
122 int seed= 0;
123 
124 long Step;
125 
126 /*
127 ****************************************************************
128 
129  Hauptfunktion
130 
131 ****************************************************************
132 */
133 int main(int argc, char * argv[])
134 {
135  strcpy(prj_name, "Project01");
136 
137  if (argc == 2)
138  {
139  strcpy(prj_name, argv[1]);
140  }
141 
142  make_project();
143 
144  read_param();
145 
146  set_coordinates();
147 
148  list_structure();
149 
150  packing_loop();
151 
152  write_sphere_data();
153 
154  end_run();
155 
156  return 0;
157 }
158 
159 /*
160 ****************************************************************
161 
162  Anlegen und �ffnen der Projektdateien
163 
164 ****************************************************************
165 */
166 void make_project(void)
167 {
168  strcpy(gen_par_name, prj_name); strcat(gen_par_name, gen_par_psfx);
169 
170  strcpy(beg_coord_name, prj_name); strcat(beg_coord_name, beg_coord_psfx);
171  strcpy(res_coord_name, prj_name); strcat(res_coord_name, res_coord_psfx);
172  strcpy(inp_coord_name, prj_name); strcat(inp_coord_name, inp_coord_psfx);
173  strcpy(res_data_name, prj_name); strcat(res_data_name, res_data_psfx);
174 
175  if (( gen_par_file= fopen(gen_par_name, "r")) == NULL)
176  {
177  printf("Keine Steuerparameterdatei %s.prj gefunden.\n", prj_name);
178  fehlerAusgabe("read_param", -3);
179  }
180 
181  if (( beg_coord_file= fopen(beg_coord_name, "w") ) == NULL)
182  {
183  printf("Konnte %s.sco nicht erzeugen.\n", prj_name);
184  fehlerAusgabe("read_param", -3);
185  }
186 
187  if (( inp_coord_file= fopen(inp_coord_name, "r") ) == NULL)
188  {
189  printf("Keine Datendatei %s.pco gefunden.\n", prj_name);
190  // fehlerAusgabe("read_param", -3);
191  }
192 
193  if (( res_data_file= fopen(res_data_name, "w") ) == NULL)
194  {
195  printf("Konnte %s.rst nicht erzeugen.\n", prj_name);
196  fehlerAusgabe("read_param", -3);
197  }
198 }
199 
200 /*
201 ****************************************************************
202 
203  Einlesen der Parameter
204 
205 ****************************************************************
206 */
207 int read_param(void)
208 {
209  int k= 0, m= 0, l= 0;
210  double sum= 0;
211  char line[80];
212 
213  //Zeile ist ein Kommentar
214  fgets(line, 80, gen_par_file);
215  l++;
216 
217  //zweite Zeile
218  fgets(line, 80, gen_par_file);
219  l++;
220 
221  if ((Flag_Random = atoi(line))== EOF)
222  {
223  printf("Fehler in Zeile %d beim Einlesen von Flag_Random.\n", l);
224  fehlerAusgabe("read_param", -1);
225  }
226 
227  //Zeile ist ein Kommentar
228  fgets(line, 80, gen_par_file);
229  l++;
230 
231  //vierte Zeile: Anzahl Kugeln
232  fgets(line, 80, gen_par_file);
233  l++;
234 
235  if (!(No_parts = atoi(line)))
236  {
237  printf("Fehler in Zeile %d.\n", l);
238  fehlerAusgabe("read_param", -1);
239  }
240 
241  //Zeile ist ein Kommentar
242  fgets(line, 80, gen_par_file);
243  l++;
244 
245  //6. 7. 8. Zeile Boxl�ngen
246  fgets(line, 80, gen_par_file);
247  l++;
248 
249  if (!(No_cells_x = atoi(line)))
250  {
251  printf("Fehler in Zeile %d.\n", l);
252  fehlerAusgabe("read_param", -1);
253  }
254 
255  fgets(line, 80, gen_par_file);
256  l++;
257 
258  if (!(No_cells_y = atoi(line)))
259  {
260  printf("Fehler in Zeile %d.\n", l);
261  fehlerAusgabe("read_param", -1);
262  }
263 
264  fgets(line, 80, gen_par_file);
265  l++;
266 
267  if (!(No_cells_z = atoi(line)))
268  {
269  printf("Fehler in Zeile %d.\n", l);
270  fehlerAusgabe("read_param", -1);
271  }
272 
273  //Zeile ist ein Kommentar
274  fgets(line, 80, gen_par_file);
275  l++;
276 
277  //Epsilon
278  fgets(line, 80, gen_par_file);
279  l++;
280 
281  if (!(Epsilon = atof(line)))
282  {
283  printf("Fehler in Zeile %d.\n", l);
284  fehlerAusgabe("read_param", -1);
285  }
286  else
287  {
288  if(Epsilon <= 0)
289  {
290  printf("Parameter Epsilon muss groesser Null sein.\n");
291  fehlerAusgabe("read_param", -1);
292  }
293  }
294 
295  //Zeile ist ein Kommentar
296  fgets(line, 80, gen_par_file);
297  l++;
298 
299  //Ntau
300  fgets(line, 80, gen_par_file);
301  l++;
302 
303  if (!(Ntau = atoi(line)))
304  {
305  printf("Fehler in Zeile %d.\n", l);
306  fehlerAusgabe("read_param", -1);
307  }
308  else
309  {
310  if(Ntau <= 0)
311  {
312  printf("Parameter Ntau muss groesser Null sein.\n");
313  fehlerAusgabe("read_param", -1);
314  }
315  }
316 
317  //Zeile ist ein Kommentar
318  fgets(line, 80, gen_par_file);
319  l++;
320 
321  //nom. Dichte
322  fgets(line, 80, gen_par_file);
323  l++;
324 
325  if (!(Pnom0 = atof(line)))
326  {
327  printf("Fehler in Zeile %d.\n", l);
328  fehlerAusgabe("read_param", -1);
329  }
330  else
331  {
332  if(Pnom0 <= 0 || Pnom0 > 1.0)
333  {
334  printf("Parameter Pnom muss zwischen Null und Eins liegen.\n");
335  fehlerAusgabe("read_param", -1);
336  }
337 
338  }
339 
340  //Zeile ist ein Kommentar
341  fgets(line, 80, gen_par_file);
342  l++;
343 
344  //min. und max. Durchmesser
345  fgets(line, 80, gen_par_file);
346  l++;
347 
348  if (!(Diam_max = atof(line)))
349  {
350  printf("Fehler in Zeile %d.\n", l);
351  fehlerAusgabe("read_param", -1);
352  }
353  else
354  {
355  if(Diam_max <= 0)
356  {
357  printf("Maximaler Durchmesser muss groesser Null sein.\n");
358  fehlerAusgabe("read_param", -1);
359  }
360  }
361 
362  fgets(line, 80, gen_par_file);
363  l++;
364 
365  if (!(Diam_min = atof(line)))
366  {
367  printf("Fehler in Zeile %d.\n", l);
368  fehlerAusgabe("read_param", -1);
369  }
370  else
371  {
372  if(Diam_min <= 0 || Diam_min > Diam_max)
373  {
374  printf("Minimaler Durchmesser muss groesser Null und kleiner als maximaler Durchmesser sein.\n");
375  fehlerAusgabe("read_param", -1);
376  }
377  }
378 
379  //Max. Anzahl Schritte
380  fgets(line, 80, gen_par_file);
381  l++;
382 
383  fgets(line, 80, gen_par_file);
384  l++;
385 
386  if (!(Max_steps = atoi(line)))
387  {
388  printf("Fehler in Zeile %d.\n", l);
389  fehlerAusgabe("read_param", -1);
390  }
391  else
392  {
393  if(Max_steps < 0)
394  {
395  printf("Maximale Anzahl Verdichtungsschritte muss positiv sein.\n");
396  fehlerAusgabe("read_param", -1);
397  }
398  }
399 
400  //Verteilungsfunktion
401  fgets(line, 80, gen_par_file);
402  l++;
403 
404  fgets(line, 80, gen_par_file);
405  l++;
406 
407  if ((Distrfunct = atoi(line))== EOF)
408  {
409  printf("Fehler in Zeile %d.\n", l);
410  fehlerAusgabe("read_param", -1);
411  }
412  else
413  {
414  if(Distrfunct < 0)
415  {
416  printf("Nummer der Verteilungsfunktion muss positiv oder 0 sein.\nEs werden gleich gro�e Kugeln erzeugt.\n");
417  fehlerAusgabe("read_param", -1);
418  }
419  }
420 
421  if (Flag_Random == 1)
422  Distrfunct= 11;
423 
424  //weitere Parameter f�r die Verteilungen
425  fgets(line, 80, gen_par_file);
426  l++;
427 
428  fgets(line, 80, gen_par_file);
429  l++;
430 
431  if (!(Mu = atof(line)))//
432  {
433  printf("Fehler in Zeile %d.\n", l);
434  fehlerAusgabe("read_param", -1);
435  }
436  else
437  {
438  if(Mu < 0)
439  {
440  printf("Mittelwert fuer Normal- oder Log-Normalvert. muss positiv sein.\n");
441  fehlerAusgabe("read_param", -1);
442  }
443  }
444 
445  fgets(line, 80, gen_par_file);
446  l++;
447 
448  if (!(GSigma = atof(line)))//
449  {
450  printf("Fehler in Zeile %d.\n", l);
451  fehlerAusgabe("read_param", -1);
452  }
453  else
454  {
455  if(GSigma < 0)
456  {
457  printf("Standardabweichung muss positiv sein.\n");
458  fehlerAusgabe("read_param", -1);
459  }
460  }
461 
462  fgets(line, 80, gen_par_file);
463  l++;
464 
465  if (!(Vv = atof(line)))//
466  {
467  printf("Fehler in Zeile %d.\n", l);
468  fehlerAusgabe("read_param", -1);
469  }
470  else
471  {
472  if(Vv < 0 || Vv > 1)
473  {
474  printf("Volumenverh�ltnis muss positiv sein und zwischen 0 und 1 liegen.\n");
475  fehlerAusgabe("read_param", -1);
476  }
477  }
478 
479  fgets(line, 80, gen_par_file);
480  l++;
481 
482  if (!(Exp = atof(line)))//
483  {
484  printf("Fehler in Zeile %d.\n", l);
485  fehlerAusgabe("read_param", -1);
486  }
487 
488  //Prozentbins
489  fgets(line, 80, gen_par_file);
490  l++;
491 
492  if (!(Perc_number = atoi(line)))//
493  {
494  printf("Fehler in Zeile %d.\n", l);
495  fehlerAusgabe("read_param", -1);
496  }
497  else
498  {
499  if(Perc_number < 1)
500  {
501  printf("Anzahl Prozentbins muss groesser 1 sein.\n");
502  fehlerAusgabe("read_param", -1);
503  }
504  }
505 
506  Percentage = (double *) calloc(Perc_number+2, sizeof(double));
507  Perc_values = (double *) calloc(Perc_number+1, sizeof(double));
508 
509  if (Distrfunct == 6)
510  {
511  *Percentage= 0;
512 
513  for (k= 1; k<= Perc_number; k++)
514  {
515  fgets(line, 80, gen_par_file);
516  l++;
517 
518  if (!(*(Percentage+k) = atof(line)))//
519  {
520  printf("Fehler in Zeile %d.\n", l);
521  fehlerAusgabe("read_param", -1);
522  }
523 
524  sum+= *(Percentage+k);
525  }
526 
527  if (sum != 1.0)
528  {
529  printf("\nProgrammausfuehrung fehlerhaft beendet.\nBitte geben Sie eine gueltige Prozentverteilung an!\n");
530  exit(1);
531  }
532 
533  for (m= 0; m < Perc_number; m++)
534  {
535  fgets(line, 80, gen_par_file);
536  l++;
537 
538  if (!(*(Perc_values+m) = atof(line)))//
539  {
540  printf("Fehler in Zeile %d.\n", l);
541  fehlerAusgabe("read_param", -1);
542  }
543  }
544  }
545 
546  fclose(gen_par_file);
547  gen_par_file = NULL;
548 
549  return 0;
550 }
551 
552 /*
553 ****************************************************************
554 
555  Belegen der x,y,z-Koordinaten und Durchmesser
556  durch Einlesen oder zuf�lliges Erzeugen
557 
558 ****************************************************************
559 */
560 int set_coordinates(void)
561 {
562  int ipart= 0, cnt= 0, r;
563  double x, y, z, d, max_x= 0, max_y= 0, max_z= 0;
564  char line[128];
565 
566  static long seed = 0;
567  time(&seed);
568  srand(seed);
569 
570  if (Flag_Random != 0)
571  {
572  while(!feof(inp_coord_file))
573  {
574  if(fscanf(inp_coord_file,"%lf %lf %lf %lf", &x, &y, &z, &d) > 1)
575  cnt++;
576  }
577 
578  No_parts = cnt;
579 
580  rewind(inp_coord_file);
581  }
582 
583  Diam = (double *) calloc(No_parts, sizeof(double));
584  X = (double *) calloc(No_parts, sizeof(double));
585  Y = (double *) calloc(No_parts, sizeof(double));
586  Z = (double *) calloc(No_parts, sizeof(double));
587 
588  Force_x = (double *) calloc(No_parts, sizeof(double));
589  Force_y = (double *) calloc(No_parts, sizeof(double));
590  Force_z = (double *) calloc(No_parts, sizeof(double));
591 
592  if (Flag_Random != 0)
593  {
594  r= fscanf(inp_coord_file,"%lf %lf %lf %lf", &x, &y, &z, &d);
595 
596  if(r != 4)
597  {
598  printf("Fehler in Koordinatendatei in Zeile 0.\n");
599  fehlerAusgabe("set_coordinates", -1);
600  }
601 
602  ipart= 0;
603 
604  while(ipart < No_parts)
605  {
606  if(max_x< x) max_x= x;
607  if(max_y< y) max_y= y;
608  if(max_z< z) max_z= z;
609 
610  *(Diam+ipart)= d;
611  *(X+ipart) = x;
612  *(Y+ipart) = y;
613  *(Z+ipart) = z;
614 
615  if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
616  *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
617  *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
618  {
619  fehlerAusgabe("set_coordinates", -1);
620  }
621 
622  r= fscanf(inp_coord_file,"%lf %lf %lf %lf", &x, &y, &z, &d);
623 
624  if(r == 4)
625  ipart++;
626  else
627  break;
628  }
629 
630  if (max_x <= 1.0 && max_y <= 1.0 && max_z <= 1.0)
631  {
632  for (ipart = 0; ipart < No_parts; ipart++)
633  {
634  *(X+ipart) *= No_cells_x;
635  *(Y+ipart) *= No_cells_y;
636  *(Z+ipart) *= No_cells_z;
637 
638  if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
639  *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
640  *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
641  {
642  fehlerAusgabe("set_coordinates", -1);
643  }
644  }
645  }
646  sort();
647 
648  fclose(inp_coord_file);
649  inp_coord_file = NULL;
650  }
651  else
652  {
653  for (ipart = 0; ipart < No_parts; ipart++)
654  {
655  *(X+ipart) = ((double)rand()/(double)RAND_MAX) * (double)No_cells_x;
656  *(Y+ipart) = ((double)rand()/(double)RAND_MAX) * (double)No_cells_y;
657  *(Z+ipart) = ((double)rand()/(double)RAND_MAX) * (double)No_cells_z;
658 
659  if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
660  *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
661  *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
662  {
663  fehlerAusgabe("set_coordinates", -1);
664  }
665  }
666  }
667 
668  return 0;
669 }
670 
671 /*
672 ****************************************************************
673 
674  Erzeugen der Durchmesserliste und Anlegen der Hilfslisten
675  f�r die Verwaltung der Zellzerlegung und
676  der periodischen Randbedingungen
677 
678  -> Verwaltungslisten progr. nach:
679  Moscinski, J., and Bargiel, M., "C-language program
680  for simulation of irregular close packing of hard spheres",
681  Comp Phys Comm, 64 (1991) 183-192
682  http://cpc.cs.qub.ac.uk/summaries/ABZF_v1_0.html
683 
684 ****************************************************************
685 */
686 int list_structure(void)
687 {
688  int ipart, i, imult;
689  int ndim_x, ndim_y, ndim_z;
690  double corr, Summe_Kubikdurchm;
691  double Power = 0.3333333333333333333333333333333;
692  double Sphere_vol = 1.9098593171027440292266051604702;
693 
694  No_cells = No_cells_x * No_cells_y * No_cells_z;
695  Ncell_min = (No_cells_x < No_cells_y) ? No_cells_x : No_cells_y;
696  Ncell_min = (No_cells_z < Ncell_min) ? No_cells_z : Ncell_min;
697  Box_x = (double) No_cells_x;
698  Box_y = (double) No_cells_y;
699  Box_z = (double) No_cells_z;
700  Half_x = 0.5 * Box_x;
701  Half_y = 0.5 * Box_y;
702  Half_z = 0.5 * Box_z;
703  Diam_dens = No_cells * Sphere_vol;
704 
705  Diam_max /= Diam_min;
706  Diam_min = 1.0;
707 
708  Dmin= Diam_min;
709  Dmax= Diam_max;
710 
711  if (Flag_Random == 0)
712  {
713  setDistribution(Distrfunct, Diam);
714  qsort ((void *) Diam, No_parts, sizeof(double), dblcmp);
715  }
716 
717  for (ipart = 0; ipart < No_parts; ipart++)
718  fprintf(beg_coord_file, "%.6f\t%.6f\t%.6f\t%.6f\n", *(X+ipart), *(Y+ipart), *(Z+ipart), *(Diam+ipart));
719 
720  fflush(beg_coord_file);
721 
722  fclose(beg_coord_file);
723  beg_coord_file = NULL;
724 
725  Summe_Kubikdurchm = 0.0;
726 
727  for (ipart = 0; ipart < No_parts; ipart++)
728  Summe_Kubikdurchm += *(Diam+ipart) * *(Diam+ipart) * *(Diam+ipart);
729 
730  Diam_dens /= Summe_Kubikdurchm;
731 
732  Pnomin = Pnom0;
733  Dout0 = pow (Diam_dens * Pnom0 , Power);
734  Dout = Dout0;
735  Relax = 0.5 * Dout / Ntau;
736  Dout2 = Dout * Dout;
737  Epsilon_scl = Epsilon / Dout0;
738 
739  Link_head = (int *) calloc(No_cells, sizeof(int));
740  Link_list = (int *) calloc(No_parts, sizeof(int));
741 
742  ndim_x = 3 * No_cells_x - 2;
743  ndim_y = 3 * No_cells_y - 2;
744  ndim_z = 3 * No_cells_z - 2;
745 
746  Ncell_bound_x = (int *) calloc(ndim_x+1, sizeof(int));
747  if (Ncell_bound_x == NULL)
748  {
749  fehlerAusgabe("list_structure", -2);
750  }
751 
752  Ncell_bound_y = (int *) calloc(ndim_y+1, sizeof(int));
753  if (Ncell_bound_y == NULL)
754  {
755  fehlerAusgabe("list_structure", -2);
756  }
757 
758  Ncell_bound_z = (int *) calloc(ndim_z+1, sizeof(int));
759  if (Ncell_bound_z == NULL)
760  {
761  fehlerAusgabe("list_structure", -2);
762  }
763 
764  Pbc_x = (double *) calloc(ndim_x+1, sizeof(double));
765  if (Pbc_x == NULL)
766  {
767  fehlerAusgabe("list_structure", -2);
768  }
769 
770  Pbc_y = (double *) calloc(ndim_y+1, sizeof(double));
771  if (Pbc_y == NULL)
772  {
773  fehlerAusgabe("list_structure", -2);
774  }
775 
776  Pbc_z = (double *) calloc(ndim_z+1, sizeof(double));
777  if (Pbc_z == NULL)
778  {
779  fehlerAusgabe("list_structure", -2);
780  }
781 
782  imult = 1;
783  corr = -No_cells_x;
784 
785  for (i= 1; i <= ndim_x; i++)
786  {
787  *(Ncell_bound_x+i) = imult * No_cells_y * No_cells_z;
788  *(Pbc_x+i) = corr;
789 
790  ++imult;
791 
792  if (imult >= No_cells_x)
793  {
794  imult = 0;
795  corr += No_cells_x;
796  }
797  }
798 
799  imult = 1;
800  corr = -No_cells_y;
801 
802  for (i = 1; i <= ndim_y; i++)
803  {
804  *(Ncell_bound_y+i) = imult * No_cells_z;
805  *(Pbc_y+i) = corr;
806 
807  ++imult;
808 
809  if (imult >= No_cells_y)
810  {
811  imult = 0;
812  corr += No_cells_y;
813  }
814  }
815 
816  imult = 1;
817  corr = -No_cells_z;
818 
819  for (i= 1; i <= ndim_z; i++)
820  {
821  *(Ncell_bound_z+i) = imult;
822  *(Pbc_z+i) = corr;
823 
824  ++imult;
825 
826  if (imult >= No_cells_z)
827  {
828  imult = 0;
829  corr += No_cells_z;
830  }
831  }
832 
833  return 0;
834 }
835 
836 /*
837 ****************************************************************
838 
839  Hauptschleife
840 
841 ****************************************************************
842 */
843 int packing_loop(void)
844 {
845  date_time (res_data_file, (long *) Time_out);
846 
847  if (Max_steps > 0)
848  {
849  interaction();
850 
851  Pactual = Din*Din*Din / Diam_dens;
852 
853  for (Step= 1; Step <= Max_steps && Flag_End == FALSE; Step++)
854  {
855  Dout -= Relax;
856 
857  Epsilon_scl = Epsilon / Dout;
858 
859  Dout2 = Dout * Dout;
860 
861  motion();
862 
863  interaction();
864 
865  stepon();
866  }
867  result(Step-1);
868  }
869 
870  return 0;
871 }
872 
873 /*
874 ****************************************************************
875 
876  Errechnen der aktuellen Dichten, Erneuern der Relaxationsrate
877 
878 ****************************************************************
879 */
880 int stepon(void)
881 {
882  int iprec, Nsf= 10;
883 
884  Pnomin = Dout * Dout * Dout / Diam_dens;
885  Pactual = Din * Din * Din / Diam_dens;
886 
887  if (Dout <= Din)
888  {
889  double dsave = Dout;
890 
891  Dout = 1.1 * Din;
892 
893  Dout2 = Dout * Dout;
894 
895  interaction();
896 
897  Dout = dsave;
898 
899  Dout2 = Dout * Dout;
900 
901  if (Din * *(Diam+No_parts-1) > Ncell_min)
902  Din = Ncell_min / *(Diam+No_parts - 1);
903 
904  Pactual = Din * Din * Din / Diam_dens;
905 
906  Flag_End= TRUE;
907 
908  return 0;
909  }
910 
911  iprec= (int) -log10 (Pnomin - Pactual);
912 
913  if (iprec >= Nsf)
914  {
915  Relax *= 0.5;
916 
917  Nsf++;
918  }
919  return 0;
920 }
921 
922 /*
923 ****************************************************************
924 
925  Interaktionen zwischen den Kugeln
926 
927 ****************************************************************
928 */
929 int interaction(void)
930 {
931  int ipart, icell;
932  double Rcut;
933 
934  Din = Dout * Dout;
935 
936  for(icell= 0; icell < No_cells; icell++)
937  {
938  *(Link_head+icell)= -1;
939  }
940 
941  for(ipart= 0; ipart < No_parts; ipart++)
942  {
943  *(Force_x+ipart)= 0.0;
944  *(Force_y+ipart)= 0.0;
945  *(Force_z+ipart)= 0.0;
946  }
947 
948  for (ipart= 0; ipart < No_parts; ipart++)
949  {
950  Rcut= ((ipart == 0) ? 0.0 : 0.5) * Dout * (*(Diam + ipart - 1) + *(Diam + ipart));
951 
952  if ((int) (2.0 * Rcut) < Ncell_min - 1)
953  force_part(ipart, Rcut);
954  else
955  force_all(ipart);
956  }
957 
958  Din = sqrt(Din);
959 
960  return 0;
961 }
962 
963 /*
964 ****************************************************************
965 
966  Neuanordnen der Zentren
967 
968 ****************************************************************
969 */
970 void motion(void)
971 {
972  int i;
973  double x_buff, y_buff, z_buff;
974 
975  for (i = 0; i < No_parts; i++)
976  {
977  x_buff = *(X+i)+(*(Force_x+i) * Epsilon_scl) / *(Diam+i);
978  y_buff = *(Y+i)+(*(Force_y+i) * Epsilon_scl) / *(Diam+i);
979  z_buff = *(Z+i)+(*(Force_z+i) * Epsilon_scl) / *(Diam+i);
980 
981  if (x_buff < 0.0)
982  x_buff+= (ceil(-(x_buff/Box_x))*Box_x);
983 
984  if (x_buff >= Box_x)
985  x_buff-= (floor(x_buff/Box_x)*Box_x);
986 
987  if (y_buff < 0.0)
988  y_buff+= (ceil(-(y_buff/Box_y))*Box_y);
989 
990  if (y_buff >= Box_y)
991  y_buff-= (floor(y_buff/Box_y)*Box_y);
992 
993  if (z_buff < 0.0)
994  z_buff+= (ceil(-(z_buff/Box_z))*Box_z);
995 
996  if (z_buff >= Box_z)
997  z_buff-= (floor(z_buff/Box_z)*Box_z);
998 
999  *(X+i) = x_buff;
1000  *(Y+i) = y_buff;
1001  *(Z+i) = z_buff;
1002  }
1003 }
1004 
1005 /*
1006 ****************************************************************
1007 
1008  Verschiebungen bei Nutzung der Nachbarschaftslisten
1009 
1010 ****************************************************************
1011 */
1012 void force_part(int ipart_p, double Rcut)
1013 {
1014  int jpart= 0, ifirst= 0;
1015  int leap_x= 0, leap_y= 0, leap_z= 0;
1016  int icell_x= 0, icell_xy= 0, icell_y= 0, icell_z= 0, icell= 0;
1017  int low_cell_x, low_cell_y, low_cell_z;
1018  int lim_cell_x, lim_cell_y, lim_cell_z;
1019  int idif_x= 0, idif_y= 0, idif_z= 0, idif= 0;
1020 
1021  double pos_x, pos_y, pos_z, dist;
1022  double diam_i, dsum, sigma;
1023  double dif_x, dif_y, dif_z, difpot;
1024  double f_x, f_y, f_z;
1025 
1026  pos_x = *(X+ipart_p);
1027  pos_y = *(Y+ipart_p);
1028  pos_z = *(Z+ipart_p);
1029  diam_i = *(Diam+ipart_p);
1030 
1031  low_cell_x = (int) (pos_x + No_cells_x - Rcut);
1032  low_cell_y = (int) (pos_y + No_cells_y - Rcut);
1033  low_cell_z = (int) (pos_z + No_cells_z - Rcut);
1034  lim_cell_x = (int) (pos_x + No_cells_x + Rcut);
1035  lim_cell_y = (int) (pos_y + No_cells_y + Rcut);
1036  lim_cell_z = (int) (pos_z + No_cells_z + Rcut);
1037  idif_x = (*(Ncell_bound_x+lim_cell_x) > *(Ncell_bound_x+low_cell_x)) ? 0 : 4;
1038  idif_y = (*(Ncell_bound_y+lim_cell_y) > *(Ncell_bound_y+low_cell_y)) ? 0 : 2;
1039  idif_z = (*(Ncell_bound_z+lim_cell_z) > *(Ncell_bound_z+low_cell_z)) ? 0 : 1;
1040  idif = idif_x + idif_y + idif_z;
1041 
1042  switch (idif)
1043  {
1044  case 0:
1045  for (leap_x= low_cell_x; leap_x<= lim_cell_x; leap_x++)
1046  {
1047  icell_x = *(Ncell_bound_x+leap_x);
1048 
1049  for (leap_y= low_cell_y; leap_y<= lim_cell_y; leap_y++)
1050  {
1051  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1052 
1053  for (leap_z= low_cell_z; leap_z<= lim_cell_z; leap_z++)
1054  {
1055  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1056  jpart = *(Link_head+icell);
1057 
1058  while (jpart != -1)
1059  {
1060  dif_x = *(X+jpart) - pos_x;
1061  dif_y = *(Y+jpart) - pos_y;
1062  dif_z = *(Z+jpart) - pos_z;
1063 
1064  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1065  dsum = 0.5 * (diam_i + *(Diam+jpart));
1066  sigma = Dout * dsum;
1067 
1068  if (dist < sigma)
1069  {
1070  if (dist < sqrt(Din)*dsum)
1071  Din = (dist*dist)/(dsum*dsum);
1072 
1073  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1074 
1075  f_x = difpot * dif_x;
1076  f_y = difpot * dif_y;
1077  f_z = difpot * dif_z;
1078 
1079  *(Force_x+ipart_p) -= f_x;
1080  *(Force_y+ipart_p) -= f_y;
1081  *(Force_z+ipart_p) -= f_z;
1082  *(Force_x+jpart) += f_x;
1083  *(Force_y+jpart) += f_y;
1084  *(Force_z+jpart) += f_z;
1085  }
1086  jpart = *(Link_list+jpart);
1087  }
1088  }
1089  }
1090  }
1091  break;
1092 
1093  case 1:
1094  for (leap_x=low_cell_x; leap_x<=lim_cell_x; leap_x++)
1095  {
1096  icell_x = *(Ncell_bound_x+leap_x);
1097 
1098  for (leap_y=low_cell_y; leap_y<=lim_cell_y; leap_y++)
1099  {
1100  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1101 
1102  for (leap_z=low_cell_z; leap_z<=lim_cell_z; leap_z++)
1103  {
1104  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1105  jpart = *(Link_head+icell);
1106 
1107  while (jpart != -1)
1108  {
1109  dif_x = *(X+jpart) - pos_x;
1110  dif_y = *(Y+jpart) - pos_y;
1111  dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1112 
1113  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1114  dsum = 0.5 * (diam_i + *(Diam+jpart));
1115  sigma = Dout * dsum;
1116 
1117  if (dist < sigma)
1118  {
1119  if (dist < sqrt(Din)*dsum)
1120  Din = (dist*dist)/(dsum*dsum);
1121 
1122  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1123 
1124  f_x = difpot * dif_x;
1125  f_y = difpot * dif_y;
1126  f_z = difpot * dif_z;
1127 
1128  *(Force_x+ipart_p) -= f_x;
1129  *(Force_y+ipart_p) -= f_y;
1130  *(Force_z+ipart_p) -= f_z;
1131  *(Force_x+jpart) += f_x;
1132  *(Force_y+jpart) += f_y;
1133  *(Force_z+jpart) += f_z;
1134  }
1135  jpart = *(Link_list+jpart);
1136  }
1137  }
1138  }
1139  }
1140  break;
1141 
1142  case 2:
1143  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1144  {
1145  icell_x = *(Ncell_bound_x+leap_x);
1146 
1147  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1148  {
1149  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1150 
1151  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1152  {
1153  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1154  jpart = *(Link_head+icell);
1155 
1156  while (jpart != -1)
1157  {
1158  dif_x = *(X+jpart) - pos_x;
1159  dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1160  dif_z = *(Z+jpart) - pos_z;
1161 
1162  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1163  dsum = 0.5 * (diam_i + *(Diam+jpart));
1164  sigma = Dout * dsum;
1165 
1166  if (dist < sigma)
1167  {
1168  if (dist < sqrt(Din)*dsum)
1169  Din = (dist*dist)/(dsum*dsum);
1170 
1171  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1172 
1173  f_x = difpot * dif_x;
1174  f_y = difpot * dif_y;
1175  f_z = difpot * dif_z;
1176 
1177  *(Force_x+ipart_p) -= f_x;
1178  *(Force_y+ipart_p) -= f_y;
1179  *(Force_z+ipart_p) -= f_z;
1180  *(Force_x+jpart) += f_x;
1181  *(Force_y+jpart) += f_y;
1182  *(Force_z+jpart) += f_z;
1183  }
1184  jpart = *(Link_list+jpart);
1185  }
1186  }
1187  }
1188  }
1189  break;
1190 
1191  case 3:
1192  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1193  {
1194  icell_x = *(Ncell_bound_x+leap_x);
1195 
1196  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1197  {
1198  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1199 
1200  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1201  {
1202  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1203  jpart = *(Link_head+icell);
1204 
1205  while (jpart != -1)
1206  {
1207  dif_x = *(X+jpart) - pos_x;
1208  dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1209  dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1210 
1211  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1212  dsum = 0.5 * (diam_i + *(Diam+jpart));
1213  sigma = Dout * dsum;
1214 
1215  if (dist < sigma)
1216  {
1217  if (dist < sqrt(Din)*dsum)
1218  Din = (dist*dist)/(dsum*dsum);
1219 
1220  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1221 
1222  f_x = difpot * dif_x;
1223  f_y = difpot * dif_y;
1224  f_z = difpot * dif_z;
1225 
1226  *(Force_x+ipart_p) -= f_x;
1227  *(Force_y+ipart_p) -= f_y;
1228  *(Force_z+ipart_p) -= f_z;
1229  *(Force_x+jpart) += f_x;
1230  *(Force_y+jpart) += f_y;
1231  *(Force_z+jpart) += f_z;
1232  }
1233  jpart = *(Link_list+jpart);
1234  }
1235  }
1236  }
1237  }
1238  break;
1239 
1240  case 4:
1241  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1242  {
1243  icell_x = *(Ncell_bound_x+leap_x);
1244 
1245  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1246  {
1247  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1248 
1249  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1250  {
1251  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1252  jpart = *(Link_head+icell);
1253 
1254  while (jpart != -1)
1255  {
1256  dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1257  dif_y = *(Y+jpart) - pos_y;
1258  dif_z = *(Z+jpart) - pos_z;
1259 
1260  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1261  dsum = 0.5 * (diam_i + *(Diam+jpart));
1262  sigma = Dout * dsum;
1263 
1264  if (dist < sigma)
1265  {
1266  if (dist < sqrt(Din)*dsum)
1267  Din = (dist*dist)/(dsum*dsum);
1268 
1269  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1270 
1271  f_x = difpot * dif_x;
1272  f_y = difpot * dif_y;
1273  f_z = difpot * dif_z;
1274 
1275  *(Force_x+ipart_p) -= f_x;
1276  *(Force_y+ipart_p) -= f_y;
1277  *(Force_z+ipart_p) -= f_z;
1278  *(Force_x+jpart) += f_x;
1279  *(Force_y+jpart) += f_y;
1280  *(Force_z+jpart) += f_z;
1281  }
1282  jpart = *(Link_list+jpart);
1283  }
1284  }
1285  }
1286  }
1287  break;
1288 
1289  case 5:
1290  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1291  {
1292  icell_x = *(Ncell_bound_x+leap_x);
1293 
1294  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1295  {
1296  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1297 
1298  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1299  {
1300  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1301  jpart = *(Link_head+icell);
1302 
1303  while (jpart != -1)
1304  {
1305  dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1306  dif_y = *(Y+jpart) - pos_y;
1307  dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1308 
1309  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1310  dsum = 0.5 * (diam_i + *(Diam+jpart));
1311  sigma = Dout * dsum;
1312 
1313  if (dist < sigma)
1314  {
1315  if (dist < sqrt(Din)*dsum)
1316  Din = (dist*dist)/(dsum*dsum);
1317 
1318  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1319 
1320  f_x = difpot * dif_x;
1321  f_y = difpot * dif_y;
1322  f_z = difpot * dif_z;
1323 
1324  *(Force_x+ipart_p) -= f_x;
1325  *(Force_y+ipart_p) -= f_y;
1326  *(Force_z+ipart_p) -= f_z;
1327  *(Force_x+jpart) += f_x;
1328  *(Force_y+jpart) += f_y;
1329  *(Force_z+jpart) += f_z;
1330  }
1331  jpart = *(Link_list+jpart);
1332  }
1333  }
1334  }
1335  }
1336  break;
1337 
1338  case 6:
1339  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1340  {
1341  icell_x = *(Ncell_bound_x+leap_x);
1342 
1343  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1344  {
1345  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1346 
1347  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1348  {
1349  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1350  jpart = *(Link_head+icell);
1351 
1352  while (jpart != -1)
1353  {
1354  dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1355  dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1356  dif_z = *(Z+jpart) - pos_z;
1357 
1358  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1359  dsum = 0.5 * (diam_i + *(Diam+jpart));
1360  sigma = Dout * dsum;
1361 
1362  if (dist < sigma)
1363  {
1364  if (dist < sqrt(Din)*dsum)
1365  Din = (dist*dist)/(dsum*dsum);
1366 
1367  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1368 
1369  f_x = difpot * dif_x;
1370  f_y = difpot * dif_y;
1371  f_z = difpot * dif_z;
1372 
1373  *(Force_x+ipart_p) -= f_x;
1374  *(Force_y+ipart_p) -= f_y;
1375  *(Force_z+ipart_p) -= f_z;
1376  *(Force_x+jpart) += f_x;
1377  *(Force_y+jpart) += f_y;
1378  *(Force_z+jpart) += f_z;
1379  }
1380  jpart= *(Link_list+jpart);
1381  }
1382  }
1383  }
1384  }
1385  break;
1386 
1387  case 7:
1388  for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1389  {
1390  icell_x = *(Ncell_bound_x+leap_x);
1391 
1392  for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1393  {
1394  icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1395 
1396  for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1397  {
1398  icell = *(Ncell_bound_z+leap_z) + icell_xy;
1399  jpart = *(Link_head+icell);
1400 
1401  while (jpart != -1)
1402  {
1403  dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1404  dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1405  dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1406 
1407  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1408  dsum = 0.5 * (diam_i + *(Diam+jpart));
1409  sigma = Dout * dsum;
1410 
1411  if (dist < sigma)
1412  {
1413  if (dist < sqrt(Din)*dsum)
1414  Din = (dist*dist)/(dsum*dsum);
1415 
1416  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1417 
1418  f_x = difpot * dif_x;
1419  f_y = difpot * dif_y;
1420  f_z = difpot * dif_z;
1421 
1422  *(Force_x+ipart_p) -= f_x;
1423  *(Force_y+ipart_p) -= f_y;
1424  *(Force_z+ipart_p) -= f_z;
1425  *(Force_x+jpart) += f_x;
1426  *(Force_y+jpart) += f_y;
1427  *(Force_z+jpart) += f_z;
1428  }
1429  jpart= *(Link_list+jpart);
1430  }
1431  }
1432  }
1433  }
1434  break;
1435  }
1436 
1437  icell_x = (int) pos_x;
1438  icell_y = (int) pos_y;
1439  icell_z = (int) pos_z;
1440  icell = No_cells_z * (No_cells_y * icell_x + icell_y) + icell_z;
1441  ifirst = *(Link_head+icell);
1442  *(Link_list+ipart_p) = ifirst;
1443  *(Link_head+icell) = ipart_p;
1444 }
1445 
1446 /*
1447 ****************************************************************
1448 
1449  Berechnen der Verschiebung f�r alle Zellen
1450 
1451 ****************************************************************
1452 */
1453 void force_all(int ipart_p)
1454 {
1455  int jpart= 0, ifirst= 0;
1456  int icell_x= 0, icell_y= 0, icell_z= 0, icell= 0;
1457 
1458  double pos_x= 0.0, pos_y= 0.0, pos_z= 0.0, dist= 0.0, sigma= 0.0;
1459  double diam_i= 0.0, dsum= 0.0;
1460  double dif_x= 0.0, dif_y= 0.0, dif_z= 0.0, difpot= 0.0;
1461  double f_x= 0, f_y= 0, f_z= 0;
1462 
1463  pos_x = *(X+ipart_p);
1464  pos_y = *(Y+ipart_p);
1465  pos_z = *(Z+ipart_p);
1466  diam_i = *(Diam+ipart_p);
1467 
1468  for (icell= 0; icell < No_cells; icell++)
1469  {
1470  jpart = *(Link_head+icell);
1471 
1472  while (jpart != -1)
1473  {
1474  dif_x = *(X+jpart) - pos_x;
1475  dif_y = *(Y+jpart) - pos_y;
1476  dif_z = *(Z+jpart) - pos_z;
1477 
1478  if (dif_x < -Half_x)
1479  dif_x += Box_x;
1480 
1481  if (dif_y < -Half_y)
1482  dif_y += Box_y;
1483 
1484  if (dif_z < -Half_z)
1485  dif_z += Box_z;
1486 
1487  if (dif_x > Half_x)
1488  dif_x -= Box_x;
1489 
1490  if (dif_y > Half_y)
1491  dif_y -= Box_y;
1492 
1493  if (dif_z > Half_z)
1494  dif_z -= Box_z;
1495 
1496  dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1497  dsum = 0.5 * (diam_i + *(Diam+jpart));
1498  sigma = Dout * dsum;
1499 
1500  if (dist < sigma)
1501  {
1502  if (dist < sqrt(Din)*dsum)
1503  Din = (dist*dist)/(dsum*dsum);
1504 
1505  difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma * sigma))/dist;
1506 
1507  f_x = difpot * dif_x;
1508  f_y = difpot * dif_y;
1509  f_z = difpot * dif_z;
1510 
1511  *(Force_x+ipart_p) -= f_x;
1512  *(Force_y+ipart_p) -= f_y;
1513  *(Force_z+ipart_p) -= f_z;
1514  *(Force_x+jpart) += f_x;
1515  *(Force_y+jpart) += f_y;
1516  *(Force_z+jpart) += f_z;
1517  }
1518  jpart = *(Link_list+jpart);
1519  }
1520  }
1521 
1522  icell_x = (int) pos_x;
1523  icell_y = (int) pos_y;
1524  icell_z = (int) pos_z;
1525  icell = No_cells_z * (No_cells_y * icell_x + icell_y) + icell_z;
1526  ifirst = *(Link_head+icell);
1527  *(Link_list+ipart_p) = ifirst;
1528  *(Link_head+icell) = ipart_p;
1529 }
1530 
1531 
1532 /*
1533 ****************************************************************
1534 
1535  Hilfsfunktion Tauschen
1536 
1537 ****************************************************************
1538 */
1539 void swap(int i, int j)
1540 {
1541  double temp;
1542  temp = *(Diam+j), *(Diam+j) = *(Diam+i), *(Diam+i) = temp;
1543  temp = *(X+j), *(X+j) = *(X+i), *(X+i) = temp;
1544  temp = *(Y+j), *(Y+j) = *(Y+i), *(Y+i) = temp;
1545  temp = *(Z+j), *(Z+j) = *(Z+i), *(Z+i) = temp;
1546 }
1547 
1548 /*
1549 ****************************************************************
1550 
1551  Hilfsfunktion Sortieren
1552 
1553 ****************************************************************
1554 */
1555 void sort(void)
1556 {
1557  int gap, i, j;
1558  for(gap = No_parts/2; gap > 0; gap /= 2)
1559  for(i = gap; i < No_parts; i++)
1560  for(j = i - gap; j >= 0 && *(Diam+j) > *(Diam+j+gap); j -= gap)
1561  swap(j, j+gap);
1562 }
1563 
1564 /*
1565 ****************************************************************
1566 
1567  Ausgabe des Packungsergebnisses
1568 
1569 ****************************************************************
1570 */
1571 void result(int nstep)
1572 {
1573  int noDist;
1574  char* DistFunct[]= {"constant","two-point","uniform","gaussian","log-normal","power","percentage","predefined (file)"};
1575 
1576  if(Distrfunct == 11)
1577  noDist= 7;
1578  else
1579  noDist= Distrfunct;
1580 
1581  fprintf(res_data_file, "----------------------------------------\n");
1582  if(Distrfunct == 11)
1583  fprintf(res_data_file,"\nDistribution : %s", "predefined (file)");
1584  else
1585  fprintf(res_data_file,"\nDistribution : %s", DistFunct[Distrfunct]);
1586 
1587  fprintf(res_data_file,"\nCells in x-dir. : %d", No_cells_x);
1588  fprintf(res_data_file,"\nCells in y-dir. : %d", No_cells_y);
1589  fprintf(res_data_file,"\nCells in z-dir. : %d", No_cells_z);
1590  fprintf(res_data_file,"\nNumber of Spheres: %d", No_parts);
1591  fprintf(res_data_file,"\nDensity : %5.13f",Pactual);
1592  fprintf(res_data_file,"\nCalculation steps: %d\n", nstep);
1593  fprintf(res_data_file, "----------------------------------------\n");
1594 
1595  date_time (res_data_file, (long *) Time_out);
1596 
1597  fflush(res_data_file);
1598 
1599  return;
1600 }
1601 
1602 /*
1603 ****************************************************************
1604 
1605  Freigeben von Speicherbereichen und Schlie�en von Dateien
1606 
1607 ****************************************************************
1608 */
1609 void end_run(void)
1610 {
1611  if (Link_head != NULL)
1612  free (Link_head);
1613 
1614  if (Link_list != NULL)
1615  free (Link_list);
1616 
1617  if (Force_x != NULL)
1618  free (Force_x);
1619  if (Force_y != NULL)
1620  free (Force_y);
1621  if (Force_z != NULL)
1622  free (Force_z);
1623 
1624  if (Ncell_bound_x != NULL)
1625  free (Ncell_bound_x);
1626  if (Ncell_bound_y != NULL)
1627  free (Ncell_bound_y);
1628  if (Ncell_bound_z != NULL)
1629  free (Ncell_bound_z);
1630 
1631  if (Pbc_x != NULL)
1632  free (Pbc_x);
1633  if (Pbc_y != NULL)
1634  free (Pbc_y);
1635  if (Pbc_z != NULL)
1636  free (Pbc_z);
1637 
1638  if (Diam != NULL)
1639  free (Diam);
1640  if (X != NULL)
1641  free (X);
1642  if (Y != NULL)
1643  free (Y);
1644  if (Z != NULL)
1645  free (Z);
1646 
1647  if (Percentage != NULL)
1648  free (Percentage);
1649  if (Perc_values != NULL)
1650  free (Perc_values);
1651 
1652  if (gen_par_file != NULL) fclose(gen_par_file);
1653  if (beg_coord_file != NULL) fclose(beg_coord_file);
1654  if (res_coord_file != NULL) fclose(res_coord_file);
1655  if (inp_coord_file != NULL) fclose(inp_coord_file);
1656  if (res_data_file != NULL) fclose(res_data_file);
1657 }
1658 
1659 /*
1660 ****************************************************************
1661 
1662  Ausgeben der x,y,z-Koordinaten und Durchmesser
1663 
1664 ****************************************************************
1665 */
1666 void write_sphere_data()
1667 {
1668  int count= 1000, i;
1669  char filename[80]= "";
1670 
1671  strcpy(filename, res_coord_name);
1672 
1673  res_coord_file= fopen(filename, "w");
1674 
1675  for (i= 0; i < No_parts; i++)
1676  {
1677  fprintf(res_coord_file, "%.14lf\t%.14lf\t%.14lf\t%.14lf\n", *(X+i), *(Y+i), *(Z+i), *(Diam+i)*Din);
1678  fflush(res_coord_file);
1679  }
1680 
1681  return;
1682 }
1683 
1684 /*
1685 ****************************************************************
1686 
1687  Hilfsfunktion Vergleich
1688 
1689 ****************************************************************
1690 */
1691 int dblcmp(const void *v1, const void *v2)
1692 {
1693  double t = *(double *)v1 - *(double *)v2;
1694 
1695  return ((t > 0) ? 1 : -1);
1696 }
1697 
1698 /*
1699 ****************************************************************
1700 
1701  Hilfsfunktion Zeitausgabe
1702 
1703 ****************************************************************
1704 */
1705 void date_time(FILE * pfile_p, long * time_p)
1706 {
1707  long clock;
1708 
1709  *time_p = time(&clock);
1710  fprintf(pfile_p, "%s", ctime(&clock));
1711 }
1712 
1713 /*
1714 ****************************************************************
1715 
1716  Fehlerbehandlung
1717 
1718 ****************************************************************
1719 */
1720 int fehlerAusgabe(char *funktion, int fehlerNummer)
1721 {
1722  int fehlerIndex= (-fehlerNummer) - 1;
1723  int fehlerart= 1; /* 0: behebbarer Fehler; 1: harter Fehler */
1724 
1725  const char fehlerNachricht[FANZ][60+1]=
1726  {
1727  " -1 Fehler in den Argumenten",
1728  " -2 Speicherplatz konnte nicht bereit gestellt werden",
1729  " -3 Datei konnte nicht zum Lesen geoeffnet werden",
1730  " -4 Datensatz konnte nicht gelesen werden",
1731  " -5 Datei konnte nicht zum Schreiben geoeffnet werden",
1732  " -6 Datensatz konnte nicht geschrieben werden",
1733  " -7 Datei konnte nicht geschlossen werden",
1734  };
1735 
1736  if (fehlerIndex >= 0 && fehlerIndex < FANZ)
1737  {
1738  printf("\nFehlermeldung: %s\n", fehlerNachricht[fehlerIndex]);
1739  }
1740 
1741  if (fehlerart)
1742  {
1743  end_run();
1744 
1745  exit(0);
1746  }
1747  return fehlerart;
1748 }
globals::temp
real(dp) temp
temperature (K)
Definition: globals.f90:49
globals::y
real(dp), dimension(:), allocatable y
state variables
Definition: globals.f90:106
globals::sigma
real(dp) sigma
interfacial tension
Definition: globals.f90:49
globals::time
real(dp) time
time (s)
Definition: globals.f90:49
main
int main(int argc, char *argv[])
Reads parameters. Creates struts and walls. Saves foam morphology to a file.
Definition: foams.cc:23
bubbleGrowth.bubbleGrowth.m
m
Bubble Growth Application Recipe.
Definition: bubbleGrowth.py:59
globals::d
real(dp), dimension(:), allocatable d
diffusion coefficients (for each dissolved gas)
Definition: globals.f90:106