io.cpp 16.1 KB
Newer Older
1
#include "common/io.h"
2
3
4
5
#include "common/config.h"
#include "common/mesh.h"
#include "common/typedef.h"

6
#include "toolboxes/errormessages.h"
jannick.wolters's avatar
jannick.wolters committed
7
#include "toolboxes/textprocessingtoolbox.h"
8

9
10
11
12
13
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
#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>;

39
void ExportVTK( const std::string fileName,
40
41
                const std::vector<std::vector<std::vector<double>>>& outputFields,
                const std::vector<std::vector<std::string>>& outputFieldNames,
42
                const Mesh* mesh ) {
Jannick Wolters's avatar
Jannick Wolters committed
43
44
45
46
47
48
49
50
51
    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();
52

Jannick Wolters's avatar
Jannick Wolters committed
53
54
        auto writer                 = vtkUnstructuredGridWriterSP::New();
        std::string fileNameWithExt = fileName;
jannick.wolters's avatar
jannick.wolters committed
55
        if( !TextProcessingToolbox::StringEndsWith( fileNameWithExt, ".vtk" ) ) {
Jannick Wolters's avatar
Jannick Wolters committed
56
            fileNameWithExt.append( ".vtk" );
57
        }
58
        writer->SetFileName( fileNameWithExt.c_str() );
Jannick Wolters's avatar
Jannick Wolters committed
59
60
61
62
63
64
65
66
67
68
69
70
71
        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 {
72
                ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
73
74
            }
        }
Jannick Wolters's avatar
Jannick Wolters committed
75
76
77
78
79
80
        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] );
81
                }
Jannick Wolters's avatar
Jannick Wolters committed
82
83
84
85
86
87
                cellArray->InsertNextCell( tri );
            }
            if( numNodesPerCell == 4 ) {
                auto quad = vtkQuad::New();
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
                    quad->GetPointIds()->SetId( j, cells[i][j] );
88
                }
Jannick Wolters's avatar
Jannick Wolters committed
89
90
                cellArray->InsertNextCell( quad );
            }
Jannick Wolters's avatar
Jannick Wolters committed
91
        }
92
        if( numNodesPerCell == 3 ) {
Jannick Wolters's avatar
Jannick Wolters committed
93
            grid->SetCells( VTK_TRIANGLE, cellArray );
94
        }
Jannick Wolters's avatar
Jannick Wolters committed
95
96
        else if( numNodesPerCell == 4 ) {
            grid->SetCells( VTK_QUAD, cellArray );
97
98
        }

99
100
        // Write the output
        for( unsigned idx_group = 0; idx_group < outputFields.size(); idx_group++ ) {
101

102
            for( unsigned idx_field = 0; idx_field < outputFields[idx_group].size(); idx_field++ ) {    // Loop over all output fields
103
104

                auto cellData = vtkDoubleArraySP::New();
105
                cellData->SetName( outputFieldNames[idx_group][idx_field].c_str() );
106
107

                for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
108
                    cellData->InsertNextValue( outputFields[idx_group][idx_field][idx_cell] );
109
110
111
                }

                grid->GetCellData()->AddArray( cellData );
Jannick Wolters's avatar
Jannick Wolters committed
112
            }
113
114
        }

Jannick Wolters's avatar
Jannick Wolters committed
115
116
        grid->SetPoints( pts );
        grid->Squeeze();
117

Jannick Wolters's avatar
Jannick Wolters committed
118
119
120
121
        auto converter = vtkCellDataToPointDataSP::New();
        converter->AddInputDataObject( grid );
        converter->PassCellDataOn();
        converter->Update();
122

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

Jannick Wolters's avatar
Jannick Wolters committed
125
        writer->SetInputData( conv_grid );
126

Jannick Wolters's avatar
Jannick Wolters committed
127
        writer->Write();
128
129
130

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

135
Mesh* LoadSU2MeshFromFile( const Config* settings ) {
136
137
138
139
140
141
142
    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;

143
144
    if( !std::filesystem::exists( settings->GetMeshFile() ) )
        ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
145
146
147
148
149
150
    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
151
                dim = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
152
153
154
155
156
157
158
                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
159
                unsigned numPoints = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
                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
183
                unsigned numBCs = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
184
185
186
187
188
189
190
191
192
193
194
195
196
                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 ) {
197
198
                                std::string errorMsg = std::string( "Invalid Boundary at marker \"" + markerTag + "\"." );
                                ErrorMessages::Error( errorMsg, CURRENT_FUNCTION );
199
200
201
                            }
                        }
                        else if( line.find( "MARKER_ELEMS", 0 ) != std::string::npos ) {
jannick.wolters's avatar
jannick.wolters committed
202
                            unsigned numMarkerElements = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
                            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 {
218
                            ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
219
220
                        }
                    }
Jannick Wolters's avatar
Jannick Wolters committed
221
222
223
                    std::sort( bnodes.begin(), bnodes.end() );
                    auto last = std::unique( bnodes.begin(), bnodes.end() );
                    bnodes.erase( last, bnodes.end() );
224
225
226
227
228
229
230
231
232
233
                    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
234
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
                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: {
252
                            ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
253
254
255
256
257
258
259
260
                        }
                    }
                }
                break;
            }
        }
        bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
        if( mixedElementMesh ) {
261
            ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
262
263
264
265
266
        }
        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
267
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
                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 {
289
        ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
290
291
292
293
294
    }
    ifs.close();
    return new Mesh( nodes, cells, boundaries );
}

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
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() ) {
315
                ErrorMessages::Error( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
316
317
318
319
320
321
322
            }
        }
    }
    return inputFile;
}

void PrintLogHeader( std::string inputFile ) {
Jannick Wolters's avatar
Jannick Wolters committed
323
    int nprocs, rank;
324
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
Jannick Wolters's avatar
Jannick Wolters committed
325
326
327
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if( rank == 0 ) {
        auto log = spdlog::get( "event" );
328

Jannick Wolters's avatar
Jannick Wolters committed
329
330
331
332
333
334
335
336
337
338
339
340
341
        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 );
342
                if( line[0] != '%' ) log->info( " {0}", line );
Jannick Wolters's avatar
Jannick Wolters committed
343
            }
344
        }
Jannick Wolters's avatar
Jannick Wolters committed
345
346
        log->info( "================================================================" );
        log->info( "" );
347
    }
Jannick Wolters's avatar
Jannick Wolters committed
348
    MPI_Barrier( MPI_COMM_WORLD );
349
}
350
351
352
353
354
355
356
357
358
359
360
361

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
362
363
    std::string pyPath = RTSN_PYTHON_PATH;

jannick.wolters's avatar
jannick.wolters committed
364
    if( !Py_IsInitialized() ) {
jannick.wolters's avatar
jannick.wolters committed
365
366
367
368
369
        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
370
    }
jannick.wolters's avatar
jannick.wolters committed
371

372
373
374
375
376
377
    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
378
379
    std::string moduleName = "mesh_from_image";
    pModule                = PyImport_ImportModule( moduleName.c_str() );
380
    if( !pModule ) {
jannick.wolters's avatar
jannick.wolters committed
381
382
        PyErr_Print();
        ErrorMessages::Error( "'" + moduleName + "' can not be imported!", CURRENT_FUNCTION );
383
384
385
386
    }

    pFunc = PyObject_GetAttrString( pModule, "generate" );
    if( !pFunc || !PyCallable_Check( pFunc ) ) {
jannick.wolters's avatar
jannick.wolters committed
387
388
389
        PyErr_Print();
        Py_DecRef( pModule );
        Py_DecRef( pFunc );
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
        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
406
407
408
    Py_DecRef( pFunc );
    Py_DecRef( pModule );
    Py_DECREF( np_ret );
409

410
    return gsImage.transpose();
411
}