1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #ifndef PVIEW_DATA_GMODEL_H
7 #define PVIEW_DATA_GMODEL_H
8 
9 #include "PViewData.h"
10 #include "GModel.h"
11 #include "SBoundingBox3d.h"
12 
13 template <class Real> class stepData {
14 private:
15   // a pointer to the underlying model
16   GModel *_model;
17   // the unrolled list of all geometrical entities in the model
18   std::vector<GEntity *> _entities;
19   // the bounding box of the view
20   SBoundingBox3d _bbox;
21   // the file the data was read from (if empty, refer to PViewData)
22   std::string _fileName;
23   // the index in the file (if negative, refer to PViewData)
24   int _fileIndex;
25   // the value of the time step and value min/max
26   double _time, _min, _max;
27   // the number of components in the data (one stepData contains only
28   // a single field type)
29   int _numComp;
30   // the values, indexed by MVertex or MElement id numbers (If the
31   // numbering is sparse, or if we only have data for high-id
32   // entities, the vector has zero entries and is thus not
33   // optimal. This is the price to pay if we want 1) rapid access to
34   // the data and 2) not to store any additional info in MVertex or
35   // MElement)
36   //
37   // FIXME: we should change this design and store a vector<int> of tags, and do
38   // indirect addressing, even if it's a bit slower...
39   std::vector<Real *> *_data;
40   // a vector containing the multiplying factor allowing to compute
41   // the number of values stored in _data for each index (number of
42   // values = getMult() * getNumComponents()). If _mult is empty, a
43   // default value of "1" is assumed
44   std::vector<int> _mult;
45   // a vector, indexed by MSH element type, of Gauss point locations
46   // in parametric space
47   std::vector<std::vector<double> > _gaussPoints;
48   // a set of all "partitions" encountered in the data
49   std::set<int> _partitions;
50 
51 public:
52   stepData(GModel *model, int numComp, const std::string &fileName = "",
53            int fileIndex = -1, double time = 0., double min = VAL_INF,
54            double max = -VAL_INF)
_model(model)55     : _model(model), _fileName(fileName), _fileIndex(fileIndex), _time(time),
56       _min(min), _max(max), _numComp(numComp), _data(0)
57   {
58   }
stepData(stepData<Real> & other)59   stepData(stepData<Real> &other) : _data(0)
60   {
61     _model = other._model;
62     _entities = other._entities;
63     _bbox = other._bbox;
64     _fileName = other._fileName;
65     _fileIndex = other._fileIndex;
66     _time = other._time;
67     _min = other._min;
68     _max = other._max;
69     _numComp = other._numComp;
70     if(other._data) {
71       std::size_t n = other.getNumData();
72       _data = new std::vector<Real *>(n, (Real *)0);
73       for(std::size_t i = 0; i < n; i++) {
74         Real *d = other.getData(i);
75         if(d) {
76           int m = other.getMult(i) * _numComp;
77           (*_data)[i] = new Real[m];
78           for(int j = 0; j < m; j++) (*_data)[i][j] = d[j];
79         }
80       }
81     }
82     _mult = other._mult;
83     _gaussPoints = other._gaussPoints;
84     _partitions = other._partitions;
85   }
~stepData()86   ~stepData() { destroyData(); }
fillEntities()87   void fillEntities() { _model->getEntities(_entities); }
computeBoundingBox()88   void computeBoundingBox() { _bbox = _model->bounds(); }
getModel()89   GModel *getModel() { return _model; }
getBoundingBox()90   SBoundingBox3d getBoundingBox() { return _bbox; }
getNumEntities()91   int getNumEntities() { return _entities.size(); }
getEntity(int ent)92   GEntity *getEntity(int ent) { return _entities[ent]; }
getNumComponents()93   int getNumComponents() { return _numComp; }
getMult(int index)94   int getMult(int index)
95   {
96     if(index < 0 || index >= (int)_mult.size()) return 1;
97     return _mult[index];
98   }
getFileName()99   std::string getFileName() { return _fileName; }
setFileName(const std::string & name)100   void setFileName(const std::string &name) { _fileName = name; }
getFileIndex()101   int getFileIndex() { return _fileIndex; }
setFileIndex(int index)102   void setFileIndex(int index) { _fileIndex = index; }
getTime()103   double getTime() { return _time; }
setTime(double time)104   void setTime(double time) { _time = time; }
getMin()105   double getMin() { return _min; }
setMin(double min)106   void setMin(double min) { _min = min; }
getMax()107   double getMax() { return _max; }
setMax(double max)108   void setMax(double max) { _max = max; }
getNumData()109   std::size_t getNumData()
110   {
111     if(!_data) return 0;
112     return _data->size();
113   }
resizeData(int n)114   void resizeData(int n)
115   {
116     if(!_data) _data = new std::vector<Real *>(n, (Real *)0);
117     if(n > (int)_data->size()) _data->resize(n, (Real *)0);
118   }
119   Real *getData(int index, bool allocIfNeeded = false, int mult = 1)
120   {
121     if(index < 0) return 0;
122     if(allocIfNeeded) {
123       if(index >= (int)getNumData()) resizeData(index + 100); // optimize this
124       if(!(*_data)[index]) {
125         (*_data)[index] = new Real[_numComp * mult];
126         for(int i = 0; i < _numComp * mult; i++) (*_data)[index][i] = 0.;
127       }
128       if(mult > 1) {
129         if(index >= (int)_mult.size())
130           _mult.resize(index + 100, 1); // optimize this
131         _mult[index] = mult;
132       }
133     }
134     else {
135       if(index >= (int)getNumData()) return 0;
136     }
137     return (*_data)[index];
138   }
destroyData()139   void destroyData()
140   {
141     if(_data) {
142       for(unsigned int i = 0; i < _data->size(); i++)
143         if((*_data)[i]) delete[](*_data)[i];
144       delete _data;
145       _data = 0;
146     }
147   }
renumberData(const std::map<int,int> & mapping)148   void renumberData(const std::map<int, int> &mapping)
149   {
150     if(!_data) return;
151     int imax = 0, imin = 0;
152     for(auto m : mapping) {
153       imax = std::max(imax, m.second);
154       imin = std::min(imin, m.second);
155     }
156     if(imin < 0) {
157       Msg::Warning("Wrong destination index %d in step data renumbering", imin);
158       return;
159     }
160     std::vector<Real *> data2(imax + 1, nullptr);
161     std::vector<int> mult2(imax + 1, 1);
162     for(auto m : mapping) {
163       if(m.first >= 0 && m.first < (int)_data->size()) {
164         data2[m.second] = (*_data)[m.first];
165       }
166       else {
167         Msg::Warning("Wrong source index %d in step data renumbering", m.first);
168         return;
169       }
170       if(m.first >= 0 && m.first < (int)_mult.size())
171         mult2[m.second] = _mult[m.first];
172     }
173     *_data = data2;
174     _mult = mult2;
175   }
getGaussPoints(int msh)176   std::vector<double> &getGaussPoints(int msh)
177   {
178     if((int)_gaussPoints.size() <= msh) _gaussPoints.resize(msh + 1);
179     return _gaussPoints[msh];
180   }
getPartitions()181   std::set<int> &getPartitions() { return _partitions; }
getMemoryInMb()182   double getMemoryInMb()
183   {
184     double b = 0.;
185     for(std::size_t i = 0; i < getNumData(); i++) b += getMult(i);
186     return b * getNumComponents() * sizeof(Real) / 1024. / 1024.;
187   }
188 };
189 
190 // The data container using elements from one or more GModel(s).
191 class PViewDataGModel : public PViewData {
192 public:
193   enum DataType {
194     NodeData = 1,
195     ElementData = 2,
196     ElementNodeData = 3,
197     GaussPointData = 4,
198     BeamData = 5
199   };
200 
201 private:
202   // the data, indexed by time step
203   std::vector<stepData<double> *> _steps;
204   // the global min/max of the view
205   double _min, _max;
206   // the type of the dataset
207   DataType _type;
208   // cache last element to speed up loops
209   MElement *_getElement(int step, int ent, int ele);
210   MVertex *_getNode(MElement *e, int nod);
211 
212 public:
213   PViewDataGModel(DataType type = NodeData);
214   ~PViewDataGModel();
215   bool finalize(bool computeMinMax = true,
216                 const std::string &interpolationScheme = "");
217   std::string getFileName(int step = -1);
218   int getNumTimeSteps();
219   int getFirstNonEmptyTimeStep(int start = 0);
220   double getTime(int step);
221   double getMin(int step = -1, bool onlyVisible = false, int tensorRep = 0,
222                 int forceNumComponents = 0, int componentMap[9] = nullptr);
223   double getMax(int step = -1, bool onlyVisible = false, int tensorRep = 0,
224                 int forceNumComponents = 0, int componentMap[9] = nullptr);
setMin(double min)225   void setMin(double min) { _min = min; }
setMax(double max)226   void setMax(double max) { _max = max; }
227   SBoundingBox3d getBoundingBox(int step = -1);
setBoundingBox(SBoundingBox3d & box)228   void setBoundingBox(SBoundingBox3d &box) {}
229   int getNumScalars(int step = -1);
230   int getNumVectors(int step = -1);
231   int getNumTensors(int step = -1);
232   int getNumPoints(int step = -1);
233   int getNumLines(int step = -1);
234   int getNumTriangles(int step = -1);
235   int getNumQuadrangles(int step = -1);
236   int getNumPolygons(int step = -1);
237   int getNumTetrahedra(int step = -1);
238   int getNumHexahedra(int step = -1);
239   int getNumPrisms(int step = -1);
240   int getNumPyramids(int step = -1);
241   int getNumTrihedra(int step = -1);
242   int getNumPolyhedra(int step = -1);
243   int getNumEntities(int step = -1);
244   int getNumElements(int step = -1, int ent = -1);
245   int getDimension(int step, int ent, int ele);
246   int getNumNodes(int step, int ent, int ele);
247   int getNode(int step, int ent, int ele, int nod, double &x, double &y,
248               double &z);
249   void setNode(int step, int ent, int ele, int nod, double x, double y,
250                double z);
251   void tagNode(int step, int ent, int ele, int nod, int tag);
252   int getNumComponents(int step, int ent, int ele);
253   int getNumValues(int step, int ent, int ele);
254   void getValue(int step, int ent, int ele, int idx, double &val);
255   void getValue(int step, int ent, int ele, int node, int comp, double &val);
256   void setValue(int step, int ent, int ele, int node, int comp, double val);
257   int getNumEdges(int step, int ent, int ele);
258   int getType(int step, int ent, int ele);
259   void reverseElement(int step, int ent, int ele);
260   void smooth();
261   double getMemoryInMb();
262   bool combineTime(nameData &nd);
263   bool skipEntity(int step, int ent);
264   bool skipElement(int step, int ent, int ele, bool checkVisibility = false,
265                    int samplingRate = 1);
266   bool hasTimeStep(int step);
267   bool hasPartition(int step, int part);
268   bool hasMultipleMeshes();
269   bool hasModel(GModel *model, int step = -1);
isNodeData()270   bool isNodeData() { return _type == NodeData; }
useGaussPoints()271   bool useGaussPoints() { return _type == GaussPointData; }
getModel(int step)272   GModel *getModel(int step) { return _steps[step]->getModel(); }
273   GEntity *getEntity(int step, int ent);
274   MElement *getElement(int step, int entity, int element);
275 
276   // get the data type
getType()277   DataType getType() { return _type; }
278   // direct access to value by index
279   bool getValueByIndex(int step, int dataIndex, int node, int comp,
280                        double &val);
281 
282   // Add some data "on the fly" (data is stored in a map, indexed by
283   // node or element number depending on the type of dataset)
284   bool addData(GModel *model, const std::map<int, std::vector<double> > &data,
285                int step, double time, int partition, int numComp);
286 
287   // Add some data "on the fly", without a map
288   bool addData(GModel *model, const std::vector<std::size_t> &tags,
289                const std::vector<std::vector<double> > &data, int step,
290                double time, int partition, int numComp);
291 
292   // Add homogeneous data "on the fly", without a map
293   bool addData(GModel *model, const std::vector<std::size_t> &tags,
294                const std::vector<double> &data, int step, double time,
295                int partition, int numComp);
296 
297   // Allow to destroy the data
298   void destroyData();
299 
300   // I/O routines
301   bool readMSH(const std::string &viewName, const std::string &fileName,
302                int fileIndex, FILE *fp, bool binary, bool swap, int step,
303                double time, int partition, int numComp, int numNodes,
304                const std::string &interpolationScheme);
305   virtual bool writeMSH(const std::string &fileName, double version = 2.2,
306                         bool binary = false, bool savemesh = true,
307                         bool multipleView = false, int partitionNum = -1,
308                         bool saveInterpolationMatrices = true,
309                         bool forceNodeData = false,
310                         bool forceElementData = false);
311   bool readCGNS(const std::pair<std::string, std::string> &solFieldName,
312                 const std::string &fileName, int index, int fileIndex,
313                 int baseIndex,
314                 const std::vector<std::vector<MVertex *> > &vertPerZone,
315                 const std::vector<std::vector<MElement *> > &eltPerZone);
316   bool readMED(const std::string &fileName, int fileIndex);
317   bool writeMED(const std::string &fileName);
318   bool readPCH(const std::string &fileName, int fileIndex);
319 
320   void importLists(int N[24], std::vector<double> *V[24]);
getStepData(int step)321   stepData<double> *getStepData(int step)
322   {
323     if(step >= 0 && step < (int)_steps.size()) return _steps[step];
324     return nullptr;
325   }
326   void sendToServer(const std::string &name);
327 };
328 
329 #endif
330