io.cpp 16.3 KB
Newer Older
1
2
3
4
5
6
/*!
 * \file io.cpp
 * \brief Set of utility io functions for rtns
 * \author S. Schotthoefer, J. Wolters et al.
 */

7
#include "common/io.h"
8
9
10
11
#include "common/config.h"
#include "common/mesh.h"
#include "common/typedef.h"

12
#include "toolboxes/errormessages.h"
jannick.wolters's avatar
jannick.wolters committed
13
#include "toolboxes/textprocessingtoolbox.h"
14

15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#include <iostream>

#include <mpi.h>
#include <omp.h>

#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
//#include <vtkPointDataToCellData.h>
#include <vtkQuad.h>
#include <vtkSmartPointer.h>
#include <vtkTriangle.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridWriter.h>

#include <Python.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>

using vtkPointsSP                 = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridSP       = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkTriangleSP               = vtkSmartPointer<vtkTriangle>;
using vtkCellArraySP              = vtkSmartPointer<vtkCellArray>;
using vtkDoubleArraySP            = vtkSmartPointer<vtkDoubleArray>;
using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
using vtkCellDataToPointDataSP    = vtkSmartPointer<vtkCellDataToPointData>;
// using vtkPointDataToCellDataSP    = vtkSmartPointer<vtkPointDataToCellData>;

45
void ExportVTK( const std::string fileName,
46
47
                const std::vector<std::vector<std::vector<double>>>& outputFields,
                const std::vector<std::vector<std::string>>& outputFieldNames,
48
                const Mesh* mesh ) {
Jannick Wolters's avatar
Jannick Wolters committed
49
50
51
52
53
54
55
56
57
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if( rank == 0 ) {
        unsigned dim             = mesh->GetDim();
        unsigned numCells        = mesh->GetNumCells();
        unsigned numNodes        = mesh->GetNumNodes();
        auto nodes               = mesh->GetNodes();
        auto cells               = mesh->GetCells();
        unsigned numNodesPerCell = mesh->GetNumNodesPerCell();
58

Jannick Wolters's avatar
Jannick Wolters committed
59
60
        auto writer                 = vtkUnstructuredGridWriterSP::New();
        std::string fileNameWithExt = fileName;
jannick.wolters's avatar
jannick.wolters committed
61
        if( !TextProcessingToolbox::StringEndsWith( fileNameWithExt, ".vtk" ) ) {
Jannick Wolters's avatar
Jannick Wolters committed
62
            fileNameWithExt.append( ".vtk" );
63
        }
64
        writer->SetFileName( fileNameWithExt.c_str() );
Jannick Wolters's avatar
Jannick Wolters committed
65
66
67
68
69
70
71
72
73
74
75
76
77
        auto grid = vtkUnstructuredGridSP::New();
        auto pts  = vtkPointsSP::New();
        pts->SetDataTypeToDouble();
        pts->SetNumberOfPoints( static_cast<int>( numNodes ) );
        unsigned nodeID = 0;
        for( const auto& node : nodes ) {
            if( dim == 2 ) {
                pts->SetPoint( nodeID++, node[0], node[1], 0.0 );
            }
            else if( dim == 3 ) {
                pts->SetPoint( nodeID++, node[0], node[1], node[2] );
            }
            else {
78
                ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
79
80
            }
        }
Jannick Wolters's avatar
Jannick Wolters committed
81
82
83
84
85
86
        vtkCellArraySP cellArray = vtkCellArraySP::New();
        for( unsigned i = 0; i < numCells; ++i ) {
            if( numNodesPerCell == 3 ) {
                auto tri = vtkTriangleSP::New();
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
                    tri->GetPointIds()->SetId( j, cells[i][j] );
87
                }
Jannick Wolters's avatar
Jannick Wolters committed
88
89
90
91
92
93
                cellArray->InsertNextCell( tri );
            }
            if( numNodesPerCell == 4 ) {
                auto quad = vtkQuad::New();
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
                    quad->GetPointIds()->SetId( j, cells[i][j] );
94
                }
Jannick Wolters's avatar
Jannick Wolters committed
95
96
                cellArray->InsertNextCell( quad );
            }
Jannick Wolters's avatar
Jannick Wolters committed
97
        }
98
        if( numNodesPerCell == 3 ) {
Jannick Wolters's avatar
Jannick Wolters committed
99
            grid->SetCells( VTK_TRIANGLE, cellArray );
100
        }
Jannick Wolters's avatar
Jannick Wolters committed
101
102
        else if( numNodesPerCell == 4 ) {
            grid->SetCells( VTK_QUAD, cellArray );
103
104
        }

105
106
        // Write the output
        for( unsigned idx_group = 0; idx_group < outputFields.size(); idx_group++ ) {
107

108
            for( unsigned idx_field = 0; idx_field < outputFields[idx_group].size(); idx_field++ ) {    // Loop over all output fields
109
110

                auto cellData = vtkDoubleArraySP::New();
111
                cellData->SetName( outputFieldNames[idx_group][idx_field].c_str() );
112
113

                for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
114
                    cellData->InsertNextValue( outputFields[idx_group][idx_field][idx_cell] );
115
116
117
                }

                grid->GetCellData()->AddArray( cellData );
Jannick Wolters's avatar
Jannick Wolters committed
118
            }
119
120
        }

Jannick Wolters's avatar
Jannick Wolters committed
121
122
        grid->SetPoints( pts );
        grid->Squeeze();
123

Jannick Wolters's avatar
Jannick Wolters committed
124
125
126
127
        auto converter = vtkCellDataToPointDataSP::New();
        converter->AddInputDataObject( grid );
        converter->PassCellDataOn();
        converter->Update();
128

Jannick Wolters's avatar
Jannick Wolters committed
129
        auto conv_grid = converter->GetOutput();
130

Jannick Wolters's avatar
Jannick Wolters committed
131
        writer->SetInputData( conv_grid );
132

Jannick Wolters's avatar
Jannick Wolters committed
133
        writer->Write();
134
135
136

        auto log = spdlog::get( "event" );
        log->info( "Result successfully exported to '{0}'!", fileNameWithExt );
Jannick Wolters's avatar
Jannick Wolters committed
137
138
    }
    MPI_Barrier( MPI_COMM_WORLD );
139
140
}

141
Mesh* LoadSU2MeshFromFile( const Config* settings ) {
142
143
144
145
146
147
148
    auto log = spdlog::get( "event" );

    unsigned dim;
    std::vector<Vector> nodes;
    std::vector<std::vector<unsigned>> cells;
    std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries;

149
150
    if( !std::filesystem::exists( settings->GetMeshFile() ) )
        ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
151
152
153
154
155
156
    std::ifstream ifs( settings->GetMeshFile(), std::ios::in );
    std::string line;

    if( ifs.is_open() ) {
        while( getline( ifs, line ) ) {
            if( line.find( "NDIME", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
157
                dim = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
158
159
160
161
162
163
164
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NPOIN", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
165
                unsigned numPoints = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
                nodes.resize( numPoints, Vector( dim, 0.0 ) );
                for( unsigned i = 0; i < numPoints; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned id = 0;
                    Vector coords( dim, 0.0 );
                    while( !ss.eof() ) {
                        for( unsigned d = 0; d < dim; ++d ) {
                            ss >> coords[d];
                        }
                        ss >> id;
                    }

                    nodes[id] = Vector( coords );
                }
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NMARK", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
189
                unsigned numBCs = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
190
191
192
193
194
195
196
197
198
199
200
201
202
                boundaries.resize( numBCs );
                for( unsigned i = 0; i < numBCs; ++i ) {
                    std::string markerTag;
                    BOUNDARY_TYPE btype;
                    std::vector<unsigned> bnodes;
                    for( unsigned k = 0; k < 2; ++k ) {
                        getline( ifs, line );
                        if( line.find( "MARKER_TAG", 0 ) != std::string::npos ) {
                            markerTag    = line.substr( line.find( "=" ) + 1 );
                            auto end_pos = std::remove_if( markerTag.begin(), markerTag.end(), isspace );
                            markerTag.erase( end_pos, markerTag.end() );
                            btype = settings->GetBoundaryType( markerTag );
                            if( btype == BOUNDARY_TYPE::INVALID ) {
203
204
                                std::string errorMsg = std::string( "Invalid Boundary at marker \"" + markerTag + "\"." );
                                ErrorMessages::Error( errorMsg, CURRENT_FUNCTION );
205
206
207
                            }
                        }
                        else if( line.find( "MARKER_ELEMS", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
208
                            unsigned numMarkerElements = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
                            for( unsigned j = 0; j < numMarkerElements; ++j ) {
                                getline( ifs, line );
                                std::stringstream ss;
                                ss << line;
                                unsigned type = 0, id = 0;
                                while( !ss.eof() ) {
                                    ss >> type;
                                    for( unsigned d = 0; d < dim; ++d ) {
                                        ss >> id;
                                        bnodes.push_back( id );
                                    }
                                }
                            }
                        }
                        else {
224
                            ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
225
226
                        }
                    }
Jannick Wolters's avatar
Jannick Wolters committed
227
228
229
                    std::sort( bnodes.begin(), bnodes.end() );
                    auto last = std::unique( bnodes.begin(), bnodes.end() );
                    bnodes.erase( last, bnodes.end() );
230
231
232
233
234
235
236
237
238
239
                    boundaries[i] = std::make_pair( btype, bnodes );
                }
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        std::vector<unsigned> numNodesPerCell;
        while( getline( ifs, line ) ) {
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
240
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
                numNodesPerCell.resize( numCells, 0u );
                for( unsigned i = 0; i < numCells; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned type = 0;
                    ss >> type;
                    switch( type ) {
                        case 5: {
                            numNodesPerCell[i] = 3;
                            break;
                        }
                        case 9: {
                            numNodesPerCell[i] = 4;
                            break;
                        }
                        default: {
258
                            ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
259
260
261
262
263
264
265
266
                        }
                    }
                }
                break;
            }
        }
        bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
        if( mixedElementMesh ) {
267
            ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
268
269
270
271
272
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
273
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
                cells.resize( numCells, std::vector<unsigned>( numNodesPerCell[0], 0u ) );
                for( unsigned i = 0; i < numCells; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned type = 0, id = 0;
                    std::vector<unsigned> buffer( numNodesPerCell[0], 0u );
                    while( !ss.eof() ) {
                        ss >> type;
                        for( unsigned d = 0; d < numNodesPerCell[i]; ++d ) {
                            ss >> buffer[d];
                        }
                        ss >> id;
                    }
                    std::copy( buffer.begin(), buffer.end(), cells[id].begin() );
                }
                break;
            }
        }
    }
    else {
295
        ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
296
297
298
299
300
    }
    ifs.close();
    return new Mesh( nodes, cells, boundaries );
}

301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
std::string ParseArguments( int argc, char* argv[] ) {
    std::string inputFile;
    std::string usage_help = "\n"
                             "Usage: " +
                             std::string( argv[0] ) + " inputfile\n";

    if( argc < 2 ) {
        std::cout << usage_help;
        exit( EXIT_FAILURE );
    }
    for( int i = 1; i < argc; i++ ) {
        std::string arg = argv[i];
        if( arg == "-h" ) {
            std::cout << usage_help;
            exit( EXIT_SUCCESS );
        }
        else {
            inputFile = std::string( argv[i] );
            std::ifstream f( inputFile );
            if( !f.is_open() ) {
321
                ErrorMessages::Error( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
322
323
324
325
326
327
328
            }
        }
    }
    return inputFile;
}

void PrintLogHeader( std::string inputFile ) {
Jannick Wolters's avatar
Jannick Wolters committed
329
    int nprocs, rank;
330
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
Jannick Wolters's avatar
Jannick Wolters committed
331
332
333
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if( rank == 0 ) {
        auto log = spdlog::get( "event" );
334

Jannick Wolters's avatar
Jannick Wolters committed
335
336
337
338
339
340
341
342
343
344
345
346
347
        log->info( "RTSN" );
        log->info( "================================================================" );
        log->info( "Git commit :\t{0}", GIT_HASH );
        log->info( "Config file:\t{0}", inputFile );
        log->info( "MPI Threads:\t{0}", nprocs );
        log->info( "OMP Threads:\t{0}", omp_get_max_threads() );
        log->info( "================================================================" );
        // print file content while omitting comments
        std::ifstream ifs( inputFile );
        if( ifs.is_open() ) {
            std::string line;
            while( !ifs.eof() ) {
                std::getline( ifs, line );
348
                if( line[0] != '%' ) log->info( " {0}", line );
Jannick Wolters's avatar
Jannick Wolters committed
349
            }
350
        }
Jannick Wolters's avatar
Jannick Wolters committed
351
352
        log->info( "================================================================" );
        log->info( "" );
353
    }
Jannick Wolters's avatar
Jannick Wolters committed
354
    MPI_Barrier( MPI_COMM_WORLD );
355
}
356
357
358
359
360
361
362
363
364
365
366
367

Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ) {
    auto log = spdlog::get( "event" );

    if( !std::filesystem::exists( imageName ) ) {
        ErrorMessages::Error( "Can not open image '" + imageName + "'!", CURRENT_FUNCTION );
    }
    std::filesystem::path outDir( std::filesystem::path( SU2Filename ).parent_path() );
    if( !std::filesystem::exists( outDir ) ) {
        ErrorMessages::Error( "Output directory '" + outDir.string() + "' does not exists!", CURRENT_FUNCTION );
    }

jannick.wolters's avatar
jannick.wolters committed
368
369
    std::string pyPath = RTSN_PYTHON_PATH;

jannick.wolters's avatar
jannick.wolters committed
370
    if( !Py_IsInitialized() ) {
jannick.wolters's avatar
jannick.wolters committed
371
372
373
374
375
        Py_InitializeEx( 0 );
        if( !Py_IsInitialized() ) {
            ErrorMessages::Error( "Python init failed!", CURRENT_FUNCTION );
        }
        PyRun_SimpleString( ( "import sys\nsys.path.append('" + pyPath + "')" ).c_str() );
jannick.wolters's avatar
jannick.wolters committed
376
    }
jannick.wolters's avatar
jannick.wolters committed
377

378
379
380
381
382
383
    PyObject *pArgs, *pReturn, *pModule, *pFunc;
    PyArrayObject* np_ret;

    auto image = PyUnicode_FromString( imageName.c_str() );
    auto su2   = PyUnicode_FromString( SU2Filename.c_str() );

jannick.wolters's avatar
jannick.wolters committed
384
385
    std::string moduleName = "mesh_from_image";
    pModule                = PyImport_ImportModule( moduleName.c_str() );
386
    if( !pModule ) {
jannick.wolters's avatar
jannick.wolters committed
387
388
        PyErr_Print();
        ErrorMessages::Error( "'" + moduleName + "' can not be imported!", CURRENT_FUNCTION );
389
390
391
392
    }

    pFunc = PyObject_GetAttrString( pModule, "generate" );
    if( !pFunc || !PyCallable_Check( pFunc ) ) {
jannick.wolters's avatar
jannick.wolters committed
393
394
395
        PyErr_Print();
        Py_DecRef( pModule );
        Py_DecRef( pFunc );
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
        ErrorMessages::Error( "'generate' is null or not callable!", CURRENT_FUNCTION );
    }

    pArgs = PyTuple_New( 2 );
    PyTuple_SetItem( pArgs, 0, reinterpret_cast<PyObject*>( image ) );
    PyTuple_SetItem( pArgs, 1, reinterpret_cast<PyObject*>( su2 ) );
    pReturn = PyObject_CallObject( pFunc, pArgs );
    np_ret  = reinterpret_cast<PyArrayObject*>( pReturn );

    size_t m{ static_cast<size_t>( PyArray_SHAPE( np_ret )[0] ) };
    size_t n{ static_cast<size_t>( PyArray_SHAPE( np_ret )[1] ) };
    double* c_out = reinterpret_cast<double*>( PyArray_DATA( np_ret ) );

    Matrix gsImage( m, n, c_out );

    // Finalizing
jannick.wolters's avatar
jannick.wolters committed
412
413
414
    Py_DecRef( pFunc );
    Py_DecRef( pModule );
    Py_DECREF( np_ret );
415

416
    return gsImage.transpose();
417
}