From 1be26a41d0826f1bcb11c4108010226c617a479b Mon Sep 17 00:00:00 2001
From: Michel Rasquin <michel.rasquin@cenaero.be>
Date: Thu, 17 Sep 2015 08:07:58 +0000
Subject: [PATCH] Implementation of a VTK interface for visualization in
parallel with ParaView, with the capacity to adapt the views: - VTK writer
for GMSH GUI - VTK static variables for GMSH external library (linked to
ParaView plugin)
---
CMakeLists.txt | 2 +-
Common/CreateFile.cpp | 2 +
Common/Gmsh.cpp | 14 +
Common/GmshDefines.h | 1 +
Fltk/fileDialogs.cpp | 162 ++++++
Fltk/fileDialogs.h | 1 +
Fltk/graphicWindow.cpp | 3 +
Post/PView.h | 5 +
Post/PViewData.cpp | 28 +
Post/PViewData.h | 14 +
Post/PViewDataGModel.cpp | 4 +
Post/PViewIO.cpp | 12 +-
Post/adaptiveData.cpp | 1040 +++++++++++++++++++++++++++++++++++++-
Post/adaptiveData.h | 268 +++++++++-
14 files changed, 1547 insertions(+), 9 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 90a6d40206..bb148b2025 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -143,7 +143,7 @@ set(GMSH_API
Solver/thermicSolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h Post/OctreePost.h
Post/PViewDataList.h Post/PViewDataGModel.h Post/PViewOptions.h Post/ColorTable.h
- Numeric/nodalBasis.h
+ Numeric/nodalBasis.h Post/adaptiveData.h
Graphics/drawContext.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
contrib/kbipack/gmp_blas.h contrib/kbipack/mpz.h
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index d12c2dd7b5..9eab419c66 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -34,6 +34,7 @@ int GetFileFormatFromExtension(const std::string &ext)
if (ext == ".geo") return FORMAT_GEO;
else if(ext == ".msh") return FORMAT_MSH;
else if(ext == ".pos") return FORMAT_POS;
+ else if(ext == ".pvtu") return FORMAT_PVTU;
else if(ext == ".opt") return FORMAT_OPT;
else if(ext == ".unv") return FORMAT_UNV;
else if(ext == ".vtk") return FORMAT_VTK;
@@ -91,6 +92,7 @@ std::string GetDefaultFileName(int format)
case FORMAT_GEO: name += ".geo_unrolled"; break;
case FORMAT_MSH: name += ".msh"; break;
case FORMAT_POS: name += ".pos"; break;
+ case FORMAT_PVTU: name += ".pvtu"; break;
case FORMAT_OPT: name += ".opt"; break;
case FORMAT_UNV: name += ".unv"; break;
case FORMAT_VTK: name += ".vtk"; break;
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index 3a8c67f289..a3df72026f 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -229,6 +229,20 @@ int GmshWriteFile(const std::string &fileName)
int GmshFinalize()
{
+ //michel.rasquin@cenaero.be:
+ // Cleaning routine critical for ParaView plugin,
+ // where GMSH can be called several time successively.
+ // For that purpose, static variables need to be cleared.
+
+ // Clear all PViews
+ PView::list.clear();
+
+ // Clear all GModels
+ GModel::list.clear();
+
+ // Clear files from CTX::instance()
+ CTX::instance()->files.clear();
+
return 1;
}
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 8a3b64770b..31c0e3780c 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -51,6 +51,7 @@
#define FORMAT_SU2 42
#define FORMAT_MPEG_PREVIEW 43
#define FORMAT_PGF 44
+#define FORMAT_PVTU 45
// Element types
#define TYPE_PNT 1
diff --git a/Fltk/fileDialogs.cpp b/Fltk/fileDialogs.cpp
index ced12c97b6..06a84e4868 100644
--- a/Fltk/fileDialogs.cpp
+++ b/Fltk/fileDialogs.cpp
@@ -1414,6 +1414,168 @@ int posFileDialog(const char *name)
return 0;
}
+static void _saveAdaptedViews(const std::string &name, int useDefaultName, int which, bool isBinary,
+ int adaptLev, double adaptErr, int npart, bool canAppend)
+{
+ if(PView::list.empty()){
+ Msg::Error("No views to save");
+ }
+ else if(which == 0){
+ int iview = FlGui::instance()->options->view.index;
+ if(iview < 0 || iview >= (int)PView::list.size()){
+ Msg::Info("No or invalid current view: saving View[0]");
+ iview = 0;
+ }
+ PView::list[iview]->writeAdapt(name, useDefaultName, isBinary, adaptLev, adaptErr, npart);
+ }
+ else if(which == 1){
+ int numVisible = 0;
+ for(unsigned int i = 0; i < PView::list.size(); i++)
+ if(PView::list[i]->getOptions()->visible)
+ numVisible++;
+ if(!numVisible){
+ Msg::Error("No visible view");
+ }
+ else{
+ bool first = true;
+ for(unsigned int i = 0; i < PView::list.size(); i++){
+ if(PView::list[i]->getOptions()->visible){
+ std::string fileName = name;
+ if(!canAppend && numVisible > 1){
+ std::ostringstream os;
+ os << "_" << i;
+ fileName += os.str();
+ }
+ PView::list[i]->writeAdapt(fileName, useDefaultName, isBinary, adaptLev, adaptErr,
+ npart, first ? false : canAppend);
+ first = false;
+ }
+ }
+ }
+ }
+ else{
+ for(unsigned int i = 0; i < PView::list.size(); i++){
+ std::string fileName = name;
+ if(!canAppend && PView::list.size() > 1){
+ std::ostringstream os;
+ os << "_" << i;
+ fileName += os.str();
+ }
+ PView::list[i]->writeAdapt(fileName, useDefaultName, isBinary, adaptLev, adaptErr,
+ npart, i ? canAppend : false);
+ }
+ }
+}
+
+int pvtuAdaptFileDialog(const char *name)
+{
+ struct _pvtuAdaptFileDialog{
+ Fl_Window *window;
+ Fl_Choice *c[2];
+ Fl_Button *ok, *cancel, *push[2];
+ Fl_Value_Input *vi[3];
+ Fl_Check_Button *defautName;
+ };
+ static _pvtuAdaptFileDialog *dialog = NULL;
+
+ static Fl_Menu_Item viewmenu[] = {
+ {"Current", 0, 0, 0},
+ {"Visible", 0, 0, 0},
+ {"All", 0, 0, 0},
+ {0}
+ };
+ static Fl_Menu_Item formatmenu[] = {
+ {"Binary", 0, 0, 0},
+ {"ASCII", 0, 0, 0},
+ {0}
+ };
+
+ int BBB = BB + 9; // labels too long
+
+ if(!dialog){
+ dialog = new _pvtuAdaptFileDialog;
+ int h = 3 * WB + 3 * BH, w = 2 * BBB + 3 * WB, y = WB;
+ dialog->window = new Fl_Double_Window(2*w, 2.5*h, "Adaptive View Options"); //FIXME: dimensions of the window
+ dialog->window->box(GMSH_WINDOW_BOX);
+ dialog->window->set_modal();
+ dialog->c[0] = new Fl_Choice(WB, y, BBB + BBB / 2, BH, "View(s)"); y += BH;
+ dialog->c[0]->menu(viewmenu);
+ dialog->c[0]->align(FL_ALIGN_RIGHT);
+ dialog->c[1] = new Fl_Choice(WB, y, BBB + BBB / 2, BH, "Format"); y += BH;
+ dialog->c[1]->menu(formatmenu);
+ dialog->c[1]->align(FL_ALIGN_RIGHT);
+
+ dialog->vi[0] = new Fl_Value_Input(WB, y, BBB + BBB / 2, BH, "Adaptive recursion level"); y += BH;
+ dialog->vi[0]->align(FL_ALIGN_RIGHT);
+ dialog->vi[0]->minimum(0);
+ dialog->vi[0]->maximum(6);
+ dialog->vi[0]->step(1);
+ dialog->vi[0]->value(1);
+ dialog->vi[0]->when(FL_WHEN_RELEASE);
+
+ dialog->vi[1] = new Fl_Value_Input(WB, y, BBB + BBB / 2, BH, "Target error"); y += BH;
+ dialog->vi[1]->align(FL_ALIGN_RIGHT);
+ dialog->vi[1]->minimum(-1.e-4);
+ dialog->vi[1]->maximum(0.1);
+ dialog->vi[1]->step(1.e-4);
+ dialog->vi[1]->value(-1.e-4);
+ dialog->vi[1]->when(FL_WHEN_RELEASE);
+
+ dialog->vi[2] = new Fl_Value_Input(WB, y, BBB + BBB / 2, BH, "Number of viz parts"); y += BH;
+ dialog->vi[2]->align(FL_ALIGN_RIGHT);
+ dialog->vi[2]->minimum(1);
+ dialog->vi[2]->maximum(262144);
+ dialog->vi[2]->step(1);
+ dialog->vi[2]->value(4);
+ dialog->vi[2]->when(FL_WHEN_RELEASE);
+
+ dialog->defautName = new Fl_Check_Button(WB, y, BBB + BBB / 2, BH, "Use default filename (recommended)"); y += BH;
+ dialog->defautName->value(1);
+
+
+ dialog->ok = new Fl_Return_Button(WB, y + WB, BBB, BH, "OK");
+ dialog->cancel = new Fl_Button(2 * WB + BBB, y + WB, BBB, BH, "Cancel");
+ dialog->window->end();
+ dialog->window->hotspot(dialog->window);
+ }
+
+ dialog->window->show();
+
+ while(dialog->window->shown()){
+ Fl::wait();
+ for (;;) {
+ Fl_Widget* o = Fl::readqueue();
+ if (!o) break;
+ if (o == dialog->ok) {
+ bool isBinary;
+ switch(dialog->c[1]->value()){
+ case 0: isBinary = true; break;
+ case 1: isBinary = false; break;
+ }
+
+ // Only one view can currently be saved at a time in a pvtu file set, with a repetition of the topology structure.
+ // Views/Fields can then be appended in ParaView using the AppendAttributes filter
+ // bool canAppend = (format == 2) ? true : false;
+
+ int adaptLev = dialog->vi[0]->value();
+ double adaptErr = dialog->vi[1]->value();
+ int npart = dialog->vi[2]->value();
+ int useDefaultName = dialog->defautName->value();
+ bool canAppend = false; // Not yet implemented for VTK format here due to a tradeoff
+ // to limit memory consumption for high levels of adaptation
+ _saveAdaptedViews(name, useDefaultName, dialog->c[0]->value(), isBinary, adaptLev, adaptErr, npart, canAppend);
+ dialog->window->hide();
+ return 1;
+ }
+ if (o == dialog->window || o == dialog->cancel){
+ dialog->window->hide();
+ return 0;
+ }
+ }
+ }
+ return 0;
+}
+
int genericViewFileDialog(const char *name, const char *title, int format)
{
struct _viewFileDialog{
diff --git a/Fltk/fileDialogs.h b/Fltk/fileDialogs.h
index 557b566d35..393c71053e 100644
--- a/Fltk/fileDialogs.h
+++ b/Fltk/fileDialogs.h
@@ -29,6 +29,7 @@ int pgfBitmapFileDialog(const char *filename, const char *title, int format);
int genericMeshFileDialog(const char *filename, const char *title, int format,
bool binary_support, bool element_tag_support);
int posFileDialog(const char *name);
+int pvtuAdaptFileDialog(const char *name);
int genericViewFileDialog(const char *name, const char *title, int format);
int gl2psFileDialog(const char *filename, const char *title, int format);
int optionsFileDialog(const char *filename);
diff --git a/Fltk/graphicWindow.cpp b/Fltk/graphicWindow.cpp
index 8a382e2383..1c43e2652b 100644
--- a/Fltk/graphicWindow.cpp
+++ b/Fltk/graphicWindow.cpp
@@ -318,6 +318,7 @@ static int _save_svg(const char *name){ return gl2psFileDialog
static int _save_yuv(const char *name){ return genericBitmapFileDialog
(name, "YUV Options", FORMAT_YUV); }
static int _save_view_pos(const char *name){ return posFileDialog(name); }
+static int _save_view_adapt_pvtu(const char *name){ return pvtuAdaptFileDialog(name); }
static int _save_view_med(const char *name){ return genericViewFileDialog
(name, "MED Options", 6); }
static int _save_view_txt(const char *name){ return genericViewFileDialog
@@ -330,6 +331,7 @@ static int _save_auto(const char *name)
switch(GuessFileFormatFromFileName(name)){
case FORMAT_MSH : return _save_msh(name);
case FORMAT_POS : return _save_view_pos(name);
+ case FORMAT_PVTU : return _save_view_adapt_pvtu(name);
case FORMAT_TXT : return _save_view_txt(name);
case FORMAT_OPT : return _save_options(name);
case FORMAT_GEO : return _save_geo(name);
@@ -413,6 +415,7 @@ static void file_save_as_cb(Fl_Widget *w, void *data)
#endif
{"Post-processing - Generic TXT" TT "*.txt", _save_view_txt},
{"Post-processing - Mesh Statistics" TT "*.pos", _save_mesh_stat},
+ {"Post-processing - Adapted data" TT "*.pvtu", _save_view_adapt_pvtu},
{"Image - Encapsulated PostScript" TT "*.eps", _save_eps},
{"Image - GIF" TT "*.gif", _save_gif},
#if defined(HAVE_LIBJPEG)
diff --git a/Post/PView.h b/Post/PView.h
index 01f3577483..64a15f3afe 100644
--- a/Post/PView.h
+++ b/Post/PView.h
@@ -129,6 +129,11 @@ class PView{
static bool writeX3D(const std::string &fileName );
// IO write routine
bool write(const std::string &fileName, int format, bool append=false);
+
+ // michel.rasquin@cenaero.be:
+ // Routines for export of adapted views to pvtu file format for parallel visualization with paraview
+ bool writeAdapt(const std::string &fileName, int useDefaultName, bool isBinary,
+ int adaptLev, double adaptErr, int npart, bool append=false);
// vertex arrays to draw the elements efficiently
VertexArray *va_points, *va_lines, *va_triangles, *va_vectors, *va_ellipses;
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index c6820feaec..3cac8a552b 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -41,6 +41,34 @@ void PViewData::initAdaptiveData(int step, int level, double tol)
}
}
+void PViewData::initAdaptiveDataLight(int step, int level, double tol)
+{
+ if(!_adaptive){
+ Msg::Info("Initializing adaptive data %p interp size=%d", this, _interpolation.size());
+ // michel.rasquin@cenaero.be:
+ // _outData in adaptive.h is only used for visualization of adapted views in the GMSH GUI.
+ // In some cases (export of adapted views under pvtu format, use of GMSH as external lib),
+ // this object is not needed so avoid its allocation in order to limit memory consumption
+ bool outDataInit = false;
+ _adaptive = new adaptiveData(this,outDataInit);
+ }
+}
+
+void PViewData::saveAdaptedViewForVTK(const std::string &guifileName, int useDefaultName,
+ int step, int level, double tol, int npart, bool isBinary)
+{
+ if(_adaptive) {
+ // _adaptiveData has already been allocated from the adaptive view panel of the GUI for instance.
+ _adaptive->changeResolutionForVTK(step, level, tol, npart, isBinary, guifileName, useDefaultName);
+ }
+ else {
+ initAdaptiveDataLight(step, level, tol);
+ _adaptive->changeResolutionForVTK(step, level, tol, npart, isBinary, guifileName, useDefaultName);
+ destroyAdaptiveData();
+ }
+}
+
+
void PViewData::destroyAdaptiveData()
{
if(_adaptive) delete _adaptive;
diff --git a/Post/PViewData.h b/Post/PViewData.h
index e9e8a45a05..bdaae1e2ea 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -46,6 +46,8 @@ class PViewData {
interpolationMatrices _interpolation;
// global map of "named" interpolation matrices
static std::map<std::string, interpolationMatrices> _interpolationSchemes;
+ // string for the name of the interpolation scheme
+ std::string _interpolationSchemeName;
public:
PViewData();
@@ -196,6 +198,15 @@ class PViewData {
// initialize/destroy adaptive data
void initAdaptiveData(int step, int level, double tol);
+
+ // michel.rasquin@cenaero.be:
+ // Routines for
+ // - export of adapted views to pvtu file format for parallel visualization with paraview,
+ // - and/or generation of VTK data structure for ParaView plugin.
+ void initAdaptiveDataLight(int step, int level, double tol);
+ void saveAdaptedViewForVTK(const std::string &guifileName, int useDefaultName,
+ int step, int level, double tol, int npart, bool isBinary);
+
void destroyAdaptiveData();
// return the adaptive data
@@ -218,6 +229,9 @@ class PViewData {
static void removeInterpolationScheme(const std::string &name);
static void addMatrixToInterpolationScheme(const std::string &name, int type,
fullMatrix<double> &mat);
+ static int getSizeInterpolationScheme() {return _interpolationSchemes.size();}
+ std::string getInterpolationSchemeName() {return _interpolationSchemeName;}
+ void setInterpolationSchemeName(std::string name) {_interpolationSchemeName = name;}
// smooth the data in the view (makes it C0)
virtual void smooth();
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 2870f23b1e..7225b35b82 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -127,6 +127,10 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
if(!haveInterpolationMatrices()){
GModel *model = _steps[0]->getModel();
+
+ // michel.rasquin@cenaero.be:
+ // Required for ParaView plugin linked to GMSH as external library.
+ setInterpolationSchemeName(interpolationScheme);
// if an interpolation scheme is explicitly provided, use it
if(interpolationScheme.size()){
diff --git a/Post/PViewIO.cpp b/Post/PViewIO.cpp
index f441b89398..0435b91ae7 100644
--- a/Post/PViewIO.cpp
+++ b/Post/PViewIO.cpp
@@ -14,6 +14,7 @@
#include "StringUtils.h"
#include "Context.h"
#include "OS.h"
+#include "adaptiveData.h"
bool PView::readPOS(const std::string &fileName, int fileIndex)
{
@@ -331,5 +332,12 @@ bool PView::write(const std::string &fileName, int format, bool append)
return ret;
}
-
-
+// michel.rasquin@cenaero.be:
+// Routines for export of adapted views to pvtu file format for parallel visualization with paraview.
+bool PView::writeAdapt(const std::string &guifileName, int useDefaultName, bool isBinary,
+ int adaptLev, double adaptErr, int npart, bool append)
+{
+ Msg::StatusBar(true, "Writing '%s'...", guifileName.c_str());
+ _data->saveAdaptedViewForVTK(guifileName, useDefaultName, _options->timeStep, adaptLev, adaptErr, npart, isBinary);
+ return true;
+}
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index dfe854b6b7..d94b795c20 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -50,6 +50,11 @@ int adaptiveTetrahedron::numEdges = 6;
int adaptiveHexahedron::numEdges = 12;
int adaptivePyramid::numEdges = 8;
+std::vector<vectInt> globalVTKData::vtkGlobalConnectivity;
+std::vector<int> globalVTKData::vtkGlobalCellType;
+std::vector<PCoords> globalVTKData::vtkGlobalCoords;
+std::vector<PValues> globalVTKData::vtkGlobalValues;
+
template <class T>
static void cleanElement()
{
@@ -1497,13 +1502,18 @@ void adaptiveElements<T>::addInView(double tol, int step,
}
}
-adaptiveData::adaptiveData(PViewData *data)
+adaptiveData::adaptiveData(PViewData *data, bool outDataInit)
: _step(-1), _level(-1), _tol(-1.), _inData(data),
_points(0), _lines(0), _triangles(0), _quadrangles(0),
_tetrahedra(0), _hexahedra(0), _prisms(0),_pyramids(0)
{
- _outData = new PViewDataList(true);
- _outData->setName(data->getName() + "_adapted");
+ if(outDataInit == true) { // For visualization of the adapted view in GMSH GUI only
+ _outData = new PViewDataList(true);
+ _outData->setName(data->getName() + "_adapted");
+ }
+ else {
+ _outData = 0; // For external used
+ }
std::vector<fullMatrix<double>*> p;
if(_inData->getNumPoints()){
_inData->getInterpolationMatrices(TYPE_PNT, p);
@@ -1537,6 +1547,8 @@ adaptiveData::adaptiveData(PViewData *data)
_inData->getInterpolationMatrices(TYPE_PYR, p);
_pyramids = new adaptiveElements<adaptivePyramid>(p);
}
+ upWriteVTK(true); // By default, write VTK data if called...
+ upBuildStaticData(false); // ... and do not generated global static data structure (only useful for ParaView plugin).
}
adaptiveData::~adaptiveData()
@@ -1549,7 +1561,7 @@ adaptiveData::~adaptiveData()
if(_prisms) delete _prisms;
if(_hexahedra) delete _hexahedra;
if(_pyramids) delete _pyramids;
- delete _outData;
+ if(_outData) delete _outData;
}
double adaptiveData::timerInit = 0.;
@@ -1591,3 +1603,1023 @@ void adaptiveData::changeResolution(int step, int level, double tol,
printf("adapt time = %g\n", timerAdapt);
#endif
}
+
+// michel.rasquin@cenaero.be:
+// Routines for
+// - export of adapted views to pvtu file format for parallel visualization with paraview,
+// - and/or generation of VTK data structure for ParaView plugin.
+
+bool VTKData::isLittleEndian()
+{
+ int num = 1;
+ if(*(char *)&num == 1)
+ return true; // Little Endian
+ else
+ return false; // Big Endian
+}
+
+void VTKData::SwapArrayByteOrder( void* array, int nbytes, int nItems )
+{
+ // This swaps the byte order for the array of nItems each of size nbytes
+ int i,j;
+ unsigned char* ucDst = (unsigned char*)array;
+
+ for(i=0; i < nItems; i++) {
+ for(j=0; j < (nbytes/2); j++)
+ std::swap( ucDst[j] , ucDst[(nbytes - 1) - j] );
+ ucDst += nbytes;
+ }
+}
+
+void VTKData::writeVTKElmData()
+{
+ // This routine writes vtu files (ascii or binary) from a elemental data base of nodes coordinates,
+ // cell connectivity, type and offset, and point data (either scalar or vector field)
+
+ // Format choice
+ if(vtkFormat == "vtu") {
+
+ if(vtkCountTotElmLev0 <= numPartMinElm*minElmPerPart) {
+ if( (vtkCountTotElmLev0-1)%minElmPerPart == 0) { //new filename
+ vtkCountFile = (vtkCountTotElmLev0-1)/minElmPerPart;
+ initVTKFile();
+ }
+ }
+ else {
+ if( (vtkCountTotElmLev0-1-numPartMinElm*minElmPerPart)%maxElmPerPart == 0) { //new filename
+ vtkCountFile = numPartMinElm+(vtkCountTotElmLev0-1-numPartMinElm*minElmPerPart)/maxElmPerPart;
+ initVTKFile();
+ }
+ }
+
+ if (vtkIsBinary == true) { // Use appended format for raw binary
+
+ // Write raw binary data to separate files first.
+ // Text headers will be added later, as wall as raw data size (needs to know the size before)
+
+ int counter;
+ uint64_t *i64array;
+ uint8_t *i8array;
+ double *darray;
+
+ // Node value
+ counter = 0;
+ darray= new double[vtkNumComp*vtkLocalValues.size()];
+ for(std::vector<PValues>::iterator it = vtkLocalValues.begin();it != vtkLocalValues.end(); ++it) {
+ for(int i=0;i<vtkNumComp;i++) {
+ darray[counter+i] = it->v[i];
+ }
+ counter+=vtkNumComp;
+ vtkCountTotVal+=vtkNumComp;
+ }
+ assert(counter==vtkNumComp* (int) vtkLocalValues.size());
+ fwrite(darray,sizeof(double),vtkNumComp*vtkLocalValues.size(),vtkFileNodVal);
+ delete[] darray;
+
+ // Points
+ int sizeArray = (int) vtkLocalCoords.size();
+ darray = new double[3*sizeArray];
+ counter = 0;
+ for(std::vector<PCoords>::iterator it = vtkLocalCoords.begin();it != vtkLocalCoords.end(); ++it) {
+ for(int i=0;i<3;i++) {
+ darray[counter+i] = (*it).c[i];
+ }
+ counter+=3;
+ vtkCountCoord+=3;
+ }
+ fwrite(darray, sizeof(double),3*sizeArray,vtkFileCoord);
+ delete[] darray;
+
+ //Cells
+
+ // First count the number of integers that described the cell data in vtkConnectivity
+ // See http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf (page 4)
+ int cellSizeData = 0;
+ for(std::vector<vectInt>::iterator it = vtkLocalConnectivity.begin();it != vtkLocalConnectivity.end(); ++it) {
+ cellSizeData+= (int) it->size(); // Contrary to vtk format, no +1 required for the number of nodes in the element
+ }
+
+ // Connectivity (and build offset at the same time)
+ counter = 0;
+ int cellcounter = 0;
+ i64array = new uint64_t[cellSizeData];
+ uint64_t *cellOffset = new uint64_t[vtkLocalConnectivity.size()];
+ for(std::vector<vectInt>::iterator it = vtkLocalConnectivity.begin();it != vtkLocalConnectivity.end(); ++it) {
+ for(vectInt::iterator jt = it->begin(); jt != it->end(); ++jt) {
+ i64array[counter] = *jt;
+ counter++;
+ vtkCountTotNodConnect++;
+ }
+ cellOffset[cellcounter] = vtkCountTotNodConnect; // build the offset
+ cellcounter++;
+ }
+ fwrite(i64array,sizeof(uint64_t),cellSizeData,vtkFileConnect);
+ delete[] i64array;
+
+ // Cell offset
+ fwrite(cellOffset,sizeof(uint64_t),vtkLocalConnectivity.size(),vtkFileCellOffset);
+ delete[] cellOffset;
+
+ // Cell type
+ counter = 0;
+ i8array = new uint8_t[vtkLocalConnectivity.size()];
+ for(std::vector<int>::iterator it = vtkLocalCellType.begin();it != vtkLocalCellType.end(); it++) {
+ i8array[counter] = *it;
+ counter++;
+ }
+ fwrite(i8array,sizeof(uint8_t),vtkLocalConnectivity.size(),vtkFileCellType);
+ delete[] i8array;
+
+ }
+ else { // ascii
+
+ // Node values
+ for(std::vector<PValues>::iterator it = vtkLocalValues.begin();it != vtkLocalValues.end(); ++it) {
+
+ for(int i=0;i<vtkNumComp;i++) {
+ fprintf(vtkFileNodVal,"%23.16e ",(*it).v[i]);
+ vtkCountTotVal++;
+ if(vtkCountTotVal%6 == 0) fprintf(vtkFileNodVal,"\n");
+ }
+ }
+
+ for(std::vector<PCoords>::iterator it = vtkLocalCoords.begin();it != vtkLocalCoords.end(); it++) {
+ fprintf(vtkFileCoord,"%23.16e %23.16e %23.16e ",(*it).c[0],(*it).c[1],(*it).c[2]);
+ vtkCountCoord+=3;
+ if(vtkCountCoord%6 == 0) fprintf(vtkFileCoord,"\n");
+ }
+
+ // Cells
+ // Connectivity
+ int *cellOffset = new int[vtkLocalConnectivity.size()];
+ int cellcounter = 0;
+ for(std::vector<vectInt>::iterator it = vtkLocalConnectivity.begin();it != vtkLocalConnectivity.end(); ++it) {
+ for(vectInt::iterator jt = it->begin(); jt != it->end(); ++jt) {
+ fprintf(vtkFileConnect,"%d ",*jt);
+ vtkCountTotNodConnect++;
+ if(vtkCountTotNodConnect%6 == 0) fprintf(vtkFileConnect,"\n");
+ }
+ cellOffset[cellcounter] = vtkCountTotNodConnect; // build the offset
+ cellcounter++;
+ }
+
+ // Cell offset
+ for(uint64_t i = 0 ; i<vtkLocalConnectivity.size(); i++) {
+ fprintf(vtkFileCellOffset,"%d ",cellOffset[i]);
+ vtkCountCellOffset++;
+ if(vtkCountCellOffset%6 == 0) fprintf(vtkFileCellOffset,"\n");
+ }
+ delete[] cellOffset;
+
+ // Cell type
+ for(std::vector<int>::iterator it = vtkLocalCellType.begin();it != vtkLocalCellType.end(); it++) {
+ fprintf(vtkFileCellType,"%d ",*it);
+ vtkCountCellType++;
+ if(vtkCountCellType%6 == 0) fprintf(vtkFileCellType,"\n");
+ }
+
+ } //if ascii
+
+ //finalize and close current vtu file
+ if(vtkCountTotElmLev0 <= numPartMinElm*minElmPerPart) {
+ if( vtkCountTotElmLev0%minElmPerPart == 0) {
+ finalizeVTKFile();
+ }
+ }
+ else {
+ if( (vtkCountTotElmLev0-numPartMinElm*minElmPerPart)%maxElmPerPart == 0) {
+ finalizeVTKFile();
+ }
+ }
+
+ } // vtu format
+ else
+ printf("ERROR: format unknown\n");
+
+ clearLocalData();
+}
+
+void VTKData::initVTKFile()
+{
+ // Temporary files
+ vtkFileCoord = fopen("vtkCoords.vtu","wb");
+ vtkFileConnect = fopen("vtkConnectivity.vtu","wb");
+ vtkFileCellOffset = fopen("vtkCellOffset.vtu","wb");
+ vtkFileCellType = fopen("vtkCellType.vtu","wb");
+ vtkFileNodVal = fopen("vtkNodeValue.vtu","wb");
+
+ if(vtkCountFile == 0) { // write the pvtu file and create the corresponding directory for vtu files
+
+ if (vtkUseDefaultName == 1) {
+ vtkDirName = vtkFieldName
+ + "_step" + ToString<int>(vtkStep)
+ + "_level" + ToString<int>(vtkLevel)
+ + "_tol" + ToString<double>(vtkTol)
+ + "_npart" + ToString<int>(vtkNpart);
+ }
+ else {
+ // Remove existing extension here to avoid duplicate
+ std::size_t found = vtkFileName.find_last_of(".");
+ if(found != std::string::npos) vtkFileName = vtkFileName.substr(0,found); //remove extension
+ vtkDirName = vtkFileName;
+ }
+
+ mkdir(vtkDirName.c_str(),S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); //755
+
+ vtkFileName = vtkDirName + ".p" + vtkFormat; // add pvtu extension to file name
+ vtkFile = fopen(vtkFileName.c_str(),"w");
+
+ bool littleEndian = isLittleEndian(); // Determine endianess
+ if (littleEndian == true)
+ fprintf(vtkFile,"<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\">\n");
+ else
+ fprintf(vtkFile,"<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" byte_order=\"BigEndian\">\n");
+
+ fprintf(vtkFile,"<PUnstructuredGrid GhostLevel=\"0\">\n");
+ fprintf(vtkFile,"<PPoints>\n");
+ fprintf(vtkFile,"<DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\"/>\n");
+ fprintf(vtkFile,"</PPoints>\n");
+
+ fprintf(vtkFile,"<PCells>\n");
+ fprintf(vtkFile,"<PDataArray type=\"Int64\" Name=\"connectivity\" NumberOfComponents=\"1\"/>\n");
+ fprintf(vtkFile,"<PDataArray type=\"Int64\" Name=\"offsets\" NumberOfComponents=\"1\"/>\n");
+ fprintf(vtkFile,"<PDataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\"/>\n");
+ fprintf(vtkFile,"</PCells>\n");
+
+ fprintf(vtkFile,"<PPointData>\n");
+ fprintf(vtkFile,"<PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"%d\"/>\n",vtkFieldName.c_str(),vtkNumComp);
+ fprintf(vtkFile,"</PPointData>\n");
+
+ fprintf(vtkFile,"<PCellData>\n");
+ fprintf(vtkFile,"</PCellData>\n");
+
+ for(int i=0; i<vtkNpart; i++)
+ fprintf(vtkFile,"<Piece Source=\"%s/data%d.vtu\"/>\n",vtkDirName.c_str(),i);
+ fprintf(vtkFile,"</PUnstructuredGrid>\n");
+ fprintf(vtkFile,"</VTKFile>\n");
+ fclose(vtkFile);
+ }
+}
+
+void VTKData::finalizeVTKFile()
+{
+ // This routine writes vtu files (ascii or binary) from a complete data base of nodes coordinates,
+ // cell connectivity, type and offset, and point data (either scalar or vector field)
+
+ // Close first temporary files.
+ // Todo: Avoid multiple open/close actions and keep the file open to write all information
+ fclose(vtkFileCoord);
+ fclose(vtkFileConnect);
+ fclose(vtkFileCellOffset);
+ fclose(vtkFileCellType);
+ fclose(vtkFileNodVal);
+
+ bool littleEndian = isLittleEndian(); // Determine endianess
+
+ // Open final file
+ std::string filename;
+ filename = vtkDirName + "/data" + ToString(vtkCountFile) + "." + vtkFormat;
+
+ Msg::StatusBar(true,"Writing VTK data in %s: fieldname = %s - numElm = %d - numNod = %d nodes\n",
+ filename.c_str(), vtkFieldName.c_str(), vtkCountTotElm, vtkCountTotNod);
+
+ assert(vtkCountTotNod == vtkCountCoord/3);
+
+ // Now concatenate headers with data files
+ if(vtkFormat == "vtu") { // Format choice
+
+ if (vtkIsBinary == true) { // Binary or ascii
+
+ vtkFile = fopen(filename.c_str(),"wb");
+ if(vtkFile == NULL) {
+ printf("Could not open file %s\n", filename.c_str());
+ return;
+ }
+
+ uint64_t byteoffset = 0;
+
+ // Headers first
+
+ if (littleEndian == true)
+ fprintf(vtkFile,"<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\n");
+ else
+ fprintf(vtkFile,"<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" byte_order=\"BigEndian\" header_type=\"UInt64\">\n");
+ fprintf(vtkFile,"<UnstructuredGrid>\n");
+ fprintf(vtkFile,"<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",vtkCountTotNod,vtkCountTotElm);
+
+
+ // Node value
+ fprintf(vtkFile,"<PointData>\n");
+ fprintf(vtkFile,"<DataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"%d\" format=\"appended\" offset=\"%" PRIu64 "\"/>\n",vtkFieldName.c_str(),vtkNumComp,byteoffset);
+ fprintf(vtkFile,"</PointData>\n");
+ byteoffset = byteoffset + (vtkCountTotNod*vtkNumComp+1)*sizeof(double); // +1 for datasize in bytes
+
+ // Cell values (none here but may change)
+ fprintf(vtkFile,"<CellData>\n");
+ fprintf(vtkFile,"</CellData>\n"); // no offset here because empty cell data
+
+ // Nodes
+ fprintf(vtkFile,"<Points>\n");
+ fprintf(vtkFile,"<DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PRIu64 "\"/>\n",byteoffset);
+ fprintf(vtkFile,"</Points>\n");
+ byteoffset = byteoffset + (vtkCountCoord+1)*sizeof(double); // +1 for datasize in bytes
+
+ // Cells
+ fprintf(vtkFile,"<Cells>\n");
+ fprintf(vtkFile,"<DataArray type=\"Int64\" Name=\"connectivity\" format=\"appended\" offset=\"%" PRIu64 "\"/>\n",byteoffset);
+ byteoffset = byteoffset + (vtkCountTotNodConnect+1)*sizeof(uint64_t);
+ fprintf(vtkFile,"<DataArray type=\"Int64\" Name=\"offsets\" format=\"appended\" offset=\"%" PRIu64 "\"/>\n",byteoffset);
+ byteoffset = byteoffset + (vtkCountTotElm+1)*sizeof(uint64_t);
+ fprintf(vtkFile,"<DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"%" PRIu64 "\"/>\n",byteoffset);
+ byteoffset = byteoffset + (vtkCountTotElm+1)*sizeof(uint8_t);
+ fprintf(vtkFile,"</Cells>\n");
+
+ fprintf(vtkFile,"</Piece>\n");
+ fprintf(vtkFile,"</UnstructuredGrid>\n");
+
+ fprintf(vtkFile,"<AppendedData encoding=\"raw\">\n");
+ fprintf(vtkFile,"_");
+
+ uint64_t datasize;
+
+ // Node values
+ datasize = vtkNumComp*vtkCountTotNod*sizeof(double);
+ fwrite(&datasize,sizeof(uint64_t),1,vtkFile);
+ fclose(vtkFile);
+
+ std::ifstream if_vtkNodeValue("vtkNodeValue.vtu", std::ios_base::binary);
+ std::ofstream of_vtkfile(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkNodeValue.rdbuf();
+ if_vtkNodeValue.close();
+ of_vtkfile.close();
+
+ // Points
+ vtkFile = fopen(filename.c_str(),"ab");
+ datasize = vtkCountTotNod*3*sizeof(double);
+ fwrite(&datasize,sizeof(uint64_t),1,vtkFile);
+ fclose(vtkFile);
+
+ std::ifstream if_vtkCoords("vtkCoords.vtu", std::ios_base::binary);
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkCoords.rdbuf();
+ if_vtkCoords.close();
+ of_vtkfile.close();
+
+ // Cells
+ // Connectivity
+ vtkFile = fopen(filename.c_str(),"ab");
+ datasize = vtkCountTotNodConnect*sizeof(uint64_t);
+ fwrite(&datasize,sizeof(uint64_t),1,vtkFile);
+ fclose(vtkFile);
+
+ std::ifstream if_vtkConnectivity("vtkConnectivity.vtu", std::ios_base::binary);
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkConnectivity.rdbuf();
+ if_vtkConnectivity.close();
+ of_vtkfile.close();
+
+ // Cell offset
+ vtkFile = fopen(filename.c_str(),"ab");
+ datasize = vtkCountTotElm*sizeof(uint64_t);
+ fwrite(&datasize,sizeof(uint64_t),1,vtkFile);
+ fclose(vtkFile);
+
+ std::ifstream if_vtkCellOffset("vtkCellOffset.vtu", std::ios_base::binary);
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkCellOffset.rdbuf();
+ if_vtkCellOffset.close();
+ of_vtkfile.close();
+
+ // Cell type
+ vtkFile = fopen(filename.c_str(),"ab");
+ datasize = vtkCountTotElm*sizeof(uint8_t);
+ fwrite(&datasize,sizeof(uint64_t),1,vtkFile);
+ fclose(vtkFile);
+
+ std::ifstream if_vtkCellType("vtkCellType.vtu", std::ios_base::binary);
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkCellType.rdbuf();
+ if_vtkCellType.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"ab");
+ fprintf(vtkFile,"\n");
+ fprintf(vtkFile,"</AppendedData>\n");
+ fprintf(vtkFile,"</VTKFile>\n"); // for both binary and ascii
+ fclose(vtkFile);
+
+ }
+ else { // ascii
+
+ vtkFile = fopen(filename.c_str(),"w");
+ if(vtkFile == NULL) {
+ printf("Could not open file %s\n", filename.c_str());
+ return;
+ }
+
+ if (littleEndian == true)
+ fprintf(vtkFile,"<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\n");
+ else
+ fprintf(vtkFile,"<VTKFile type=\"PUnstructuredGrid\" version=\"1.0\" byte_order=\"BigEndian\" header_type=\"UInt64\">\n");
+ fprintf(vtkFile,"<UnstructuredGrid>\n");
+ fprintf(vtkFile,"<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",vtkCountTotNod,vtkCountTotElm);
+
+ // Node values
+ fprintf(vtkFile,"<PointData>\n");
+ fprintf(vtkFile,"<DataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"%d\" format=\"ascii\">\n",vtkFieldName.c_str(),vtkNumComp);
+ fclose(vtkFile); //close file for binary concatenation
+
+ std::ifstream if_vtkNodeValue("vtkNodeValue.vtu", std::ios_base::binary);
+ std::ofstream of_vtkfile(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ of_vtkfile << if_vtkNodeValue.rdbuf();
+ if_vtkNodeValue.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"a");
+ fprintf(vtkFile,"</DataArray>\n");
+ fprintf(vtkFile,"</PointData>\n");
+
+ //Cell values
+ fprintf(vtkFile,"<CellData>\n");
+ fprintf(vtkFile,"</CellData>\n");
+
+ //Nodes
+ fprintf(vtkFile,"<Points>\n");
+ fprintf(vtkFile,"<DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\n");
+ fclose(vtkFile); //close file for binary concatenation
+
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ std::ifstream if_vtkCoords("vtkCoords.vtu", std::ios_base::binary);
+ of_vtkfile << if_vtkCoords.rdbuf();
+ if_vtkCoords.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"a");
+ fprintf(vtkFile,"</DataArray>\n");
+ fprintf(vtkFile,"</Points>\n");
+
+ // Cells
+ fprintf(vtkFile,"<Cells>\n");
+ fprintf(vtkFile,"<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n");
+ fclose(vtkFile); //close file for binary concatenation
+
+ // Connectivity
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ std::ifstream if_vtkConnectivity("vtkConnectivity.vtu", std::ios_base::binary);
+ of_vtkfile << if_vtkConnectivity.rdbuf();
+ if_vtkConnectivity.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"a");
+ fprintf(vtkFile,"</DataArray>\n");
+
+ // Cell offset
+ fprintf(vtkFile,"<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n");
+ fclose(vtkFile); //close file for binary concatenation
+
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ std::ifstream if_vtkCellOffset("vtkCellOffset.vtu", std::ios_base::binary);
+ of_vtkfile << if_vtkCellOffset.rdbuf();
+ if_vtkCellOffset.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"a");
+ fprintf(vtkFile,"</DataArray>\n");
+
+ // Cell type
+ fprintf(vtkFile,"<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
+ fclose(vtkFile); //close file for binary concatenation
+
+ of_vtkfile.open(filename.c_str(), std::ios_base::binary | std::ios_base::app);
+ std::ifstream if_vtkCellType("vtkCellType.vtu", std::ios_base::binary);
+ of_vtkfile << if_vtkCellType.rdbuf();
+ if_vtkCellType.close();
+ of_vtkfile.close();
+
+ vtkFile = fopen(filename.c_str(),"a");
+ fprintf(vtkFile,"</DataArray>\n");
+ fprintf(vtkFile,"</Cells>\n");
+
+ fprintf(vtkFile,"</Piece>\n");
+ fprintf(vtkFile,"</UnstructuredGrid>\n");
+
+
+ fprintf(vtkFile,"</VTKFile>\n"); // for both binary and ascii
+ fclose(vtkFile);
+ } // if binary/ascii
+
+ // Remove temporary files now
+ if(remove("vtkCoords.vtu") !=0 ) printf("ERROR: Could not remove vtkCoords.vtu\n");
+ if(remove("vtkConnectivity.vtu") !=0 ) printf("ERROR: Could not remove vtkConnectivity.vtu\n");
+ if(remove("vtkCellOffset.vtu") !=0 ) printf("ERROR: Could not remove vtkCellOffset.vtu\n");
+ if(remove("vtkCellType.vtu") !=0 ) printf("ERROR: Could not remove vtkCellType.vtu\n");
+ if(remove("vtkNodeValue.vtu") !=0 ) printf("ERROR: Could not remove vtkNodeValue.vtu\n");
+
+ // Reset counters for next file
+ vtkCountTotNod = 0;
+ vtkCountTotElm = 0;
+ vtkCountCoord = 0;
+ vtkCountTotNodConnect = 0;
+ vtkCountTotVal = 0;
+ vtkCountCellOffset = 0;
+ vtkCountCellType = 0;
+ }
+
+ else
+ Msg::Error("File format unknown: %s\n",vtkFormat.c_str());
+}
+
+int VTKData::getPVCellType(int numEdges)
+{
+ int cellType; // Convention for cell types in ParaView
+ switch(numEdges){
+ case 0:
+ printf("WARNING: Trying to write a node to the ParaView data base and file\n");
+ cellType = -1;
+ break;
+ case 1:
+ printf("WARNING: Trying to write a node to the ParaView data base and file\n");
+ cellType = -2;
+ break;
+ case 3:
+ cellType = 5; // 2D VTK triangle
+ break;
+ case 4:
+ cellType = 9; // 2D VTK quadrangle
+ break;
+ case 6:
+ cellType = 10; // 3D VTK tetrahedron
+ break;
+ case 9:
+ cellType = 13; // 3D VTK prism/wedge
+ break;
+ case 8:
+ cellType = 14; // 3D VTK pyramid
+ break;
+ case 12:
+ cellType = 12; // 3D VTK hexahedron
+ break;
+ default:
+ printf("ERROR: No cell type was detected\n");
+ cellType = -1;
+ break;
+ }
+
+ return cellType;
+}
+
+
+
+template <class T>
+void adaptiveElements<T>::adaptForVTK(double tol,
+ int numComp,
+ std::vector<PCoords> &coords,
+ std::vector<PValues> &values)
+{
+ int numVertices = T::allVertices.size();
+
+ if(!numVertices){
+ Msg::Error("No adapted vertices to interpolate");
+ return;
+ }
+
+ int numVals = _coeffsVal ? _coeffsVal->size1() : T::numNodes;
+ if(numVals != (int)values.size()){
+ Msg::Error("Wrong number of values in adaptation %d != %i",
+ numVals, values.size());
+ return;
+ }
+
+ #ifdef TIMER
+ double t1 = GetTimeInSeconds();
+ #endif
+
+ fullVector<double> val(numVals), res(numVertices);
+ switch (numComp) {
+ case 1:
+ {
+ for(int i = 0; i < numVals; i++) val(i) = values[i].v[0];
+ break;
+ }
+ case 3:
+ case 9:
+ {
+ for(int i = 0; i < numVals; i++) {
+ val(i) = 0;
+ for (int k=0; k < numComp; k++) val(i) += values[i].v[k] * values[i].v[k];
+ }
+ break;
+ }
+ default:
+ {
+ Msg::Error("Can only adapt scalar, vector or tensor data");
+ return;
+ }
+ }
+
+ _interpolVal->mult(val, res);
+
+ double minVal = VAL_INF;
+ double maxVal = -VAL_INF;
+ for(int i = 0; i < numVertices; i++){
+ minVal = std::min(minVal, res(i));
+ maxVal = std::max(maxVal, res(i));
+ }
+
+ fullMatrix<double> *resxyz = 0;
+ if(numComp == 3 || numComp == 9){
+ fullMatrix<double> valxyz(numVals, numComp);
+ resxyz = new fullMatrix<double>(numVertices, numComp);
+ for(int i = 0; i < numVals; i++){
+ for (int k = 0; k < numComp; k++) {
+ valxyz(i,k) = values[i].v[k];
+ }
+ }
+ _interpolVal->mult(valxyz, *resxyz);
+ }
+
+ int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
+ if(numNodes != (int)coords.size()){
+ Msg::Error("Wrong number of nodes in adaptation %d != %i",
+ numNodes, coords.size());
+ if(resxyz) delete resxyz;
+ return;
+ }
+
+ fullMatrix<double> xyz(numNodes, 3), XYZ(numVertices, 3);
+ for(int i = 0; i < numNodes; i++){
+ xyz(i, 0) = coords[i].c[0];
+ xyz(i, 1) = coords[i].c[1];
+ xyz(i, 2) = coords[i].c[2];
+ }
+ _interpolGeom->mult(xyz, XYZ);
+
+ #ifdef TIMER
+ adaptiveData::timerAdapt += GetTimeInSeconds() - t1;
+ return;
+ #endif
+
+ int i = 0;
+ for(std::set<adaptiveVertex>::iterator it = T::allVertices.begin(); it != T::allVertices.end(); ++it) {
+ // ok because we know this will not change the set ordering
+ adaptiveVertex *p = (adaptiveVertex*)&(*it);
+ p->val = res(i);
+ if(resxyz){
+ p->val = (*resxyz)(i, 0);
+ p->valy = (*resxyz)(i, 1);
+ p->valz = (*resxyz)(i, 2);
+ if (numComp == 9) {
+ p->valyx = (*resxyz)(i,3);
+ p->valyy = (*resxyz)(i,4);
+ p->valyz = (*resxyz)(i,5);
+ p->valzx = (*resxyz)(i,6);
+ p->valzy = (*resxyz)(i,7);
+ p->valzz = (*resxyz)(i,8);
+ }
+ }
+ p->X = XYZ(i, 0);
+ p->Y = XYZ(i, 1);
+ p->Z = XYZ(i, 2);
+ i++;
+ }
+
+ if(resxyz) delete resxyz;
+
+ for(typename std::list<T*>::iterator it = T::all.begin();it != T::all.end(); it++)
+ (*it)->visible = false;
+
+ if(tol != 0.){
+ double avg = fabs(maxVal - minVal);
+ if(tol < 0) avg = 1.; // force visibility to the smallest subdivision
+ T::error(avg, tol);
+ }
+
+ coords.clear();
+ values.clear();
+ for(typename std::list<T*>::iterator it = T::all.begin();it != T::all.end(); it++){
+ if((*it)->visible){
+ adaptiveVertex **p = (*it)->p;
+ for(int i = 0; i < T::numNodes; i++) {
+ coords.push_back(PCoords(p[i]->X, p[i]->Y, p[i]->Z));
+ switch (numComp) {
+ case 1:
+ values.push_back(PValues(p[i]->val));
+ break;
+ case 3:
+ values.push_back(PValues(p[i]->val, p[i]->valy, p[i]->valz));
+ break;
+ case 9:
+ values.push_back(PValues(p[i]->val,
+ p[i]->valy, p[i]->valz,
+ p[i]->valyx,p[i]->valyy,p[i]->valyz,
+ p[i]->valzx,p[i]->valzy,p[i]->valzz));
+ break;
+ }
+ }
+ }
+ }
+}
+
+
+template <class T>
+void adaptiveElements<T>::buildMapping(nodMap<T> &myNodMap, double tol, int &numNodInsert)
+{
+
+ if (tol > 0.0 || myNodMap.getSize() == 0) {
+ // Either this is not a uniform refinement and we need to rebuild the whole mapping
+ // for each canonical element, or this is the first time we try to build the mapping
+
+ myNodMap.cleanMapping(); // Required if tol > 0 (local error based adaptation)
+
+ for(typename std::list<T*>::iterator itleaf = T::all.begin();itleaf != T::all.end(); itleaf++) {
+ // Visit all the leaves of the refined canonical element
+
+ if ((*itleaf)->visible == true) {
+ // Find the leaves that are flagged for visibility
+
+ for(int i=0;i<T::numNodes;i++) {
+ // Visit each nodes of the leaf (3 for triangles, 4 for quadrangle, etc)
+ adaptiveVertex pquery;
+ pquery.x = (*itleaf)->p[i]->x;
+ pquery.y = (*itleaf)->p[i]->y;
+ pquery.z = (*itleaf)->p[i]->z;
+ std::set<adaptiveVertex>::iterator it = T::allVertices.find(pquery);
+ if(it == T::allVertices.end()){
+ printf("ERROR: Could not find adaptive Vertex in adaptiveElements<T>::buildMapping %f %f %f\n",
+ pquery.x,pquery.y,pquery.z);
+ }
+ else{
+ // Compute the distance in the list to get the mapping for
+ // the canonical element (note std:distance returns long int
+ int dist = (int) std::distance(T::allVertices.begin(),it);
+ myNodMap.mapping.push_back(dist);
+ }
+ // quit properly if vertex not found - Should not happen though
+ assert(it != T::allVertices.end());
+ } //for
+ } //if
+ }//for
+
+ // Count number of unique nodes from the mapping
+ // Use an ordered set for efficiency
+ // This set is also used in case of partiel refinement
+ std::set<int> uniqueNod;
+ for(std::vector<int>::iterator it = myNodMap.mapping.begin(); it != myNodMap.mapping.end(); it++) {
+ uniqueNod.insert(*it);
+ }
+ numNodInsert = (int) uniqueNod.size();
+
+ // Renumber the elm in the mapping in case of partial refinement (when vis tolerance > 0)
+ // so that we have a continuous numbering starting from 0 with no missing node id in the connectivity
+ // This require a new local and temporary mapping, based on uniqueNod already generated above
+ if(tol > 0.0) {
+ std::set<int>::iterator jt;
+ for(std::vector<int>::iterator it = myNodMap.mapping.begin();it != myNodMap.mapping.end();++it) {
+ jt = uniqueNod.find(*it);
+ *it = std::distance(uniqueNod.begin(), jt);
+ }
+ }
+ }
+}
+
+template <class T>
+void adaptiveElements<T>::addInViewForVTK(int step,
+ PViewData *in,
+ VTKData &myVTKData,
+ bool writeVTK,
+ bool buildStaticData)
+{
+ int numComp = in->getNumComponents(0, 0, 0);
+ if(numComp != 1 && numComp != 3 && numComp != 9) return;
+
+ int numEle = 0;
+ switch(T::numEdges){
+ case 0:
+ numEle = in->getNumPoints();
+ break;
+ case 1:
+ numEle = in->getNumLines();
+ break;
+ case 3:
+ numEle = in->getNumTriangles();
+ break;
+ case 4:
+ numEle = in->getNumQuadrangles();
+ break;
+ case 6:
+ numEle = in->getNumTetrahedra();
+ break;
+ case 9:
+ numEle = in->getNumPrisms();
+ break;
+ case 8:
+ numEle = in->getNumPyramids();
+ break;
+ case 12:
+ numEle = in->getNumHexahedra();
+ break;
+ }
+ if(!numEle) return;
+
+ // New variables for high order visualiztion through vtk files
+ int numNodInsert;
+ nodMap<T> myNodMap;
+
+ for(int ent = 0; ent < in->getNumEntities(step); ent++){
+ for(int ele = 0; ele < in->getNumElements(step, ent); ele++){
+ if(in->skipElement(step, ent, ele) ||
+ in->getNumEdges(step, ent, ele) != T::numEdges) continue;
+ int numNodes = in->getNumNodes(step, ent, ele);
+ std::vector<PCoords> coords;
+ for(int i = 0; i < numNodes; i++){
+ double x, y, z;
+ in->getNode(step, ent, ele, i, x, y, z);
+ coords.push_back(PCoords(x, y, z));
+ }
+ int numVal = in->getNumValues(step, ent, ele);
+ std::vector<PValues> values;
+
+ switch (numComp) {
+ case 1:
+ for(int i = 0; i < numVal; i++){
+ double val;
+ in->getValue(step, ent, ele, i, val);
+ values.push_back(PValues(val));
+ }
+ break;
+ case 3:
+ for(int i = 0; i < numVal / 3; i++){
+ double vx, vy, vz;
+ in->getValue(step, ent, ele, 3 * i, vx);
+ in->getValue(step, ent, ele, 3 * i + 1, vy);
+ in->getValue(step, ent, ele, 3 * i + 2, vz);
+ values.push_back(PValues(vx, vy, vz));
+ }
+ break;
+ case 9:
+ for(int i = 0; i < numVal / 9; i++){
+ double vxx, vxy, vxz,vyx, vyy, vyz,vzx, vzy, vzz ;
+ in->getValue(step, ent, ele, 9 * i + 0, vxx);
+ in->getValue(step, ent, ele, 9 * i + 1, vxy);
+ in->getValue(step, ent, ele, 9 * i + 2, vxz);
+ in->getValue(step, ent, ele, 9 * i + 3, vyx);
+ in->getValue(step, ent, ele, 9 * i + 4, vyy);
+ in->getValue(step, ent, ele, 9 * i + 5, vyz);
+ in->getValue(step, ent, ele, 9 * i + 6, vzx);
+ in->getValue(step, ent, ele, 9 * i + 7, vzy);
+ in->getValue(step, ent, ele, 9 * i + 8, vzz);
+ values.push_back(PValues(vxx,vxy,vxz,
+ vyx,vyy,vyz,
+ vzx,vzy,vzz));
+ }
+ break;
+ }
+
+ adaptForVTK(myVTKData.vtkTol, numComp, coords, values);//, out->Min, out->Max, plug);
+
+ // Inside initial element, after adapt() has been called
+
+ // Build the mapping of the canonical element,
+ // or recycle existing one in case of uniform refinement
+ buildMapping(myNodMap, myVTKData.vtkTol, numNodInsert);
+
+ // Pre-allocate some space for the local coordinates and connectivity
+ // in order to write to any component of the vector later through vec[i]
+
+ myVTKData.vtkLocalCoords.resize(numNodInsert,PCoords(0.0,0.0,0.0));
+ myVTKData.vtkLocalValues.resize(numNodInsert,PValues(numComp));
+
+ for(unsigned int i = 0; i < coords.size() / T::numNodes; i++){
+
+ // Loop over
+ // - all refined elements if refinement level > 0
+ // - single initial element when refinement box is checked for the first time (ref level =0)
+
+ // local connectivity for the considered sub triangle
+ vectInt vtkElmConnectivity;
+
+ for(int k = 0; k < T::numNodes; ++k) {
+
+ // Connectivity of the considered sub-element
+ int countTotNodloc = T::numNodes * i + k; // Nodes are duplicate here
+ int vtkNodeId = myVTKData.vtkCountTotNod + myNodMap.mapping[countTotNodloc];
+ vtkElmConnectivity.push_back(vtkNodeId);
+
+ // Coordinates of the nodes of the considered sub-element
+ double px, py, pz;
+ px = coords[T::numNodes * i + k].c[0];
+ py = coords[T::numNodes * i + k].c[1];
+ pz = coords[T::numNodes * i + k].c[2];
+ PCoords tmpCoords = PCoords(px,py,pz);
+ myVTKData.vtkLocalCoords[myNodMap.mapping[countTotNodloc]] = tmpCoords;
+
+ // Value associated with each nodes of the sub-element
+ myVTKData.vtkLocalValues[myNodMap.mapping[countTotNodloc]] = values[T::numNodes * i + k];
+
+ }
+
+ // Add elm connectivity to vector
+ myVTKData.vtkLocalConnectivity.push_back(vtkElmConnectivity);
+
+ // Increment global elm number
+ myVTKData.incrementTotElm(1);
+
+ // Save element type
+ myVTKData.vtkLocalCellType.push_back(myVTKData.getPVCellType(T::numEdges));
+
+
+ // Global variables
+ if(buildStaticData == true) {
+ globalVTKData::vtkGlobalConnectivity.push_back(vtkElmConnectivity);
+ globalVTKData::vtkGlobalCellType.push_back(myVTKData.getPVCellType(T::numEdges));
+ }
+
+ // Clear existing structure (safer)
+ vtkElmConnectivity.clear();
+ }
+
+ // Increment global node and elm-lev0 number
+ myVTKData.incrementTotNod(numNodInsert);
+ myVTKData.incrementTotElmLev0(1);
+
+ // Write the VTK data structure of the consider element to vtu file
+
+ if(writeVTK == true) {
+ myVTKData.writeVTKElmData();
+ }
+
+ if(buildStaticData == true) {
+
+ for(int i=0;i<numNodInsert; i++) {
+ globalVTKData::vtkGlobalCoords.push_back(myVTKData.vtkLocalCoords[i]);
+ }
+
+ for(int i=0;i<numNodInsert; i++) {
+ globalVTKData::vtkGlobalValues.push_back(myVTKData.vtkLocalValues[i]);
+ }
+ }
+
+ } // loop over mesh element
+ }
+
+}
+
+template <class T>
+int adaptiveElements<T>::countElmLev0(int step, PViewData *in)
+{
+ int sum = 0;
+ for(int ent = 0; ent < in->getNumEntities(step); ent++){
+ for(int ele = 0; ele < in->getNumElements(step, ent); ele++){
+ if(in->skipElement(step, ent, ele) ||
+ in->getNumEdges(step, ent, ele) != T::numEdges) continue;
+ else
+ sum++;
+ }
+ }
+ return sum;
+}
+
+int adaptiveData::countTotElmLev0(int step, PViewData *in)
+{
+
+ int sumElm = 0;
+
+ if(_triangles) sumElm+=_triangles->countElmLev0(step, in);
+ if(_quadrangles) sumElm+=_quadrangles->countElmLev0(step, in);
+ if(_tetrahedra) sumElm+=_tetrahedra->countElmLev0(step, in);
+ if(_prisms) sumElm+=_prisms->countElmLev0(step, in);
+ if(_hexahedra) sumElm+=_hexahedra->countElmLev0(step, in);
+ if(_pyramids) sumElm+=_pyramids->countElmLev0(step, in);
+
+ return sumElm;
+}
+
+
+void adaptiveData::changeResolutionForVTK(int step, int level, double tol, int npart, bool isBinary,
+ const std::string &guiFileName, int useDefaultName)
+{
+ //clean global VTK data structure before (re)generating it
+ if(buildStaticData == true) globalVTKData::clearGlobalData();
+
+ VTKData myVTKData(_inData->getName(), _inData->getNumComponents(0, 0, 0),
+ step, level, tol, guiFileName, useDefaultName, npart, isBinary);
+ myVTKData.vtkTotNumElmLev0 = countTotElmLev0(step, _inData);
+ myVTKData.setFileDistribution();
+
+ // Views of 2D and 3D elements only supported for VTK. _points and _lines are currently ignored.
+ if(_triangles) _triangles->init(myVTKData.vtkLevel);
+ if(_quadrangles) _quadrangles->init(myVTKData.vtkLevel);
+ if(_tetrahedra) _tetrahedra->init(myVTKData.vtkLevel);
+ if(_prisms) _prisms->init(myVTKData.vtkLevel);
+ if(_hexahedra) _hexahedra->init(myVTKData.vtkLevel);
+ if(_pyramids) _pyramids->init(myVTKData.vtkLevel);
+
+ if(_triangles) _triangles->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+ if(_quadrangles) _quadrangles->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+ if(_tetrahedra) _tetrahedra->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+ if(_prisms) _prisms->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+ if(_hexahedra) _hexahedra->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+ if(_pyramids) _pyramids->addInViewForVTK(step, _inData, myVTKData, writeVTK, buildStaticData);
+
+ Msg::StatusBar(true,"Done writing VTK data");
+}
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index adf026f952..80a3399a92 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -11,12 +11,32 @@
#include <vector>
#include <cstdlib>
#include <algorithm>
+#include <sys/stat.h>
+#include <assert.h>
+#include <fstream>
+#include <stdio.h>
+#include <string>
+#include <sstream>
#include "fullMatrix.h"
+
+#define __STDC_FORMAT_MACROS
+#include <inttypes.h>
+typedef std::vector<int> vectInt;
+
class PViewData;
class PViewDataList;
class GMSH_PostPlugin;
+// For old compilers that do not support yet std::to_string()
+template <class T>
+std::string ToString(const T& val)
+{
+ std::stringstream stream;
+ stream << val;
+ return stream.str();
+}
+
class adaptiveVertex {
public:
float x, y, z; //!< parametric coordinates
@@ -38,6 +58,23 @@ class adaptiveVertex {
}
};
+template <class T>
+class nodMap {
+public:
+ std::vector<int> mapping;
+
+public:
+ void cleanMapping()
+ {
+ mapping.clear();
+ }
+ ~nodMap()
+ {
+ cleanMapping();
+ }
+ int getSize() {return (int) mapping.size();}
+};
+
class adaptivePoint {
public:
bool visible;
@@ -352,23 +389,210 @@ class PCoords {
class PValues{
public:
- double v[9];
+ short int sizev; //acceptable values: 1 (scalar), 3 (vector), 9 (tensor)
+ double *v;
+ PValues(const PValues& obj) {
+ sizev = obj.sizev;
+ v = new double[sizev];
+ for(int i=0;i<sizev;i++) {
+ v[i] = obj.v[i];
+ }
+ }
+ PValues(int size)
+ {
+ sizev = size;
+ v = new double[sizev];
+ for(int i=0;i<sizev;i++) {
+ v[i] = 0.0;
+ }
+ }
PValues(double vx)
{
+ sizev = 1;
+ v = new double[sizev];
v[0] = vx;
}
PValues(double vx, double vy, double vz)
{
+ sizev = 3;
+ v = new double[sizev];
v[0] = vx; v[1] = vy; v[2] = vz;
}
PValues(double vxx, double vxy, double vxz,
double vyx, double vyy, double vyz,
double vzx, double vzy, double vzz)
{
+ sizev = 9;
+ v = new double[sizev];
v[0] = vxx; v[1] = vxy; v[2] = vxz;
v[3] = vyx; v[4] = vyy; v[5] = vyz;
v[6] = vzx; v[7] = vzy; v[8] = vzz;
}
+ ~PValues() {
+ delete[] v;
+ }
+ void operator = (const PValues& obj) {
+ // Assume PValues object has already been generated
+ // and v allocated when the operator = is called
+ if(sizev != obj.sizev) Msg::Error("In PValues overlodaing operator: size mistmatch %d %d",sizev);
+ for(int i=0;i<sizev;i++) {
+ v[i] = obj.v[i];
+ }
+ }
+};
+
+class globalVTKData {
+ public:
+
+ static std::vector<vectInt> vtkGlobalConnectivity; // conectivity (vector of vector)
+ static std::vector<int> vtkGlobalCellType; // topology
+ static std::vector<PCoords> vtkGlobalCoords; // coordinates
+ static std::vector<PValues> vtkGlobalValues; // nodal values (either scalar or vector)
+
+ globalVTKData();
+
+ static void clearGlobalConnectivity() {
+ for(std::vector<vectInt>::iterator it = vtkGlobalConnectivity.begin();it != vtkGlobalConnectivity.end(); ++it) {
+ it->clear();
+ }
+ vtkGlobalConnectivity.clear();
+ }
+
+ static void clearGlobalCellType() {
+ vtkGlobalCellType.clear();
+ }
+
+ static void clearGlobalCoords() {
+ vtkGlobalCoords.clear();
+ }
+
+ static void clearGlobalValues() {
+ vtkGlobalValues.clear();
+ }
+
+ static void clearGlobalData() {
+ clearGlobalConnectivity();
+ clearGlobalCellType();
+ clearGlobalCoords();
+ clearGlobalValues();
+ }
+
+ ~globalVTKData() {
+ clearGlobalData();
+ }
+
+};
+
+class VTKData {
+ public:
+ // Data container to write output files readable for ParaView
+ // vtk legacy and vtu for now
+ std::string vtkFieldName;
+ std::string vtkFileName;
+ std::string vtkFormat;
+ std::string vtkDirName;
+
+ int vtkStep;
+ int vtkLevel;
+ int vtkNumComp;
+ double vtkTol;
+ int vtkNpart;
+
+ bool vtkIsBinary;
+ int vtkUseDefaultName;
+ int minElmPerPart, maxElmPerPart, numPartMinElm, numPartMaxElm;
+
+ // File variables
+ FILE *vtkFile;
+ FILE *vtkFileCoord;
+ FILE *vtkFileConnect;
+ FILE *vtkFileCellOffset;
+ FILE *vtkFileCellType;
+ FILE *vtkFileNodVal;
+ int vtkCountFile;
+
+ int vtkTotNumElmLev0;
+ int vtkCountTotElmLev0;
+ int vtkCountTotNod;
+ int vtkCountTotElm;
+ int vtkCountCoord;
+ int vtkCountTotNodConnect;
+ int vtkCountTotVal;
+ int vtkCountCellOffset; //used only for ascii output
+ int vtkCountCellType; //used only for ascii output
+
+ std::vector<vectInt> vtkLocalConnectivity; // conectivity (vector of vector)
+ std::vector<int> vtkLocalCellType; // topology
+ std::vector<PCoords> vtkLocalCoords; // coordinates
+ std::vector<PValues> vtkLocalValues; // nodal values (either scalar or vector)
+
+
+public:
+ VTKData(std::string fieldName="unknown", int numComp = -1, int step = -1, int level = -1, double tol=0.0,
+ std::string filename="unknown", int useDefaultName = 1, int npart = -1, bool isBinary = true) {
+
+ vtkIsBinary = isBinary; // choice: true, false
+ vtkFormat = std::string("vtu"); // choice: vtk (VTK legacy), vtu (XML appended)
+
+ vtkFieldName = fieldName;
+ vtkFileName = filename;
+ vtkUseDefaultName = useDefaultName;
+ vtkNumComp = numComp;
+ vtkStep = step;
+ vtkLevel = level;
+ vtkTol = tol;
+ vtkNpart = npart;
+
+ vtkCountFile = 0;
+ vtkTotNumElmLev0 = 0;
+ vtkCountTotElmLev0 = 0;
+ vtkCountTotNod = 0;
+ vtkCountTotElm = 0;
+ vtkCountCoord = 0;
+ vtkCountTotNodConnect = 0;
+ vtkCountTotVal = 0;
+ vtkCountCellOffset = 0; //used only for ascii output
+ vtkCountCellType = 0;
+ }
+
+ void clearLocalData()
+ {
+ for(std::vector<vectInt>::iterator it = vtkLocalConnectivity.begin();it != vtkLocalConnectivity.end(); ++it) {
+ it->clear();
+ }
+ vtkLocalConnectivity.clear();
+ vtkLocalCellType.clear();
+ vtkLocalCoords.clear();
+ vtkLocalValues.clear();
+ }
+
+ ~VTKData()
+ {
+ clearLocalData();
+ }
+
+ void incrementTotNod(int increment) {vtkCountTotNod+=increment;}
+ void incrementTotElm(int increment) {vtkCountTotElm+=increment;}
+ void incrementTotElmLev0(int increment) {vtkCountTotElmLev0+=increment;}
+
+ bool isLittleEndian();
+ void SwapArrayByteOrder(void* array, int nbytes, int nItems); // used only for VTK
+ int getPVCellType(int numEdges);
+// void writeParaViewData();
+ void writeVTKElmData();
+ void initVTKFile();
+ void finalizeVTKFile();
+ void setFileDistribution() {
+ int tmpmod = vtkTotNumElmLev0 % vtkNpart;
+ minElmPerPart = (vtkTotNumElmLev0-tmpmod)/vtkNpart;
+ numPartMinElm = vtkNpart - tmpmod;
+
+ if(tmpmod == 0 ) maxElmPerPart = minElmPerPart;
+ else maxElmPerPart = minElmPerPart+1;
+ numPartMaxElm = tmpmod;
+ assert(vtkTotNumElmLev0 == minElmPerPart*numPartMinElm+maxElmPerPart*numPartMaxElm);
+
+ }
};
template <class T>
@@ -393,6 +617,26 @@ class adaptiveElements {
// switch to true on-the-fly local refinement in drawPost())
void addInView(double tol, int step, PViewData *in, PViewDataList *out,
GMSH_PostPlugin *plug=0);
+
+ // michel.rasquin@cenaero.be:
+ // Routines for
+ // - export of adapted views to pvtu file format for parallel visualization with paraview,
+ // - and/or generation of VTK data structure for ParaView plugin.
+
+ // Clone of adapt for VTK output files
+ void adaptForVTK(double tol, int numComp,
+ std::vector<PCoords> &coords, std::vector<PValues> &values);
+
+ // Clone of addInView for VTK output files
+ void addInViewForVTK(int step, PViewData *in, VTKData &myVTKData,
+ bool writeVtk=true, bool buildStaticData=false);
+
+ int countElmLev0(int step, PViewData *in);
+
+ // Build a mapping between all the nodes of the refined element
+ // and the node of the canonical refined element in order to
+ // generate a connectivity related to the canonical element
+ void buildMapping(nodMap<T> &myNodMap, double tol, int &numNodInsert);
};
class adaptiveData {
@@ -409,12 +653,32 @@ class adaptiveData {
adaptiveElements<adaptiveHexahedron> *_hexahedra;
adaptiveElements<adaptivePrism> *_prisms;
adaptiveElements<adaptivePyramid> *_pyramids;
+
+ // When set to true, this builds a global VTK data structure (connectivity, coords, etc) for the adaptive views.
+ // This can be very memory consuming for high adaptation levels. Use with caution.
+ // Usefull when GMSH is used as an external library to provide for instance a GMSH reader in a ParaView plugin.
+ // By default, set to false in the constructor.
+ bool buildStaticData;
+
+ // This variable helps limit memory consumption (no global data structure) when GMSH is requested to
+ // write the data structure of adapted view under pvtu format
+ // In this case, one adapted element is considered at a time so that it can generate
+ // billions of adapted elements on a single core, as long as disk space allows it.
+ // This variable is set to true by default in the constructor.
+ bool writeVTK;
+
public:
static double timerInit, timerAdapt;
- adaptiveData(PViewData *data);
+ adaptiveData(PViewData *data, bool outDataInit=true);
~adaptiveData();
PViewData *getData(){ return (PViewData*)_outData; }
void changeResolution(int step, int level, double tol, GMSH_PostPlugin *plug=0);
+ int countTotElmLev0(int step, PViewData *in);
+ void changeResolutionForVTK(int step, int level, double tol, int npart = 1, bool isBinary = true,
+ const std::string &guifileName = "unknown", int useDefaultName = 1);
+ void upBuildStaticData(bool newValue) { buildStaticData = newValue; }
+ void upWriteVTK(bool newValue) { writeVTK = newValue; }
+
};
#endif
--
GitLab