1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35
36 Christoph Kloss (DCS Computing GmbH, Linz)
37 Arno Mayrhofer (CFDEMresearch GmbH, Linz)
38
39 Copyright 2016- CFDEMresearch GmbH, Linz
40 Copyright 2012- DCS Computing GmbH, Linz
41 Copyright 2009-2012 JKU Linz
42 ------------------------------------------------------------------------- */
43
44 #include <mpi.h>
45 #include <stdio.h>
46 #include <stdlib.h>
47 #include <string.h>
48 #include "ctype.h"
49 #include "memory.h"
50 #include "input.h"
51 #include "modify.h"
52 #include "update.h"
53 #include "error.h"
54 #include "region.h"
55 #include "domain.h"
56 #include <cmath>
57 #include "vector_liggghts.h"
58 #include "input_mesh_tri.h"
59 #include "tri_mesh.h"
60
61 using namespace LAMMPS_NS;
62
InputMeshTri(LAMMPS * lmp,int argc,char ** argv)63 InputMeshTri::InputMeshTri(LAMMPS *lmp, int argc, char **argv) : Input(lmp, argc, argv),
64 verbose_(false),
65 i_exclusion_list_(0),
66 size_exclusion_list_(0),
67 exclusion_list_(0)
68 {}
69
~InputMeshTri()70 InputMeshTri::~InputMeshTri()
71 {}
72
73 /* ----------------------------------------------------------------------
74 process all input from filename
75 ------------------------------------------------------------------------- */
76
meshtrifile(const char * filename,class TriMesh * mesh,bool verbose,const int size_exclusion_list,int * exclusion_list,class Region * region)77 void InputMeshTri::meshtrifile(const char *filename, class TriMesh *mesh,bool verbose,
78 const int size_exclusion_list, int *exclusion_list,
79 class Region *region)
80 {
81 verbose_ = verbose;
82 size_exclusion_list_ = size_exclusion_list;
83 exclusion_list_ = exclusion_list;
84
85 if(strlen(filename) < 5)
86 error->all(FLERR,"Illegal command, file name too short for input of triangular mesh");
87 const char *ext = &(filename[strlen(filename)-3]);
88
89 // read file
90 // case STL file or VTK file
91
92 bool is_stl = (strcmp(ext,"stl") == 0) || (strcmp(ext,"STL") == 0);
93 bool is_vtk = (strcmp(ext,"vtk") == 0) || (strcmp(ext,"VTK") == 0);
94
95 // error if another nested file still open
96 // if single open file is not stdin, close it
97 // open new filename and set stl___file
98
99 if (me == 0)
100 {
101 nonlammps_file = fopen(filename,"r");
102 if (nonlammps_file == NULL) {
103 char str[512];
104 sprintf(str,"Cannot open mesh file %s",filename);
105 error->one(FLERR,str);
106 }
107 } else nonlammps_file = NULL;
108
109 if(is_stl)
110 {
111 if (comm->me == 0) fprintf(screen,"\nReading STL file '%s' (mesh processing step 1/3) \n",filename);
112 meshtrifile_stl(mesh,region,filename);
113
114 }
115 else if(is_vtk)
116 {
117 if (comm->me == 0) fprintf(screen,"\nReading VTK file '%s' (mesh processing step 1/3) \n",filename);
118 meshtrifile_vtk(mesh,region);
119 }
120 else error->all(FLERR,"Illegal command, need either an STL file or a VTK file as input for triangular mesh.");
121
122 if(nonlammps_file) fclose(nonlammps_file);
123 }
124
125 /* ----------------------------------------------------------------------
126 process VTK file
127 ------------------------------------------------------------------------- */
128
meshtrifile_vtk(class TriMesh * mesh,class Region * region)129 void InputMeshTri::meshtrifile_vtk(class TriMesh *mesh,class Region *region)
130 {
131 int n,m;
132
133 double **points = NULL;
134 int ipoint = 0,npoints = 0;
135
136 int **cells = NULL, *lines = NULL;
137 int icell = 0,ncells = 0;
138
139 int ntris = 0;
140 int iLine = 0;
141
142 int nLines = 0;
143
144 while (1)
145 {
146 // read a line from input script
147 // n = length of line including str terminator, 0 if end of file
148 // if line ends in continuation char '&', concatenate next line
149
150 if (me == 0) {
151 m = 0;
152 while (1) {
153 if (maxline-m < 2) reallocate(line,maxline,0);
154 if (fgets(&line[m],maxline-m,nonlammps_file) == NULL) {
155 if (m) n = strlen(line) + 1;
156 else n = 0;
157 break;
158 }
159 m = strlen(line);
160 if (line[m-1] != '\n') continue;
161
162 m--;
163 while (m >= 0 && isspace(line[m])) m--;
164 if (m < 0 || line[m] != '&') {
165 line[m+1] = '\0';
166 n = m+2;
167 break;
168 }
169 }
170 }
171
172 // bcast the line
173 // if n = 0, end-of-file
174 // error if label_active is set, since label wasn't encountered
175 // if original input file, code is done
176 // else go back to previous input file
177
178 MPI_Bcast(&n,1,MPI_INT,0,world);
179 if (n == 0) {
180 break;
181 }
182
183 if (n > maxline) reallocate(line,maxline,n);
184 MPI_Bcast(line,n,MPI_CHAR,0,world);
185
186 // lines start with 1 (not 0)
187 nLines++;
188
189 if (0 == (nLines % 100000) && comm->me == 0)
190 fprintf(screen," successfully read a chunk of 100000 elements in VTK file\n");
191
192 // parse one line from the file
193 parse_nonlammps();
194 // skip empty lines
195 if(narg == 0){
196 if (me == 0 && verbose_)
197 fprintf(screen,"Note: Skipping empty line in VTK mesh file\n");
198 continue;
199 }
200
201 //increase line counter
202 iLine++;
203
204 if(iLine < 3) continue;
205
206 if(iLine == 3)
207 {
208 if(strcmp(arg[0],"ASCII"))
209 error->all(FLERR,"Expecting ASCII VTK mesh file, cannot continue");
210 continue;
211 }
212
213 if(iLine == 4)
214 {
215 if(strcmp(arg[0],"DATASET UNSTRUCTURED_GRID"))
216 error->all(FLERR,"Expecting ASCII VTK unstructured grid mesh file, cannot continue");
217 continue;
218 }
219
220 if(iLine == 5)
221 {
222 if(strcmp(arg[0],"POINTS"))
223 error->all(FLERR,"Expecting 'POINTS' section in ASCII VTK mesh file, cannot continue");
224 npoints = atoi(arg[1]);
225 memory->create<double>(points,npoints,3,"input_mesh:points");
226 continue;
227 }
228
229 if(iLine <= 5+npoints)
230 {
231 if(narg != 3)
232 error->all(FLERR,"Expecting 3 values for each point in 'POINTS' section of ASCII VTK mesh file, cannot continue");
233
234 points[ipoint][0] = atof(arg[0]);
235 points[ipoint][1] = atof(arg[1]);
236 points[ipoint][2] = atof(arg[2]);
237 ipoint++;
238 continue;
239 }
240
241 if(iLine == 6+npoints)
242 {
243 if(strcmp(arg[0],"CELLS"))
244 error->all(FLERR,"Expecting 'CELLS' section in ASCII VTK mesh file, cannot continue");
245 ncells = atoi(arg[1]);
246 memory->create<int>(cells,ncells,3,"input_mesh:cells");
247 memory->create<int>(lines,ncells,"input_mesh:lines");
248 continue;
249 }
250
251 //copy data of all which have 3 values - can be tri, polygon etc
252 if(iLine <= 6+npoints+ncells)
253 {
254 if(narg == 4)
255 for (int j=0;j<3;j++) cells[icell][j] = atoi(arg[1+j]);
256 else
257 cells[icell][0] = -1;
258
259 lines[icell] = nLines;
260
261 icell++;
262 continue;
263 }
264
265 if(iLine == 7+npoints+ncells)
266 {
267 if(strcmp(arg[0],"CELL_TYPES"))
268 error->all(FLERR,"Expecting 'CELL_TYPES' section in ASCII VTK mesh file, cannot continue");
269 if(ncells != atoi(arg[1]))
270 error->all(FLERR,"Inconsistency in 'CELL_TYPES' section in ASCII VTK mesh file, cannot continue");
271 icell = 0;
272 continue;
273 }
274
275 //only take triangles (cell type 5 according to VTK standard) - count them
276 if(iLine <= 7+npoints+2*ncells)
277 {
278 if(strcmp(arg[0],"5")) cells[icell][0] = -1; //remove if not a tet
279 else ntris++;
280 icell++;
281 continue;
282 }
283
284 }
285
286 //now that everything is parsed, write the data into the mesh
287 for(int i = 0; i < ncells; i++)
288 {
289 if(cells[i][0] == -1) continue;
290 if(size_exclusion_list_ > 0 && lines[i] == exclusion_list_[i_exclusion_list_])
291 {
292 if(i_exclusion_list_ < size_exclusion_list_-1)
293 i_exclusion_list_++;
294 continue;
295 }
296 if(!region || ( region->match(points[cells[i][0]]) && region->match(points[cells[i][1]]) && region->match(points[cells[i][2]]) ) )
297 addTriangle(mesh,points[cells[i][0]],points[cells[i][1]],points[cells[i][2]],lines[i]);
298 }
299
300 memory->destroy<double>(points);
301 memory->destroy<int>(cells);
302 }
303
304 /* ----------------------------------------------------------------------
305 process STL file
306 ------------------------------------------------------------------------- */
307
meshtrifile_stl(class TriMesh * mesh,class Region * region,const char * filename)308 void InputMeshTri::meshtrifile_stl(class TriMesh *mesh,class Region *region, const char *filename)
309 {
310 int n,m;
311 int iVertex = 0;
312 double vertices[3][3];
313 bool insideSolidObject = false;
314 bool insideFacet = false;
315 bool insideOuterLoop = false;
316
317 int nLines = 0, nLinesTri = 0;
318
319 while (1)
320 {
321
322 // read a line from input script
323 // n = length of line including str terminator, 0 if end of file
324 // if line ends in continuation char '&', concatenate next line
325
326 if (me == 0) {
327 m = 0;
328 while (1) {
329 if (maxline-m < 2) reallocate(line,maxline,0);
330 if (fgets(&line[m],maxline-m,nonlammps_file) == NULL) {
331 if (m) n = strlen(line) + 1;
332 else n = 0;
333 break;
334 }
335 m = strlen(line);
336 if (line[m-1] != '\n') continue;
337
338 m--;
339 while (m >= 0 && isspace(line[m])) m--;
340 if (m < 0 || line[m] != '&') {
341 line[m+1] = '\0';
342 n = m+2;
343 break;
344 }
345 }
346 }
347
348 // bcast the line
349 // if n = 0, end-of-file
350 // error if label_active is set, since label wasn't encountered
351 // if original input file, code is done
352 // else go back to previous input file
353
354 MPI_Bcast(&n,1,MPI_INT,0,world);
355 if (n == 0) {
356 break;
357 }
358
359 if (n > maxline) reallocate(line,maxline,n);
360 MPI_Bcast(line,n,MPI_CHAR,0,world);
361
362 // lines start with 1 (not 0)
363 nLines++;
364
365 if (0 == (nLines % 100000) && comm->me == 0)
366 fprintf(screen," successfully read a chunk of 100000 elements in STL file '%s'\n",filename);
367
368 // parse one line from the stl file
369
370 parse_nonlammps();
371
372 // skip empty lines
373 if(narg==0){
374 if (me == 0 && verbose_)
375 fprintf(screen,"Note: Skipping empty line in STL file\n");
376 continue;
377 }
378
379 if (strcmp(arg[0],"solid") != 0 && nLines == 1)
380 {
381 if (me == 0)
382 {
383 fclose(nonlammps_file);
384 if (verbose_)
385 fprintf(screen,"Note: solid keyword not found, assuming binary stl file\n");
386 }
387 nonlammps_file = NULL;
388 meshtrifile_stl_binary(mesh, region, filename);
389 break;
390 }
391
392 // detect begin and end of a solid object, facet and vertices
393 if (strcmp(arg[0],"solid") == 0)
394 {
395 if (insideSolidObject)
396 error->all(FLERR,"Corrupt or unknown STL file: New solid object begins without closing prior solid object.");
397 insideSolidObject=true;
398 if (me == 0 && verbose_){
399 fprintf(screen,"Solid body detected in STL file\n");
400 }
401 }
402 else if (strcmp(arg[0],"endsolid") == 0)
403 {
404 if (!insideSolidObject)
405 error->all(FLERR,"Corrupt or unknown STL file: End of solid object found, but no begin.");
406 insideSolidObject=false;
407 if (me == 0 && verbose_) {
408 fprintf(screen,"End of solid body detected in STL file.\n");
409 }
410 }
411
412 // detect begin and end of a facet within a solids object
413 else if (strcmp(arg[0],"facet") == 0)
414 {
415 if (insideFacet)
416 error->all(FLERR,"Corrupt or unknown STL file: New facet begins without closing prior facet.");
417 if (!insideSolidObject)
418 error->all(FLERR,"Corrupt or unknown STL file: New facet begins outside solid object.");
419 insideFacet = true;
420
421 nLinesTri = nLines;
422
423 // check for keyword normal belonging to facet
424 if (strcmp(arg[1],"normal") != 0)
425 error->all(FLERR,"Corrupt or unknown STL file: Facet normal not defined.");
426
427 // do not import facet normal (is calculated later)
428 }
429 else if (strcmp(arg[0],"endfacet") == 0)
430 {
431 if (!insideFacet)
432 error->all(FLERR,"Corrupt or unknown STL file: End of facet found, but no begin.");
433 insideFacet = false;
434 if (iVertex != 3)
435 error->all(FLERR,"Corrupt or unknown STL file: Number of vertices not equal to three (no triangle).");
436
437 // add triangle to mesh
438 //printVec3D(screen,"vertex",vertices[0]);
439 //printVec3D(screen,"vertex",vertices[1]);
440 //printVec3D(screen,"vertex",vertices[2]);
441
442 // nLinesTri is the line
443
444 if(size_exclusion_list_ > 0 && nLinesTri == exclusion_list_[i_exclusion_list_])
445 {
446
447 if(i_exclusion_list_ < size_exclusion_list_-1)
448 {
449 i_exclusion_list_++;
450
451 while((exclusion_list_[i_exclusion_list_-1] == exclusion_list_[i_exclusion_list_]) && (i_exclusion_list_ < size_exclusion_list_-1))
452 i_exclusion_list_++;
453 }
454 }
455 else if(!region || ( region->match(vertices[0]) && region->match(vertices[1]) && region->match(vertices[2]) ) )
456 {
457 addTriangle(mesh,vertices[0],vertices[1],vertices[2],nLinesTri);
458 }
459
460 if (me == 0) {
461 //fprintf(screen," End of facet detected in in solid body.\n");
462 }
463 }
464
465 //detect begin and end of an outer loop within a facet
466 else if (strcmp(arg[0],"outer") == 0)
467 {
468 if (insideOuterLoop)
469 error->all(FLERR,"Corrupt or unknown STL file: New outer loop begins without closing prior outer loop.");
470 if (!insideFacet)
471 error->all(FLERR,"Corrupt or unknown STL file: New outer loop begins outside facet.");
472 insideOuterLoop = true;
473 iVertex = 0;
474
475 if (me == 0){
476 //fprintf(screen," Outer loop detected in facet.\n");
477 }
478 }
479 else if (strcmp(arg[0],"endloop") == 0)
480 {
481 if (!insideOuterLoop)
482 error->all(FLERR,"Corrupt or unknown STL file: End of outer loop found, but no begin.");
483 insideOuterLoop=false;
484 if (me == 0) {
485 //fprintf(screen," End of outer loop detected in facet.\n");
486 }
487 }
488
489 else if (strcmp(arg[0],"vertex") == 0)
490 {
491 if (!insideOuterLoop)
492 error->all(FLERR,"Corrupt or unknown STL file: Vertex found outside a loop.");
493
494 if (me == 0) {
495 //fprintf(screen," Vertex found.\n");
496 }
497
498 // read the vertex
499 for (int j=0;j<3;j++)
500 vertices[iVertex][j]=atof(arg[1+j]);
501
502 iVertex++;
503 if (iVertex > 3)
504 error->all(FLERR,"Corrupt or unknown STL file: Can not have more than 3 vertices "
505 "in a facet (only triangular meshes supported).");
506 }
507 }
508 }
509
meshtrifile_stl_binary(class TriMesh * mesh,class Region * region,const char * filename)510 void InputMeshTri::meshtrifile_stl_binary(class TriMesh *mesh, class Region *region, const char *filename)
511 {
512 unsigned int num_of_facets;
513 std::ifstream stl_file;
514
515 if (me == 0) {
516 // open file for reading
517 stl_file.open(filename, std::ifstream::in | std::ifstream::binary);
518
519 // read 80 byte header into nirvana
520 for (int i=0; i<20; i++){
521 float dum;
522 stl_file.read((char *)&dum, sizeof(float));
523 }
524
525 // read number of triangles
526 stl_file.read((char *)&num_of_facets, sizeof(int));
527 }
528
529 // communicate number of triangles
530 MPI_Bcast(&num_of_facets,1,MPI_INT,0,world);
531
532 unsigned int count = 0;
533 int nElems = 0;
534 float tri_data[12];
535 while(1) {
536 if (me == 0) {
537 // read one triangle dataset (normal + 3 vertices)
538 for (int i=0; i<12; i++)
539 stl_file.read((char *)&tri_data[i], sizeof(float));
540 // read triangle attribute into nirvana
541 short dum;
542 stl_file.read((char *)&dum, sizeof(short));
543 // error handling
544 if (!stl_file)
545 error->one(FLERR,"Corrupt STL file: Error in reading binary STL file.");
546 }
547 // communicate triangle data
548 MPI_Bcast(&tri_data,12,MPI_FLOAT,0,world);
549 // increase triangle counter
550 count++;
551 if(size_exclusion_list_ > 0 && count == (unsigned int)exclusion_list_[i_exclusion_list_]) {
552 if(i_exclusion_list_ < size_exclusion_list_-1)
553 i_exclusion_list_++;
554 }
555 else
556 {
557 double vert[9];
558 // copy float vertices into double variables (ignore normal at the beginning of tri_data)
559 for (int i=0; i<9; i++)
560 vert[i] = (double) tri_data[i+3];
561 if(!region || ( region->match(&vert[0]) && region->match(&vert[3]) && region->match(&vert[6]) ) )
562 {
563 addTriangle(mesh, &vert[0], &vert[3], &vert[6], count);
564 nElems++;
565 if (0 == (nElems % 100000) && comm->me == 0)
566 fprintf(screen," successfully read a chunk of 100000 elements in STL file '%s'\n",filename);
567 }
568 }
569
570 if (count == num_of_facets)
571 break;
572 }
573
574 if (me == 0)
575 stl_file.close();
576 }
577
578 /* ----------------------------------------------------------------------
579 add a triangle to the mesh
580 ------------------------------------------------------------------------- */
581
addTriangle(TriMesh * mesh,double * a,double * b,double * c,int lineNumber)582 void InputMeshTri::addTriangle(TriMesh *mesh,double *a, double *b, double *c,int lineNumber)
583 {
584 double **nodeTmp = create<double>(nodeTmp,3,3);
585 for(int i=0;i<3;i++){
586 nodeTmp[0][i] = a[i];
587 nodeTmp[1][i] = b[i];
588 nodeTmp[2][i] = c[i];
589 }
590 mesh->addElement(nodeTmp,lineNumber);
591 destroy<double>(nodeTmp);
592 }
593