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