MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
skeleton.cc
Go to the documentation of this file.
1 
6 #include "globals.hh"
7 #include <iostream>
8 #include <fstream>
9 #include <sstream>
10 #include "voro++.hh"
11 #include "geometry.hh"
12 using namespace std;
13 using namespace voro;
14 using namespace globals;
16 void makeFoamSkeleton(\
17  string filename ,\
18  int ncell ,\
19  int *center_x ,\
20  int *center_y ,\
21  int *center_z ,\
22  bool report )
23 {
24  // Uses Voro++ to create foam skeleton - tessellation on provided seeds.
25  // Skeleton is saved to gnuplot file.
26  int i;
27  double xmin=0,xmax=nx;
28  double ymin=0,ymax=ny;
29  double zmin=0,zmax=nz;
30  const int n_x=6,n_y=6,n_z=6; // Set up the number of blocks that
31  // the container is divided into (voro++)
32  if (report) {
33  cout << "voro++ is creating Voronoi tesselation..." << endl;
34  }
35  // Create a container with the geometry given above, and make it
36  // non-periodic in each of the three coordinates. Allocate space for
37  // eight particles within each computational block
38  container con(\
39  xmin,xmax,ymin,ymax,zmin,zmax,n_x,n_y,n_z,true,true,true,8);
40  for (i=0; i<ncell; i++) { //put cell centers in container
41  con.put(i,center_x[i],center_y[i],center_z[i]);
42  }
43  // Save the Voronoi network of all the particles to text files
44  con.draw_cells_gnuplot(filename.c_str());
45 }
47 void importFoamSkeleton(\
48  string filename ,\
49  double **vert ,\
50  int **vinc ,\
51  int vmax ,\
52  int incmax ,\
53  int &sv ,\
54  bool report )
55 {
56  // Imports foam skeleton from gnuplot file. Creates `vert` and `vinc` and
57  // calculates `sv`.
58  int i,j;
59  double xmin=0,xmax=nx;
60  double ymin=0,ymax=ny;
61  double zmin=0,zmax=nz;
62  double x[6],y[6],z[6]; //vertex coordinates
63  int vind[4]; //indices of incident vertices
64  ifstream fin;
65  string line;
66  bool found;
67  double eps=1e-2; //tolerance for close voronoi vertices
68  if (report) {
69  cout << "loading cell vertices and edges" << endl;
70  }
71  fin.open(filename);
72  if (!fin.is_open()) {
73  cout << "can't open gnuplot file with foam skeleton" << endl;
74  exit(1);
75  }
76  for (i=0; i<vmax; i++) { //initialize vert
77  vert[i][0]=0;
78  vert[i][1]=0;
79  vert[i][2]=0;
80  for (j=0; j<incmax; j++) vinc[i][j]=-1;
81  }
82  sv=0; //number of stored vertices
83  getline(fin,line);
84  istringstream tmp(line);
85  if (!import_vtk) {
86  tmp >> x[0]; tmp >> y[0]; tmp >> z[0];
87  } else {
88  tmp >> z[0]; tmp >> y[0]; tmp >> x[0];
89  x[0]=x[0]*nx;
90  y[0]=y[0]*ny;
91  z[0]=z[0]*nz;
92  }
93  while (getline(fin,line)) {
94  //store coordinates of vertices and their incidence to each other
95  if (!line.empty()) {
96  istringstream tmp(line);
97  if (!import_vtk) {
98  tmp >> x[1]; tmp >> y[1]; tmp >> z[1];
99  } else {
100  tmp >> z[1]; tmp >> y[1]; tmp >> x[1];
101  x[1]=x[1]*nx;
102  y[1]=y[1]*ny;
103  z[1]=z[1]*nz;
104  }
105  if (in_domain(x[0],y[0],z[0]) && in_domain(x[1],y[1],z[1])) {
106  //both vertices are in the domain
107  for (i=0; i<2; i++) {
108  found=false;
109  for (j=0; j<sv; j++) {
110  if (abs(x[i]-vert[j][0])<eps && \
111  abs(y[i]-vert[j][1])<eps && \
112  abs(z[i]-vert[j][2])<eps) {
113  found=true;
114  vind[i]=j; //store the vertex index
115  break;
116  }
117  }
118  if (found == false) {
119  //store the vertex coordinates to the end of matrix
120  vert[sv][0]=x[i];
121  vert[sv][1]=y[i];
122  vert[sv][2]=z[i];
123  vind[i]=sv; //store the vertex index
124  sv++;
125  }
126  }
127  for (i=0; i<incmax; i++) {
128  //vertex 1 is incident to vertex 0
129  if (vinc[vind[0]][i] == vind[1]) break;
130  if (vinc[vind[0]][i] == -1) {
131  vinc[vind[0]][i] = vind[1];
132  break;
133  }
134  }
135  for (i=0; i<incmax; i++) {
136  //vertex 0 is incident to vertex 1
137  if (vinc[vind[1]][i] == vind[0]) break;
138  if (vinc[vind[1]][i] == -1) {
139  vinc[vind[1]][i] = vind[0];
140  break;
141  }
142  }
143  } else if ((in_domain(x[0],y[0],z[0]) && \
144  !in_domain(x[1],y[1],z[1])) || (!in_domain(x[0],y[0],z[0]) \
145  && in_domain(x[1],y[1],z[1]))) {
146  //only one of the vertices is in the domain
147  if (x[0]<0 || x[1]<0) { //mirror the edge over the domain
148  x[2]=x[0]+xmax;
149  x[3]=x[1]+xmax;
150  } else if (x[0]>xmax || x[1]>xmax) {
151  x[2]=x[0]-xmax;
152  x[3]=x[1]-xmax;
153  } else {
154  x[2]=x[0];
155  x[3]=x[1];
156  }
157  if (y[0]<0 || y[1]<0) {
158  y[2]=y[0]+ymax;
159  y[3]=y[1]+ymax;
160  } else if (y[0]>ymax || y[1]>ymax) {
161  y[2]=y[0]-ymax;
162  y[3]=y[1]-ymax;
163  } else {
164  y[2]=y[0];
165  y[3]=y[1];
166  }
167  if (z[0]<0 || z[1]<0) {
168  z[2]=z[0]+zmax;
169  z[3]=z[1]+zmax;
170  } else if (z[0]>zmax || z[1]>zmax) {
171  z[2]=z[0]-zmax;
172  z[3]=z[1]-zmax;
173  } else {
174  z[2]=z[0];
175  z[3]=z[1];
176  }
177  for (i=0; i<4; i++) {
178  found=false;
179  for (j=0; j<sv; j++) {
180  if (abs(x[i]-vert[j][0])<eps && \
181  abs(y[i]-vert[j][1])<eps && \
182  abs(z[i]-vert[j][2])<eps) {
183  found=true;
184  vind[i]=j; //store the vertex index
185  break;
186  }
187  }
188  if (found == false) {
189  //store the vertex coordinates to the end of matrix
190  vert[sv][0]=x[i];
191  vert[sv][1]=y[i];
192  vert[sv][2]=z[i];
193  vind[i]=sv; //store the vertex index
194  sv++;
195  }
196  }
197  for (i=0; i<incmax; i++) {
198  //vertex 1 is incident to vertex 0
199  if (vinc[vind[0]][i] == vind[1]) break;
200  if (vinc[vind[0]][i] == -1) {
201  vinc[vind[0]][i] = vind[1];
202  break;
203  }
204  }
205  for (i=0; i<incmax; i++) {
206  //vertex 0 is incident to vertex 1
207  if (vinc[vind[1]][i] == vind[0]) break;
208  if (vinc[vind[1]][i] == -1) {
209  vinc[vind[1]][i] = vind[0];
210  break;
211  }
212  }
213 
214  for (i=0; i<incmax; i++) {
215  //vertex 3 is incident to vertex 2
216  if (vinc[vind[2]][i] == vind[3]) break;
217  if (vinc[vind[2]][i] == -1) {
218  vinc[vind[2]][i] = vind[3];
219  break;
220  }
221  }
222  for (i=0; i<incmax; i++) {
223  //vertex 2 is incident to vertex 3
224  if (vinc[vind[3]][i] == vind[2]) break;
225  if (vinc[vind[3]][i] == -1) {
226  vinc[vind[3]][i] = vind[2];
227  break;
228  }
229  }
230  } else {
231  if (x[0]<0) {
232  //mirror the edge over the domain for the vertex 0
233  x[2]=x[0]+xmax;
234  x[3]=x[1]+xmax;
235  } else if (x[0]>xmax) {
236  x[2]=x[0]-xmax;
237  x[3]=x[1]-xmax;
238  } else {
239  x[2]=x[0];
240  x[3]=x[1];
241  }
242  if (y[0]<0) {
243  y[2]=y[0]+ymax;
244  y[3]=y[1]+ymax;
245  } else if (y[0]>ymax) {
246  y[2]=y[0]-ymax;
247  y[3]=y[1]-ymax;
248  } else {
249  y[2]=y[0];
250  y[3]=y[1];
251  }
252  if (z[0]<0) {
253  z[2]=z[0]+zmax;
254  z[3]=z[1]+zmax;
255  } else if (z[0]>zmax) {
256  z[2]=z[0]-zmax;
257  z[3]=z[1]-zmax;
258  } else {
259  z[2]=z[0];
260  z[3]=z[1];
261  }
262  if (x[1]<0) { //mirror the edge over the domain
263  x[4]=x[0]+xmax;
264  x[5]=x[1]+xmax;
265  } else if (x[1]>xmax) {
266  x[4]=x[0]-xmax;
267  x[5]=x[1]-xmax;
268  } else {
269  x[4]=x[0];
270  x[5]=x[1];
271  }
272  if (y[1]<0) {
273  y[4]=y[0]+ymax;
274  y[5]=y[1]+ymax;
275  } else if (y[1]>ymax) {
276  y[4]=y[0]-ymax;
277  y[5]=y[1]-ymax;
278  } else {
279  y[4]=y[0];
280  y[5]=y[1];
281  }
282  if (z[1]<0) {
283  z[4]=z[0]+zmax;
284  z[5]=z[1]+zmax;
285  } else if (z[1]>zmax) {
286  z[4]=z[0]-zmax;
287  z[5]=z[1]-zmax;
288  } else {
289  z[4]=z[0];
290  z[5]=z[1];
291  }
292  for (i=2; i<6; i++) {
293  found=false;
294  for (j=0; j<sv; j++) {
295  if (abs(x[i]-vert[j][0])<eps && \
296  abs(y[i]-vert[j][1])<eps && \
297  abs(z[i]-vert[j][2])<eps) {
298  found=true;
299  vind[i-2]=j; //store the vertex index
300  break;
301  }
302  }
303  if (found == false) {
304  //store the vertex coordinates to the end of matrix
305  vert[sv][0]=x[i];
306  vert[sv][1]=y[i];
307  vert[sv][2]=z[i];
308  vind[i-2]=sv; //store the vertex index
309  sv++;
310  }
311  }
312  for (i=0; i<incmax; i++) {
313  //vertex 1 is incident to vertex 0
314  if (vinc[vind[0]][i] == vind[1]) break;
315  if (vinc[vind[0]][i] == -1) {
316  vinc[vind[0]][i] = vind[1];
317  break;
318  }
319  }
320  for (i=0; i<incmax; i++) {
321  //vertex 0 is incident to vertex 1
322  if (vinc[vind[1]][i] == vind[0]) break;
323  if (vinc[vind[1]][i] == -1) {
324  vinc[vind[1]][i] = vind[0];
325  break;
326  }
327  }
328  for (i=0; i<incmax; i++) {
329  //vertex 3 is incident to vertex 2
330  if (vinc[vind[2]][i] == vind[3]) break;
331  if (vinc[vind[2]][i] == -1) {
332  vinc[vind[2]][i] = vind[3];
333  break;
334  }
335  }
336  for (i=0; i<incmax; i++) {
337  //vertex 2 is incident to vertex 3
338  if (vinc[vind[3]][i] == vind[2]) break;
339  if (vinc[vind[3]][i] == -1) {
340  vinc[vind[3]][i] = vind[2];
341  break;
342  }
343  }
344  }
345  x[0]=x[1]; y[0]=y[1]; z[0]=z[1];
346  } else {
347  getline(fin,line); //skip the line,
348  //there are always two empty lines together
349  getline(fin,line); //refresh first vertex
350  istringstream tmp(line);
351  if (!import_vtk) {
352  tmp >> x[0]; tmp >> y[0]; tmp >> z[0];
353  } else {
354  tmp >> z[0]; tmp >> y[0]; tmp >> x[0];
355  x[0]=x[0]*nx;
356  y[0]=y[0]*ny;
357  z[0]=z[0]*nz;
358  }
359  }
360  }
361  fin.close();
362 }
importFoamSkeleton
void importFoamSkeleton(string filename, double **vert, int **vinc, int vmax, int incmax, int &sv, bool report)
Load geometric tessellation, move it to the domain and store it.
Definition: skeleton.cc:47
globals::y
real(dp), dimension(:), allocatable y
state variables
Definition: globals.f90:106
globals::eps
real(dp) eps
foam porosity
Definition: globals.f90:30
globals.hh
Defines global variables, macros, templates and namespace.
in_domain
bool in_domain(double x, double y, double z)
True if the point is in the domain, False otherwise.
Definition: geometry.cc:9
globals::nx
int nx
domain size in X
Definition: globals.cc:19
globals::ncell
integer ncell
number of cells
Definition: globals.f90:15
constants::nz
integer nz
spatial discretization
Definition: constants.f90:32
globals
namespace with global variables
Definition: globals.f90:8
globals::import_vtk
bool import_vtk
import morphology from vtk
Definition: globals.cc:27
makeFoamSkeleton
void makeFoamSkeleton(string filename, int ncell, int *center_x, int *center_y, int *center_z, bool report)
Make geometric tessellation using voro++ and store it to a file.
Definition: skeleton.cc:16
globals::ny
int ny
domain size in Y
Definition: globals.cc:20