FreeFOAM The Cross-Platform CFD Toolkit
fieldview9Reader.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  Reader module for Fieldview9 to read OpenFOAM mesh and data.
26 
27  Creates new 'fvbin' type executable which needs to be installed in place
28  of bin/fvbin.
29 
30  Implements a reader for combined mesh&results on an unstructured mesh.
31 
32  See Fieldview Release 9 Reference Manual and coding in user/ directory
33  of the Fieldview release.
34 
35 \*---------------------------------------------------------------------------*/
36 
37 #include <finiteVolume/fvCFD.H>
38 #include <OpenFOAM/IOobjectList.H>
40 #include <OpenFOAM/pointFields.H>
42 #include "readerDatabase.H"
43 #include <OpenFOAM/wallPolyPatch.H>
44 #include <OpenFOAM/ListOps.H>
45 #include <OpenFOAM/cellModeller.H>
46 
47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
48 
49 extern "C"
50 {
51 // Define various Fieldview constants and prototypes
52 
53 #include "fv_reader_tags.h"
54 
55 static const int FVHEX = 2;
56 static const int FVPRISM = 3;
57 static const int FVPYRAMID = 4;
58 static const int FVTET = 1;
59 static const int ITYPE = 1;
60 
61 unsigned int fv_encode_elem_header(int elem_type, int wall_info[]);
62 void time_step_get_value(float*, int, int*, float*, int*);
63 void fv_error_msg(const char*, const char*);
64 
65 void reg_single_unstruct_reader
66 (
67  char *,
68  void
69  (
70  char*, int*, int*, int*, int*, int*, int*,
71  int[], int*, char[][80], int[], int*, char[][80], int*
72  ),
73  void
74  (
75  int*, int*, int*, float[], int*, float[], int*
76  )
77 );
78 
79 int create_tet(const int, const int[8], const int[]);
80 int create_pyramid(const int, const int[8], const int[]);
81 int create_prism(const int, const int[8], const int[]);
82 int create_hex(const int, const int[8], const int[]);
83 
84 typedef unsigned char uChar;
85 extern uChar create_bndry_face_buffered
86 (
87  int bndry_type,
88  int num_verts,
89  int verts[],
90  int *normals_flags,
91  int num_grid_nodes
92 );
93 
94 /*
95  * just define empty readers here for ease of linking.
96  * Comment out if you have doubly defined linking error on this symbol
97  */
98 void ftn_register_data_readers()
99 {}
100 
101 /*
102  * just define empty readers here for ease of linking.
103  * Comment out if you have doubly defined linking error on this symbol
104  */
105 void ftn_register_functions()
106 {}
107 
108 /*
109  * just define empty readers here for ease of linking.
110  * Comment out if you have doubly defined linking error on this symbol
111  */
112 void register_functions()
113 {}
114 
115 
116 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
117 
118 
119 //
120 // Storage for all Foam state (mainly database & mesh)
121 //
122 static readerDatabase db_;
123 
124 
125 // Write fv error message.
126 static void errorMsg(const string& msg)
127 {
128  fv_error_msg("Foam Reader", msg.c_str());
129 }
130 
131 
132 // Simple check if directory is valid case directory.
133 static bool validCase(const fileName& rootAndCase)
134 {
135  //if (isDir(rootAndCase/"system") && isDir(rootAndCase/"constant"))
136  if (isDir(rootAndCase/"constant"))
137  {
138  return true;
139  }
140  else
141  {
142  return false;
143  }
144 }
145 
146 
147 // Check whether case has topo changes by looking back from last time
148 // to first directory with polyMesh/cells.
149 static bool hasTopoChange(const instantList& times)
150 {
151  label lastIndex = times.size()-1;
152 
153  const Time& runTime = db_.runTime();
154 
155  // Only set time; do not update mesh.
156  runTime.setTime(times[lastIndex], lastIndex);
157 
158  fileName facesInst(runTime.findInstance(polyMesh::meshSubDir, "faces"));
159 
160  // See if cellInst is constant directory. Note extra .name() is for
161  // parallel cases when runTime.constant is '../constant'
162  if (facesInst != runTime.constant().name())
163  {
164  Info<< "Found cells file in " << facesInst << " so case has "
165  << "topological changes" << endl;
166 
167  return true;
168  }
169  else
170  {
171  Info<< "Found cells file in " << facesInst << " so case has "
172  << "no topological changes" << endl;
173 
174  return false;
175  }
176 }
177 
178 
179 static bool selectTime(const instantList& times, int* iret)
180 {
181  List<float> fvTimes(2*times.size());
182 
183  forAll(times, timeI)
184  {
185  fvTimes[2*timeI] = float(timeI);
186  fvTimes[2*timeI+1] = float(times[timeI].value());
187  }
188 
189  int istep;
190 
191  float time;
192 
193  *iret=0;
194 
195  time_step_get_value(fvTimes.begin(), times.size(), &istep, &time, iret);
196 
197  if (*iret == -5)
198  {
199  errorMsg("Out of memory.");
200 
201  return false;
202  }
203  if (*iret == -15)
204  {
205  // Cancel action.
206  return false;
207  }
208  if (*iret != 0)
209  {
210  errorMsg("Unspecified error.");
211 
212  return false;
213  }
214  Info<< "Selected timeStep:" << istep << " time:" << scalar(time) << endl;
215 
216  // Set time and load mesh.
217  db_.setTime(times[istep], istep);
218 
219  return true;
220 }
221 
222 
223 // Gets (names of) all fields in all timesteps.
224 static void createFieldNames
225 (
226  const Time& runTime,
227  const instantList& Times,
228  const word& setName
229 )
230 {
231  // From foamToFieldView9/getFieldNames.H:
232 
237 
238  if (setName.empty())
239  {
240  forAll(Times, timeI)
241  {
242  const word& timeName = Times[timeI].name();
243 
244  // Add all fields to hashtable
245  IOobjectList objects(runTime, timeName);
246 
247  wordList vsNames(objects.names(volScalarField::typeName));
248 
249  forAll(vsNames, fieldI)
250  {
251  volScalarHash.insert(vsNames[fieldI]);
252  }
253 
254  wordList vvNames(objects.names(volVectorField::typeName));
255 
256  forAll(vvNames, fieldI)
257  {
258  volVectorHash.insert(vvNames[fieldI]);
259  }
260  }
261  }
262  db_.setFieldNames(volScalarHash.toc(), volVectorHash.toc());
263 }
264 
265 
266 // Appends interpolated values of fieldName to vars array.
267 static void storeScalarField
268 (
269  const volPointInterpolation& pInterp,
270  const word& fieldName,
271  float vars[],
272  label& pointI
273 )
274 {
275  label nPoints = db_.mesh().nPoints();
276  label nTotPoints = nPoints + db_.polys().size();
277 
278  // Check if present
279  IOobject ioHeader
280  (
281  fieldName,
282  db_.runTime().timeName(),
283  db_.runTime(),
284  IOobject::MUST_READ,
285  IOobject::NO_WRITE
286  );
287 
288  if (ioHeader.headerOk())
289  {
290  Info<< "Storing " << nTotPoints << " of interpolated " << fieldName
291  << endl;
292 
293  volScalarField field(ioHeader, db_.mesh());
294 
295  pointScalarField psf(pInterp.interpolate(field));
296 
297  forAll(psf, i)
298  {
299  vars[pointI++] = float(psf[i]);
300  }
301 
302  const labelList& polys = db_.polys();
303 
304  forAll(polys, i)
305  {
306  label cellI = polys[i];
307 
308  vars[pointI++] = float(field[cellI]);
309  }
310  }
311  else
312  {
313  Info<< "Storing " << nTotPoints << " of dummy values of " << fieldName
314  << endl;
315 
316  for(label i = 0; i < nPoints; i++)
317  {
318  vars[pointI++] = 0.0;
319  }
320 
321  const labelList& polys = db_.polys();
322 
323  forAll(polys, i)
324  {
325  vars[pointI++] = 0.0;
326  }
327  }
328 }
329 
330 
331 // Appends interpolated values of fieldName to vars array.
332 static void storeVectorField
333 (
334  const volPointInterpolation& pInterp,
335  const word& fieldName,
336  float vars[],
337  label& pointI
338 )
339 {
340  label nPoints = db_.mesh().nPoints();
341 
342  label nTotPoints = nPoints + db_.polys().size();
343 
344  // Check if present
345  IOobject ioHeader
346  (
347  fieldName,
348  db_.runTime().timeName(),
349  db_.runTime(),
350  IOobject::MUST_READ,
351  IOobject::NO_WRITE
352  );
353 
354  if (ioHeader.headerOk())
355  {
356  Info<< "Storing " << nTotPoints << " of interpolated " << fieldName
357  << endl;
358 
359  volVectorField field(ioHeader, db_.mesh());
360 
361  for (direction d = 0; d < vector::nComponents; d++)
362  {
363  tmp<volScalarField> tcomp = field.component(d);
364  const volScalarField& comp = tcomp();
365 
366  pointScalarField psf(pInterp.interpolate(comp));
367 
368  forAll(psf, i)
369  {
370  vars[pointI++] = float(psf[i]);
371  }
372 
373  const labelList& polys = db_.polys();
374 
375  forAll(polys, i)
376  {
377  label cellI = polys[i];
378 
379  vars[pointI++] = float(comp[cellI]);
380  }
381  }
382  }
383  else
384  {
385  Info<< "Storing " << nTotPoints << " of dummy values of " << fieldName
386  << endl;
387 
388  for (direction d = 0; d < vector::nComponents; d++)
389  {
390  for(label i = 0; i < nPoints; i++)
391  {
392  vars[pointI++] = 0.0;
393  }
394 
395  const labelList& polys = db_.polys();
396 
397  forAll(polys, i)
398  {
399  vars[pointI++] = 0.0;
400  }
401  }
402  }
403 }
404 
405 
406 // Returns Fieldview face_type of mesh face faceI.
407 static label getFvType(const polyMesh& mesh, const label faceI)
408 {
409  return mesh.boundaryMesh().whichPatch(faceI) + 1;
410 }
411 
412 
413 // Returns Fieldview face_type of face f.
414 static label getFaceType
415 (
416  const polyMesh& mesh,
417  const labelList& faceLabels,
418  const face& f
419 )
420 {
421  // Search in subset faceLabels of faces for index of face f.
422  const faceList& faces = mesh.faces();
423 
424  forAll(faceLabels, i)
425  {
426  label faceI = faceLabels[i];
427 
428  if (f == faces[faceI])
429  {
430  // Convert patch to Fieldview face_type.
431  return getFvType(mesh, faceI);
432  }
433  }
434 
435  FatalErrorIn("getFaceType")
436  << "Cannot find face " << f << " in mesh face subset " << faceLabels
437  << abort(FatalError);
438 
439  return -1;
440 }
441 
442 
443 // Returns Fieldview face_types for set of faces
444 static labelList getFaceTypes
445 (
446  const polyMesh& mesh,
447  const labelList& cellFaces,
448  const faceList& cellShapeFaces
449 )
450 {
451  labelList faceLabels(cellShapeFaces.size());
452 
453  forAll(cellShapeFaces, i)
454  {
455  faceLabels[i] = getFaceType(mesh, cellFaces, cellShapeFaces[i]);
456  }
457  return faceLabels;
458 }
459 
460 
461 /*
462  * Callback for querying file contents. Taken from user/user_unstruct_combined.f
463  */
464 void user_query_file_function
465 (
466  /* input */
467  char* fname, /* filename */
468  int* lenf, /* length of fName */
469  int* iunit, /* fortran unit to use */
470  int* max_grids, /* maximum number of grids allowed */
471  int* max_face_types,/* maximum number of face types allowed */
472  int* max_vars, /* maximum number of result variables allowed per */
473  /* grid point*/
474 
475  /* output */
476  int* num_grids, /* number of grids that will be read from data file */
477  int num_nodes[], /* (array of node counts for all grids) */
478  int* num_face_types, /* number of special face types */
479  char face_type_names[][80], /* array of face-type names */
480  int wall_flags[], /* array of flags for the "wall" behavior */
481  int* num_vars, /* number of result variables per grid point */
482  char var_names[][80], /* array of variable names */
483  int* iret /* return value */
484 )
485 {
486  fprintf(stderr, "\n** user_query_file_function\n");
487 
488  string rootAndCaseString(fname, *lenf);
489 
490  fileName rootAndCase(rootAndCaseString);
491 
492  word setName("");
493 
494  if (!validCase(rootAndCase))
495  {
496  setName = rootAndCase.name();
497 
498  rootAndCase = rootAndCase.path();
499 
500  word setDir = rootAndCase.name();
501 
502  rootAndCase = rootAndCase.path();
503 
504  word meshDir = rootAndCase.name();
505 
506  rootAndCase = rootAndCase.path();
507  rootAndCase = rootAndCase.path();
508 
509  if
510  (
511  setDir == "sets"
512  && meshDir == polyMesh::typeName
513  && validCase(rootAndCase)
514  )
515  {
516  // Valid set (hopefully - cannot check contents of setName yet).
517  }
518  else
519  {
520  errorMsg
521  (
522  "Could not find system/ and constant/ directory in\n"
523  + rootAndCase
524  + "\nPlease select a Foam case directory."
525  );
526 
527  *iret = 1;
528 
529  return;
530  }
531 
532  }
533 
534  fileName rootDir(rootAndCase.path());
535 
536  fileName caseName(rootAndCase.name());
537 
538  // handle trailing '/'
539  if (caseName.empty())
540  {
541  caseName = rootDir.name();
542  rootDir = rootDir.path();
543  }
544 
545  Info<< "rootDir : " << rootDir << endl
546  << "caseName : " << caseName << endl
547  << "setName : " << setName << endl;
548 
549  //
550  // Get/reuse database and mesh
551  //
552 
553  bool caseChanged = db_.setRunTime(rootDir, caseName, setName);
554 
555 
556  //
557  // Select time
558  //
559 
560  instantList Times = db_.runTime().times();
561 
562  // If topo case set database time and update mesh.
563  if (hasTopoChange(Times))
564  {
565  if (!selectTime(Times, iret))
566  {
567  return;
568  }
569  }
570  else if (caseChanged)
571  {
572  // Load mesh (if case changed) to make sure we have nPoints etc.
573  db_.loadMesh();
574  }
575 
576 
577  //
578  // Set output variables
579  //
580 
581  *num_grids = 1;
582 
583  const fvMesh& mesh = db_.mesh();
584 
585  label nTotPoints = mesh.nPoints() + db_.polys().size();
586 
587  num_nodes[0] = nTotPoints;
588 
589  Info<< "setting num_nodes:" << num_nodes[0] << endl;
590 
591  Info<< "setting num_face_types:" << mesh.boundary().size() << endl;
592 
593  *num_face_types = mesh.boundary().size();
594 
595  if (*num_face_types > *max_face_types)
596  {
597  errorMsg("Too many patches. FieldView limit:" + name(*max_face_types));
598 
599  *iret = 1;
600 
601  return;
602  }
603 
604 
605  forAll(mesh.boundaryMesh(), patchI)
606  {
607  const polyPatch& patch = mesh.boundaryMesh()[patchI];
608 
609  strcpy(face_type_names[patchI], patch.name().c_str());
610 
611  if (isA<wallPolyPatch>(patch))
612  {
613  wall_flags[patchI] = 1;
614  }
615  else
616  {
617  wall_flags[patchI] = 0;
618  }
619  Info<< "Patch " << patch.name() << " is wall:"
620  << wall_flags[patchI] << endl;
621  }
622 
623  //- Find all volFields and add them to database
624  createFieldNames(db_.runTime(), Times, setName);
625 
626  *num_vars = db_.volScalarNames().size() + 3*db_.volVectorNames().size();
627 
628  if (*num_vars > *max_vars)
629  {
630  errorMsg("Too many variables. FieldView limit:" + name(*max_vars));
631 
632  *iret = 1;
633 
634  return;
635  }
636 
637 
638  int nameI = 0;
639 
640  forAll(db_.volScalarNames(), i)
641  {
642  const word& fieldName = db_.volScalarNames()[i];
643 
644  const word& fvName = db_.getFvName(fieldName);
645 
646  strcpy(var_names[nameI++], fvName.c_str());
647  }
648 
649  forAll(db_.volVectorNames(), i)
650  {
651  const word& fieldName = db_.volVectorNames()[i];
652 
653  const word& fvName = db_.getFvName(fieldName);
654 
655  strcpy(var_names[nameI++], (fvName + "x;" + fvName).c_str());
656  strcpy(var_names[nameI++], (fvName + "y").c_str());
657  strcpy(var_names[nameI++], (fvName + "z").c_str());
658  }
659 
660  *iret = 0;
661 }
662 
663 
664 /*
665  * Callback for reading timestep. Taken from user/user_unstruct_combined.f
666  */
667 void user_read_one_grid_function
668 (
669  int* iunit, /* in: fortran unit number */
670  int* igrid, /* in: grid number to read */
671  int* nodecnt, /* in: number of nodes to read */
672  float xyz[], /* out: coordinates of nodes: x1..xN y1..yN z1..zN */
673  int* num_vars, /* in: number of results per node */
674  float vars[], /* out: values per node */
675  int* iret /* out: return value */
676 )
677 {
678  fprintf(stderr, "\n** user_read_one_grid_function\n");
679 
680  if (*igrid != 1)
681  {
682  errorMsg("Illegal grid number " + Foam::name(*igrid));
683 
684  *iret = 1;
685 
686  return;
687  }
688 
689  // Get current time
690  instantList Times = db_.runTime().times();
691 
692  // Set database time and update mesh.
693  // Note: this should not be nessecary here. We already have the correct
694  // time set and mesh loaded. This is only nessecary because Fieldview
695  // otherwise thinks the case is non-transient.
696  if (!selectTime(Times, iret))
697  {
698  return;
699  }
700 
701 
702  const fvMesh& mesh = db_.mesh();
703 
704  // With mesh now loaded check for change in number of points.
705  label nTotPoints = mesh.nPoints() + db_.polys().size();
706 
707  if (*nodecnt != nTotPoints)
708  {
709  errorMsg
710  (
711  "nodecnt differs from number of points in mesh.\nnodecnt:"
712  + Foam::name(*nodecnt)
713  + " mesh:"
714  + Foam::name(nTotPoints)
715  );
716 
717  *iret = 1;
718 
719  return;
720  }
721 
722 
723  if
724  (
725  *num_vars
726  != (db_.volScalarNames().size() + 3*db_.volVectorNames().size())
727  )
728  {
729  errorMsg("Illegal number of variables " + name(*num_vars));
730 
731  *iret = 1;
732 
733  return;
734  }
735 
736  //
737  // Set coordinates
738  //
739 
740  const pointField& points = mesh.points();
741 
742  int xIndex = 0;
743  int yIndex = xIndex + nTotPoints;
744  int zIndex = yIndex + nTotPoints;
745 
746  // Add mesh points first.
747  forAll(points, pointI)
748  {
749  xyz[xIndex++] = points[pointI].x();
750  xyz[yIndex++] = points[pointI].y();
751  xyz[zIndex++] = points[pointI].z();
752  }
753 
754  // Add cell centres of polys
755  const pointField& ctrs = mesh.cellCentres();
756 
757  const labelList& polys = db_.polys();
758 
759  forAll(polys, i)
760  {
761  label cellI = polys[i];
762 
763  xyz[xIndex++] = ctrs[cellI].x();
764  xyz[yIndex++] = ctrs[cellI].y();
765  xyz[zIndex++] = ctrs[cellI].z();
766  }
767 
768 
769  //
770  // Define elements by calling fv routines
771  //
772 
773  static const cellModel& tet = *(cellModeller::lookup("tet"));
774  static const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
775  static const cellModel& pyr = *(cellModeller::lookup("pyr"));
776  static const cellModel& prism = *(cellModeller::lookup("prism"));
777  static const cellModel& wedge = *(cellModeller::lookup("wedge"));
778  static const cellModel& hex = *(cellModeller::lookup("hex"));
779  //static const cellModel& splitHex = *(cellModeller::lookup("splitHex"));
780 
781  int tetVerts[4];
782  int pyrVerts[5];
783  int prismVerts[6];
784  int hexVerts[8];
785 
786  int tetFaces[4];
787  int pyrFaces[5];
788  int prismFaces[5];
789  int hexFaces[6];
790 
791  const cellShapeList& cellShapes = mesh.cellShapes();
792  const faceList& faces = mesh.faces();
793  const cellList& cells = mesh.cells();
794  const labelList& owner = mesh.faceOwner();
795  label nPoints = mesh.nPoints();
796 
797  // Fieldview face_types array with all faces marked as internal.
798  labelList internalFaces(6, 0);
799 
800  // Mark all cells next to boundary so we don't have to calculate
801  // wall_types for internal cells and can just pass internalFaces.
802  boolList wallCell(mesh.nCells(), false);
803 
804  label nWallCells = 0;
805 
806  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
807  {
808  label cellI = owner[faceI];
809 
810  if (!wallCell[cellI])
811  {
812  wallCell[cellI] = true;
813 
814  nWallCells++;
815  }
816  }
817 
818  label nPolys = 0;
819 
820  forAll(cellShapes, cellI)
821  {
822  const cellShape& cellShape = cellShapes[cellI];
823  const cellModel& cellModel = cellShape.model();
824  const cell& cellFaces = cells[cellI];
825 
826  int istat = 0;
827 
828  if (cellModel == tet)
829  {
830  tetVerts[0] = cellShape[3] + 1;
831  tetVerts[1] = cellShape[0] + 1;
832  tetVerts[2] = cellShape[1] + 1;
833  tetVerts[3] = cellShape[2] + 1;
834 
835  if (wallCell[cellI])
836  {
837  labelList faceTypes =
838  getFaceTypes(mesh, cellFaces, cellShape.faces());
839 
840  tetFaces[0] = faceTypes[2];
841  tetFaces[1] = faceTypes[3];
842  tetFaces[2] = faceTypes[0];
843  tetFaces[3] = faceTypes[1];
844 
845  istat = create_tet(ITYPE, tetVerts, tetFaces);
846  }
847  else
848  {
849  // All faces internal so use precalculated zero.
850  istat = create_tet(ITYPE, tetVerts, internalFaces.begin());
851  }
852  }
853  else if (cellModel == tetWedge)
854  {
855  prismVerts[0] = cellShape[0] + 1;
856  prismVerts[1] = cellShape[3] + 1;
857  prismVerts[2] = cellShape[4] + 1;
858  prismVerts[3] = cellShape[1] + 1;
859  prismVerts[4] = cellShape[4] + 1;
860  prismVerts[5] = cellShape[2] + 1;
861 
862  if (wallCell[cellI])
863  {
864  labelList faceTypes =
865  getFaceTypes(mesh, cellFaces, cellShape.faces());
866 
867  prismFaces[0] = faceTypes[1];
868  prismFaces[1] = faceTypes[2];
869  prismFaces[2] = faceTypes[3];
870  prismFaces[3] = faceTypes[0];
871  prismFaces[4] = faceTypes[3];
872 
873  istat = create_prism(ITYPE, prismVerts, prismFaces);
874  }
875  else
876  {
877  istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
878  }
879  }
880  else if (cellModel == pyr)
881  {
882  pyrVerts[0] = cellShape[0] + 1;
883  pyrVerts[1] = cellShape[1] + 1;
884  pyrVerts[2] = cellShape[2] + 1;
885  pyrVerts[3] = cellShape[3] + 1;
886  pyrVerts[4] = cellShape[4] + 1;
887 
888  if (wallCell[cellI])
889  {
890  labelList faceTypes =
891  getFaceTypes(mesh, cellFaces, cellShape.faces());
892 
893  pyrFaces[0] = faceTypes[0];
894  pyrFaces[1] = faceTypes[3];
895  pyrFaces[2] = faceTypes[2];
896  pyrFaces[3] = faceTypes[1];
897  pyrFaces[4] = faceTypes[4];
898 
899  istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
900  }
901  else
902  {
903  istat = create_pyramid(ITYPE, pyrVerts, internalFaces.begin());
904  }
905  }
906  else if (cellModel == prism)
907  {
908  prismVerts[0] = cellShape[0] + 1;
909  prismVerts[1] = cellShape[3] + 1;
910  prismVerts[2] = cellShape[4] + 1;
911  prismVerts[3] = cellShape[1] + 1;
912  prismVerts[4] = cellShape[5] + 1;
913  prismVerts[5] = cellShape[2] + 1;
914 
915  if (wallCell[cellI])
916  {
917  labelList faceTypes =
918  getFaceTypes(mesh, cellFaces, cellShape.faces());
919 
920  prismFaces[0] = faceTypes[4];
921  prismFaces[1] = faceTypes[2];
922  prismFaces[2] = faceTypes[3];
923  prismFaces[3] = faceTypes[0];
924  prismFaces[4] = faceTypes[1];
925 
926  istat = create_prism(ITYPE, prismVerts, prismFaces);
927  }
928  else
929  {
930  istat = create_prism(ITYPE, prismVerts, internalFaces.begin());
931  }
932  }
933  else if (cellModel == wedge)
934  {
935  hexVerts[0] = cellShape[0] + 1;
936  hexVerts[1] = cellShape[1] + 1;
937  hexVerts[2] = cellShape[0] + 1;
938  hexVerts[3] = cellShape[2] + 1;
939  hexVerts[4] = cellShape[3] + 1;
940  hexVerts[5] = cellShape[4] + 1;
941  hexVerts[6] = cellShape[6] + 1;
942  hexVerts[7] = cellShape[5] + 1;
943 
944  if (wallCell[cellI])
945  {
946  labelList faceTypes =
947  getFaceTypes(mesh, cellFaces, cellShape.faces());
948 
949  hexFaces[0] = faceTypes[2];
950  hexFaces[1] = faceTypes[3];
951  hexFaces[2] = faceTypes[0];
952  hexFaces[3] = faceTypes[1];
953  hexFaces[4] = faceTypes[4];
954  hexFaces[5] = faceTypes[5];
955 
956  istat = create_hex(ITYPE, hexVerts, hexFaces);
957  }
958  else
959  {
960  istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
961  }
962  }
963  else if (cellModel == hex)
964  {
965  hexVerts[0] = cellShape[0] + 1;
966  hexVerts[1] = cellShape[1] + 1;
967  hexVerts[2] = cellShape[3] + 1;
968  hexVerts[3] = cellShape[2] + 1;
969  hexVerts[4] = cellShape[4] + 1;
970  hexVerts[5] = cellShape[5] + 1;
971  hexVerts[6] = cellShape[7] + 1;
972  hexVerts[7] = cellShape[6] + 1;
973 
974  if (wallCell[cellI])
975  {
976  labelList faceTypes =
977  getFaceTypes(mesh, cellFaces, cellShape.faces());
978 
979  hexFaces[0] = faceTypes[0];
980  hexFaces[1] = faceTypes[1];
981  hexFaces[2] = faceTypes[4];
982  hexFaces[3] = faceTypes[5];
983  hexFaces[4] = faceTypes[2];
984  hexFaces[5] = faceTypes[3];
985 
986  istat = create_hex(ITYPE, hexVerts, hexFaces);
987  }
988  else
989  {
990  istat = create_hex(ITYPE, hexVerts, internalFaces.begin());
991  }
992  }
993  else
994  {
995  forAll(cellFaces, cFaceI)
996  {
997  label faceI = cellFaces[cFaceI];
998 
999  // Get Fieldview facetype (internal/on patch)
1000  label fvFaceType = getFvType(mesh, faceI);
1001 
1002  const face& f = faces[faceI];
1003 
1004  // Indices into storage
1005  label nQuads = 0;
1006  label nTris = 0;
1007 
1008  // Storage for triangles and quads created by face
1009  // decomposition (sized for worst case)
1010  faceList quadFaces((f.size() - 2)/2);
1011  faceList triFaces(f.size() - 2);
1012 
1013  f.trianglesQuads
1014  (
1015  points,
1016  nTris,
1017  nQuads,
1018  triFaces,
1019  quadFaces
1020  );
1021 
1022  // Label of cell centre in fv point list.
1023  label polyCentrePoint = nPoints + nPolys;
1024 
1025  for (label i=0; i<nTris; i++)
1026  {
1027  if (cellI == owner[faceI])
1028  {
1029  tetVerts[0] = triFaces[i][0] + 1;
1030  tetVerts[1] = triFaces[i][1] + 1;
1031  tetVerts[2] = triFaces[i][2] + 1;
1032  tetVerts[3] = polyCentrePoint + 1;
1033  }
1034  else
1035  {
1036  tetVerts[0] = triFaces[i][2] + 1;
1037  tetVerts[1] = triFaces[i][1] + 1;
1038  tetVerts[2] = triFaces[i][0] + 1;
1039  tetVerts[3] = polyCentrePoint + 1;
1040  }
1041 
1042  if (wallCell[cellI])
1043  {
1044  // Outside face is one without polyCentrePoint
1045  tetFaces[0] = fvFaceType;
1046  tetFaces[1] = 0;
1047  tetFaces[2] = 0;
1048  tetFaces[3] = 0;
1049 
1050  istat = create_tet(ITYPE, tetVerts, tetFaces);
1051  }
1052  else
1053  {
1054  istat =
1055  create_tet
1056  (
1057  ITYPE,
1058  tetVerts,
1059  internalFaces.begin()
1060  );
1061  }
1062  }
1063 
1064  for (label i=0; i<nQuads; i++)
1065  {
1066  if (cellI == owner[faceI])
1067  {
1068  pyrVerts[0] = quadFaces[i][3] + 1;
1069  pyrVerts[1] = quadFaces[i][2] + 1;
1070  pyrVerts[2] = quadFaces[i][1] + 1;
1071  pyrVerts[3] = quadFaces[i][0] + 1;
1072  pyrVerts[4] = polyCentrePoint + 1;
1073 
1074  }
1075  else
1076  {
1077  pyrVerts[0] = quadFaces[i][0] + 1;
1078  pyrVerts[1] = quadFaces[i][1] + 1;
1079  pyrVerts[2] = quadFaces[i][2] + 1;
1080  pyrVerts[3] = quadFaces[i][3] + 1;
1081  pyrVerts[4] = polyCentrePoint + 1;
1082  }
1083 
1084  if (wallCell[cellI])
1085  {
1086  // Outside face is one without polyCentrePoint
1087  pyrFaces[0] = fvFaceType;
1088  pyrFaces[1] = 0;
1089  pyrFaces[2] = 0;
1090  pyrFaces[3] = 0;
1091  pyrFaces[4] = 0;
1092 
1093  istat = create_pyramid(ITYPE, pyrVerts, pyrFaces);
1094  }
1095  else
1096  {
1097  istat =
1098  create_pyramid
1099  (
1100  ITYPE,
1101  pyrVerts,
1102  internalFaces.begin()
1103  );
1104  }
1105  }
1106 
1107  if (istat != 0)
1108  {
1109  errorMsg("Error during adding cell " + name(cellI));
1110 
1111  *iret = 1;
1112 
1113  return;
1114  }
1115  }
1116 
1117  nPolys++;
1118  }
1119 
1120  if (istat != 0)
1121  {
1122  errorMsg("Error during adding cell " + name(cellI));
1123 
1124  *iret = 1;
1125 
1126  return;
1127  }
1128  }
1129 
1130 
1131 
1132  //
1133  // Set fieldvalues
1134  //
1135 
1136  pointMesh pMesh(mesh);
1137 
1138  volPointInterpolation pInterp(mesh, pMesh);
1139 
1140  int pointI = 0;
1141 
1142  forAll(db_.volScalarNames(), i)
1143  {
1144  const word& fieldName = db_.volScalarNames()[i];
1145 
1146  storeScalarField(pInterp, fieldName, vars, pointI);
1147  }
1148 
1149 
1150  forAll(db_.volVectorNames(), i)
1151  {
1152  const word& fieldName = db_.volVectorNames()[i];
1153 
1154  storeVectorField(pInterp, fieldName, vars, pointI);
1155  }
1156 
1157  // Return without failure.
1158  *iret = 0;
1159 }
1160 
1161 
1162 void register_data_readers()
1163 {
1164  /*
1165  **
1166  ** You should edit this file to "register" a user-defined data
1167  ** reader with FIELDVIEW, if the user functions being registered
1168  ** here are written in C.
1169  ** You should edit "ftn_register_data_readers.f" if the user functions
1170  ** being registered are written in Fortran.
1171  ** In either case, the user functions being registered may call other
1172  ** functions written in either language (C or Fortran); only the
1173  ** language of the "query" and "read" functions referenced here matters
1174  ** to FIELDVIEW.
1175  **
1176  ** The following shows a sample user-defined data reader being
1177  ** "registered" with FIELDVIEW.
1178  **
1179  ** The "extern void" declarations should match the names of the
1180  ** query and grid-reading functions you are providing. It is
1181  ** strongly suggested that all such names begin with "user" so
1182  ** as not to conflict with global names in FIELDVIEW.
1183  **
1184  ** You may call any combination of the data reader registration
1185  ** functions shown below ("register_two_file_reader" and/or
1186  ** "register_single_file_reader" and/or "register_single_unstruct_reader"
1187  ** and/or "register_double_unstruct_reader") as many times as you like,
1188  ** in order to create several different data readers. Each data reader
1189  ** should of course have different "query" and "read" functions, all of
1190  ** which should also appear in "extern void" declarations.
1191  **
1192  */
1193 
1194  /*
1195  ** like this for combined unstructured grids & results in a single file
1196  */
1197  reg_single_unstruct_reader (
1198  "Foam Reader", /* title you want for data reader */
1199  user_query_file_function, /* whatever you called this */
1200  user_read_one_grid_function /* whatever you called this */
1201  );
1202 }
1203 
1204 
1205 
1206 
1207 }
1208 
1209 
1210 // ************************ vim: set sw=4 sts=4 et: ************************ //