MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Private Attributes | Static Private Attributes | List of all members
VtkIo< nDim, SysEqn > Class Template Reference

#include <iovtk.h>

Collaboration diagram for VtkIo< nDim, SysEqn >:
[legend]

Public Types

using SolverCell = FvCell
 

Public Member Functions

 VtkIo (FvCartesianSolverXD< nDim, SysEqn > *)
 
 ~VtkIo ()
 
SysEqn sysEqn () const
 
SysEqn & sysEqn ()
 
MInt solverCount () const
 
MInt noDomains () const
 
MInt domainId () const
 
MPI_Comm mpiComm () const
 
MInt writeVtuArrayParallel (MPI_File &, void *, MPI_Offset, MPI_Offset, MPI_Offset)
 
void insertDataHeader (char *data, uint64_t &memsize, uint64_t &memsizeGlobal, uint64_t &offset)
 The zeroth domain stores the header for an uncompressed vtu file appended data array specifying its numer of bytes. More...
 
void updateDataOffsets (uint64_t memsize, uint64_t &memsizeGlobal, uint64_t &offset, uint64_t &oldMemsizeGlobal)
 Recomputes global offsets and data size given the local memsize. More...
 
MInt estimateMemorySizeSolverwise (uint64_t, ScratchSpace< uint64_t > &, uint64_t)
 
MBool initializeVtkXmlOutput (const MChar *, MString, MInt, MBool, MBool)
 
void writeVtkXmlOutput (const MChar *fileName)
 Copy of parallel single-file VTU output using MPI I/O. More...
 
template<typename uint_t = uint32_t>
void writeVtuOutputParallel (const MString fileName, const MString geomFileName, const MInt noSolverSpecificVars=0, const MFloatScratchSpace &solverSpecificVars=MFloatScratchSpace(1, "writeVtuOutputParallel", "defaultParameter1"), const MInt noUserVars=0, const MFloatScratchSpace &userVars=MFloatScratchSpace(1, "writeVtuOutputParallel", "defaultParameter2"))
 
template<typename T >
uint64_t vtuAssembleFaceStream (std::vector< T > *&, std::vector< T > *&, uint64_t &, ScratchSpace< MInt > &, MInt &, const MBool)
 

Private Attributes

FvCartesianSolverXD< nDim, SysEqn > * m_solver = nullptr
 
Collector< FvBndryCell< nDim, SysEqn > > * m_bndryCells = nullptr
 
Collector< PointBasedCell< nDim > > * m_extractedCells = nullptr
 
Collector< CartesianGridPoint< nDim > > * m_gridPoints = nullptr
 
List< MInt > * m_bndryCellIds = nullptr
 
MInt m_noSolver
 

Static Private Attributes

static constexpr MInt m_noDirs = 2 * nDim
 

Detailed Description

template<MInt nDim, class SysEqn>
class VtkIo< nDim, SysEqn >

Definition at line 24 of file iovtk.h.

Member Typedef Documentation

◆ SolverCell

template<MInt nDim, class SysEqn >
using VtkIo< nDim, SysEqn >::SolverCell = FvCell

Definition at line 37 of file iovtk.h.

Constructor & Destructor Documentation

◆ VtkIo()

template<MInt nDim, class SysEqn >
VtkIo< nDim, SysEqn >::VtkIo ( FvCartesianSolverXD< nDim, SysEqn > *  solver)

Definition at line 15 of file iovtk.cpp.

15 : m_solver(solver) {
16 DEBUG("VtkIo::VtkIo: entry", MAIA_DEBUG_ALLOCATION);
17 DEBUG("VtkIo::VtkIo: return", MAIA_DEBUG_ALLOCATION);
18}
VtkIo::m_solver
FvCartesianSolverXD< nDim, SysEqn > * m_solver
Definition: iovtk.h:26

◆ ~VtkIo()

template<MInt nDim, class SysEqn >
VtkIo< nDim, SysEqn >::~VtkIo

Definition at line 22 of file iovtk.cpp.

22 {
23 DEBUG("VtkIo::~VtkIo: entry", MAIA_DEBUG_ALLOCATION);
24 DEBUG("VtkIo::~VtkIo: return", MAIA_DEBUG_ALLOCATION);
25}

Member Function Documentation

◆ domainId()

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::domainId ( ) const
inline

Definition at line 47 of file iovtk.h.

47{ return m_solver->domainId(); }
Solver::domainId
virtual MInt domainId() const
Return the domainId (rank)
Definition: solver.h:383

◆ estimateMemorySizeSolverwise()

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::estimateMemorySizeSolverwise ( uint64_t  noElements,
ScratchSpace< uint64_t > &  noElementsPerDomain,
uint64_t  factor 
)

Definition at line 1142 of file iovtk.cpp.

1144 {
1145 ASSERT(noElements == noElementsPerDomain[domainId()], "");
1146 MInt ret = noElements;
1147 if(domainId() == 0) ret += (MInt)((sizeof(uint64_t) + factor - 1) / factor);
1148 return mMax(1, ret);
1149}
VtkIo::domainId
MInt domainId() const
Definition: iovtk.h:47
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
MInt
int32_t MInt
Definition: maiatypes.h:62

◆ initializeVtkXmlOutput()

template<MInt nDim, class SysEqn >
MBool VtkIo< nDim, SysEqn >::initializeVtkXmlOutput ( const MChar *  fileName,
MString  outputDir_,
MInt  solutionInterval,
MBool  restart,
MBool  firstUseInitializeVtkXmlOutput 
)

Definition at line 29 of file iovtk.cpp.

30 {
31 TRACE();
32 if(domainId() != 0) { // only root process executes the following code
33 return true;
34 }
35
36 m_log << "VtkIo::initializeVtkXmlOutput" << endl;
37
38 MString outputFileName = "QOUT";
39 MString pvd = ".pvd";
40
41 MString pvdPath = outputDir_ + outputFileName + pvd; // e.g. out/QOUT.pvd
42 MString pvdTmpPath = pvdPath + ".tmp";
43 MString pvdAllPath = outputDir_ + outputFileName + "_all" + pvd;
44 MString buPath = outputDir_ + outputFileName + "_BU" + pvd;
45 MString vtuPath = MString("./") + outputFileName + "_";
46
47 if(firstUseInitializeVtkXmlOutput) {
48 if(fileExists(pvdPath.c_str())) {
49 rename(pvdPath.c_str(), buPath.c_str());
50 }
51 if(fileExists(pvdTmpPath.c_str())) {
52 remove(pvdTmpPath.c_str());
53 }
54 ofstream ofile(pvdTmpPath.c_str(), ios_base::out | ios_base::trunc);
55 if(ofile.is_open() && ofile.good()) {
56 ofile << R"(<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">)" << endl;
57 ofile << "<Collection>" << endl;
58 if(restart) {
59 for(MInt t = 0; t <= globalTimeStep; t += solutionInterval) {
60 for(MInt p = 0; p < noDomains(); p++) {
61 stringstream tmp;
62 tmp << outputDir_ << fileName << "_D" << p << "_00" << t << ".vtu";
63 if(fileExists((tmp.str()).c_str())) {
64 ofile << "<DataSet part=\"" << p << "\" timestep=\"" << t << "\" file=\"" << vtuPath << t << "_D" << p
65 << ".vtu\"/>\n";
66 }
67 }
68 }
69 }
70 ofile.close();
71 ofile.clear();
72 } else {
73 cerr << "Error opening file " << pvdTmpPath << endl;
74 }
75 }
76 if(firstUseInitializeVtkXmlOutput) {
77 ofstream ofile(pvdAllPath.c_str(), ios_base::out | ios_base::trunc);
78 if(ofile.is_open() && ofile.good()) {
79 ofile << R"(<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">)" << endl;
80 ofile << "<Collection>" << endl;
81 MInt t = 0;
82 MBool found = true;
83 while(found) {
84 for(MInt p = 0; p < noDomains(); p++) {
85 stringstream tmp;
86 tmp << outputDir_ << fileName << "_D" << p << "_00" << t << ".vtu";
87 if(fileExists(tmp.str())) {
88 ofile << "<DataSet part=\"" << p << "\" timestep=\"" << t << "\" file=\"" << vtuPath << t << "_D" << p
89 << ".vtu\"/>\n";
90 } else {
91 found = false;
92 }
93 }
94 t += solutionInterval;
95 }
96 ofile << "</Collection>" << endl;
97 ofile << "</VTKFile>" << endl;
98 ofile.close();
99 ofile.clear();
100 } else {
101 cerr << "Error opening file " << pvdAllPath << endl;
102 }
103 }
104 ofstream ofile(pvdTmpPath.c_str(), ios_base::out | ios_base::app);
105 if(ofile.is_open() && ofile.good()) {
106 for(MInt p = 0; p < noDomains(); p++) {
107 ofile << "<DataSet part=\"" << p << "\" timestep=\"" << globalTimeStep << "\" file=\"" << vtuPath
108 << globalTimeStep << "_D" << p << ".vtu\"/>\n";
109 }
110 ofile.close();
111 ofile.clear();
112 } else {
113 cerr << "Error opening file " << pvdTmpPath << endl;
114 }
115
116 if(fileExists(pvdPath)) {
117 remove(pvdPath.c_str());
118 }
119 copyFile(pvdTmpPath.c_str(), pvdPath);
120 ofstream ofile2(pvdPath.c_str(), ios_base::out | ios_base::app);
121 if(ofile2.is_open() && ofile2.good()) {
122 ofile2 << "</Collection>" << endl;
123 ofile2 << "</VTKFile>" << endl;
124 ofile2.close();
125 } else {
126 cerr << "Error opening file " << pvdPath << endl;
127 }
128
129 return true;
130}
VtkIo::noDomains
MInt noDomains() const
Definition: iovtk.h:46
copyFile
MInt copyFile(const MString &fromName, const MString &toName)
Copies file fromName to file toName.
Definition: functions.cpp:83
fileExists
MBool fileExists(const MString &fileName)
Returns true if the file fileName exists, false otherwise.
Definition: functions.cpp:73
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
MBool
bool MBool
Definition: maiatypes.h:58
maia::fc::cell::p
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Definition: fccellproperties.h:32

◆ insertDataHeader()

template<MInt nDim, class SysEqn >
void VtkIo< nDim, SysEqn >::insertDataHeader ( char *  data,
uint64_t &  memsize,
uint64_t &  memsizeGlobal,
uint64_t &  offset 
)
Author
Lennart Scheiders

Definition at line 1157 of file iovtk.cpp.

1157 {
1158 if(domainId() == 0) {
1159 if(memsize > 0) {
1160 memmove(&data[sizeof(uint64_t)], &data[0], memsize);
1161 }
1162 uint64_t dataSize = memsizeGlobal;
1163 memcpy(&data[0], reinterpret_cast<void*>(&dataSize), sizeof(uint64_t));
1164 memsize += sizeof(uint64_t);
1165 } else {
1166 offset += sizeof(uint64_t);
1167 }
1168 memsizeGlobal += sizeof(uint64_t);
1169}

◆ mpiComm()

template<MInt nDim, class SysEqn >
MPI_Comm VtkIo< nDim, SysEqn >::mpiComm ( ) const
inline

Definition at line 48 of file iovtk.h.

48{ return m_solver->mpiComm(); }
Solver::mpiComm
MPI_Comm mpiComm() const
Return the MPI communicator used by this solver.
Definition: solver.h:380

◆ noDomains()

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::noDomains ( ) const
inline

Definition at line 46 of file iovtk.h.

46{ return m_solver->noDomains(); }
Solver::noDomains
virtual MInt noDomains() const
Definition: solver.h:387

◆ solverCount()

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::solverCount ( ) const
inline

Definition at line 45 of file iovtk.h.

45{ return m_noSolver; };
VtkIo::m_noSolver
MInt m_noSolver
Definition: iovtk.h:32

◆ sysEqn() [1/2]

template<MInt nDim, class SysEqn >
SysEqn & VtkIo< nDim, SysEqn >::sysEqn ( )
inline

Definition at line 44 of file iovtk.h.

44{ return m_solver->m_sysEqn; };

◆ sysEqn() [2/2]

template<MInt nDim, class SysEqn >
SysEqn VtkIo< nDim, SysEqn >::sysEqn ( ) const
inline

Definition at line 43 of file iovtk.h.

43{ return m_solver->m_sysEqn; };

◆ updateDataOffsets()

template<MInt nDim, class SysEqn >
void VtkIo< nDim, SysEqn >::updateDataOffsets ( uint64_t  memsize,
uint64_t &  memsizeGlobal,
uint64_t &  offset,
uint64_t &  oldMemsizeGlobal 
)
Author
Lennart Scheiders

Definition at line 1177 of file iovtk.cpp.

1178 {
1179 ScratchSpace<uint64_t> tmpCntPerDomain(noDomains(), AT_, "tmpCntPerDomain");
1180 MPI_Allgather(&memsize, 1, MPI_UINT64_T, &tmpCntPerDomain[0], 1, MPI_UINT64_T, mpiComm(), AT_, "memsize",
1181 "tmpCntPerDomain[0]");
1182 oldMemsizeGlobal += memsizeGlobal;
1183 memsizeGlobal = 0;
1184 offset = 0;
1185 for(MInt d = 0; d < noDomains(); d++) {
1186 memsizeGlobal += tmpCntPerDomain[d];
1187 }
1188 for(MInt d = 0; d < domainId(); d++) {
1189 offset += tmpCntPerDomain[d];
1190 }
1191}
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
VtkIo::mpiComm
MPI_Comm mpiComm() const
Definition: iovtk.h:48
MPI_Allgather
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
Definition: mpioverride.cpp:1274

◆ vtuAssembleFaceStream()

template<MInt nDim, class SysEqn >
template<typename T >
uint64_t VtkIo< nDim, SysEqn >::vtuAssembleFaceStream ( std::vector< T > *&  facestream,
std::vector< T > *&  conn,
uint64_t &  facesSize,
ScratchSpace< MInt > &  polyMap,
MInt &  polyCnt,
const MBool  getInternalPolyhedra 
)

Definition at line 1195 of file iovtk.cpp.

1197 {
1198 TRACE();
1199
1200 static constexpr MInt pointCode[6][4] = {{0, 2, 6, 4}, {1, 3, 7, 5}, {0, 1, 5, 4},
1201 {2, 3, 7, 6}, {0, 1, 3, 2}, {4, 5, 7, 6}};
1202 static constexpr MInt edgeCode[6][4] = {{0, 10, 4, 8}, {1, 11, 5, 9}, {2, 9, 6, 8},
1203 {3, 11, 7, 10}, {2, 1, 3, 0}, {6, 5, 7, 4}};
1204 static constexpr MInt faceStencil[2][12] = {{0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1},
1205 {4, 4, 4, 4, 5, 5, 5, 5, 2, 2, 3, 3}};
1206 // static constexpr MInt revDir[6] = {1,0,3,2,5,4};
1207 const MFloat eps0 = 0.0001;
1208 const MFloat eps = 1e-12;
1209 const MInt noNodes = IPOW2(nDim);
1210 const MInt noBndryCells = m_solver->m_bndryCells->size();
1211 MIntScratchSpace cellMap(m_solver->a_noCells(), FUN_, "cellMap");
1212 // MIntScratchSpace polyMap(m_solver->m_extractedCells->size(), FUN_, "polyMap");
1213 MIntScratchSpace revMap(m_solver->a_noCells(), FUN_, "revMap");
1214 MIntScratchSpace edgeCheck(50, 50, FUN_, "edgeCheck");
1215 MIntScratchSpace nghbrList(200, AT_, "nghbrList");
1216 MIntScratchSpace layerId(200, AT_, "layerId");
1217 vector<MInt> noInternalNodes(noBndryCells, 0);
1218 // noInternalNodes.resize( noBndryCells );
1219 // facestream.resize( noBndryCells );
1220 // conn.resize( noBndryCells );
1221 // facestream = new vector<T> [ mMax(1,noBndryCells) ];
1222 // conn = new vector<T> [ mMax(1,noBndryCells) ];
1223 uint64_t connSize = 0;
1224 facesSize = 0;
1225 cellMap.fill(-1);
1226 polyMap.fill(-1);
1227 revMap.fill(-1);
1228
1229 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1230 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1231 cellMap[cellId] = c;
1232 }
1233 polyCnt = 0;
1234 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1235 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1236 MInt bndryId = m_solver->a_bndryId(cellId);
1237 if(bndryId > -1) {
1238 polyMap[c] = polyCnt;
1239 revMap[polyCnt] = cellId;
1240 polyCnt++;
1241 }
1242 }
1243 if(getInternalPolyhedra) {
1244 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1245 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1246 MInt bndryId = m_solver->a_bndryId(cellId);
1247 if(bndryId > -1) continue;
1248 MBool poly = false;
1249 MBool polyEdge = false;
1250 for(MInt i = 0; i < m_noDirs; i++) {
1251 if(m_solver->checkNeighborActive(cellId, i) && m_solver->a_hasNeighbor(cellId, i) > 0) {
1252 if(m_solver->c_noChildren(m_solver->c_neighborId(cellId, i)) > 0) {
1253 poly = true;
1254 }
1255 }
1256 for(MInt j = 0; j < m_noDirs; j++) {
1257 if((i / 2) == (j / 2)) continue;
1258 MInt nb0 = (m_solver->checkNeighborActive(cellId, i) && m_solver->a_hasNeighbor(cellId, i) > 0)
1259 ? m_solver->c_neighborId(cellId, i)
1260 : cellId;
1261 MInt nb1 = (m_solver->checkNeighborActive(cellId, j) && m_solver->a_hasNeighbor(cellId, j) > 0)
1262 ? m_solver->c_neighborId(cellId, j)
1263 : cellId;
1264 MInt nb01 = (m_solver->checkNeighborActive(nb0, j) && m_solver->a_hasNeighbor(nb0, j) > 0)
1265 ? m_solver->c_neighborId(nb0, j)
1266 : cellId;
1267 MInt nb10 = (m_solver->checkNeighborActive(nb1, i) && m_solver->a_hasNeighbor(nb1, i) > 0)
1268 ? m_solver->c_neighborId(nb1, i)
1269 : cellId;
1270 polyEdge = polyEdge || (m_solver->c_noChildren(nb0) > 0) || (m_solver->c_noChildren(nb1) > 0)
1271 || (m_solver->c_noChildren(nb01) > 0) || (m_solver->c_noChildren(nb10) > 0);
1272 }
1273 }
1274 if(poly || polyEdge) {
1275 polyMap[c] = polyCnt;
1276 revMap[polyCnt] = cellId;
1277 polyCnt++;
1278 }
1279 }
1280 }
1281
1282 if(facestream != nullptr) {
1283 delete[] facestream;
1284 }
1285 if(conn != nullptr) {
1286 delete[] conn;
1287 }
1288 facestream = new vector<T>[mMax(1, polyCnt)];
1289 conn = new vector<T>[mMax(1, polyCnt)];
1290 // find duplicate points
1291 /*
1292 for ( MInt c = 0; c < m_solver->m_extractedCells->size(); c++ ) {
1293 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1294 //MInt bndryId = a_bndryId( cellId );
1295 //if ( bndryId < 0 ) {
1296 //const MInt counter = getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
1297 const MInt counter = getAdjacentLeafCells_d2( cellId, 1, nghbrList, layerId );
1298 for ( MInt n = 0; n < counter; n++ ) {
1299 MInt nghbrCell = nghbrList[n];
1300 if ( nghbrCell < 0 ) continue;
1301 if ( a_isHalo(nghbrCell) ) continue;
1302 MInt ncell = cellMap[nghbrCell];
1303 if( ncell < 0 ) continue;
1304 //if (a_bndryId( nghbrCell )> -1 ) continue;
1305 for ( MInt p = 0; p < noNodes; p++ ) {
1306 for ( MInt q = 0; q < noNodes; q++ ) {
1307 MInt gridPointId = m_solver->m_extractedCells->a[c].m_pointIds[p];
1308 MInt nghbrGridPointId = m_solver->m_extractedCells->a[ncell].m_pointIds[q];
1309 if ( gridPointId == nghbrGridPointId ) continue;
1310 if ( nghbrGridPointId < 0 ) cerr << "negative grid point id " << endl;
1311 MFloat delta = sqrt( POW2( m_solver->m_gridPoints->a[ nghbrGridPointId ].m_coordinates[0] -
1312 m_solver->m_gridPoints->a[ gridPointId ].m_coordinates[0] )
1313 + POW2( m_solver->m_gridPoints->a[ nghbrGridPointId ].m_coordinates[1] -
1314 m_solver->m_gridPoints->a[ gridPointId ].m_coordinates[1] )
1315 + POW2( m_solver->m_gridPoints->a[ nghbrGridPointId ].m_coordinates[2] -
1316 m_solver->m_gridPoints->a[ gridPointId ].m_coordinates[2] ) ); if ( delta < eps0*c_cellLengthAtCell(cellId) ) {
1317 //cerr << "found another match 0" << endl;
1318 MInt pid = ( gridPointId < nghbrGridPointId ) ? gridPointId : nghbrGridPointId;
1319 MInt pid2 = ( gridPointId < nghbrGridPointId ) ? nghbrGridPointId : gridPointId;
1320 MInt cid = ( gridPointId < nghbrGridPointId ) ? ncell : c;
1321 m_solver->m_extractedCells->a[cid].m_pointIds[p] = pid;
1322 if ( m_solver->m_gridPoints->a[ pid ].m_noAdjacentCells < noNodes ) {
1323 m_solver->m_gridPoints->a[ pid ].m_cellIds[ m_solver->m_gridPoints->a[ pid ].m_noAdjacentCells ] =
1324 m_solver->m_extractedCells->a[cid].m_cellId; m_solver->m_gridPoints->a[ pid ].m_noAdjacentCells++;
1325 }
1326 for ( MInt k = 0; k < m_solver->m_gridPoints->a[ pid2 ].m_noAdjacentCells; k++ ) {
1327 if ( m_solver->m_gridPoints->a[ pid2 ].m_cellIds[ k ] == m_solver->m_extractedCells->a[cid].m_cellId ) {
1328 m_solver->m_gridPoints->a[ pid2 ].m_cellIds[ k ] = -1;
1329 m_solver->m_gridPoints->a[ pid2 ].m_noAdjacentCells--;
1330 if ( k < m_solver->m_gridPoints->a[ pid2 ].m_noAdjacentCells )
1331 m_solver->m_gridPoints->a[ pid2 ].m_cellIds[ k ] = m_solver->m_gridPoints->a[ pid2 ].m_cellIds[
1332 m_solver->m_gridPoints->a[ pid2
1333 ].m_noAdjacentCells ];
1334 }
1335 }
1336 if ( m_solver->m_gridPoints->a[ pid2 ].m_noAdjacentCells == 0 ) {
1337 for ( MInt j = 0; j < nDim; j++ ) {
1338 m_solver->m_gridPoints->a[ pid2 ].m_coordinates[ j ] = std::numeric_limits<float>::max();
1339 }
1340 }
1341 }
1342 }
1343 }
1344 }
1345 //}
1346 }
1347 */
1348
1349
1350 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1351 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1352 MInt bndryId = m_solver->a_bndryId(cellId);
1353 /*MBool poly = false;
1354 MBool polyEdge = false;
1355 if ( noEntries > polyCnt0 ) {
1356 for ( MInt i = 0; i < m_noDirs; i++ ) {
1357 if ( a_hasNeighbor( cellId , i ) > 0 ) {
1358 if ( c_noChildren( c_neighborId( cellId , i ) ) > 0 ) {
1359 poly = true;
1360 }
1361 }
1362 for ( MInt j = 0; j < m_noDirs; j++ ) {
1363 if( (i/2) == (j/2) ) continue;
1364 MInt nb0 = ( a_hasNeighbor( cellId , i ) > 0 ) ? c_neighborId( cellId , i ) : cellId;
1365 MInt nb1 = ( a_hasNeighbor( cellId , j ) > 0 ) ? c_neighborId( cellId , j ) : cellId;
1366 MInt nb01 = ( a_hasNeighbor( nb0 , j ) > 0 ) ? c_neighborId( nb0 , j ) : cellId;
1367 MInt nb10 = ( a_hasNeighbor( nb1 , i ) > 0 ) ? c_neighborId( nb1 , i ) : cellId;
1368 polyEdge = polyEdge || ( c_noChildren( nb0 ) > 0 ) || ( c_noChildren( nb1 ) > 0 ) || ( c_noChildren( nb01 ) >
1369 0 ) || ( c_noChildren( nb10 ) > 0 );
1370 }
1371 }
1372 }
1373 if ( bndryId < 0 ) {
1374 if ( !poly && !polyEdge ) {
1375 connSize += (uint64_t)noNodes;
1376 facesSize++;
1377 }
1378 continue;
1379 }*/
1380 MInt polyId = polyMap[c];
1381 if(polyId < 0) {
1382 connSize += (uint64_t)noNodes;
1383 facesSize++;
1384 continue;
1385 }
1386 facestream[polyId].clear();
1387 conn[polyId].clear();
1388 for(MInt i = 0; i < noNodes; i++) {
1389 const MInt gridPointId = m_solver->m_extractedCells->a[c].m_pointIds[i];
1390 if(gridPointId < 0) continue;
1391 if(bndryId > -1 && m_solver->gridPointIsInside(c, i)) {
1392 /*for ( MInt k = 0; k < m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells; k++ ) {
1393 if ( m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ k ] == m_solver->m_extractedCells->a[c].m_cellId ) {
1394 m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ k ] = -1;
1395 m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells--;
1396 if ( k < m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells )
1397 m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ k ] = m_solver->m_gridPoints->a[ gridPointId
1398 ].m_cellIds[ m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells ];
1399 //break;
1400 }
1401 }
1402 if ( m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells == 0 ) {
1403 for ( MInt j = 0; j < nDim; j++ ) {
1404 m_solver->m_gridPoints->a[ gridPointId ].m_coordinates[ j ] = std::numeric_limits<float>::max();
1405 }
1406 }*/
1407 m_solver->m_extractedCells->a[c].m_pointIds[i] = -1;
1408
1409 /*
1410 for ( MInt k = 0; k < m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells; k++ ) {
1411 MInt ncellId = m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ k ];
1412 if ( ncellId < 0 ) continue;
1413 MInt ncell = cellMap[ncellId];
1414 if ( ncell < 0 ) continue;
1415 for ( MInt p = 0; p < noNodes; p++ ) {
1416 const MInt gridPointId2 = m_solver->m_extractedCells->a[ ncell ].m_pointIds[ p ];
1417 if ( gridPointId2 < 0 ) continue;
1418 if ( gridPointId == gridPointId2 ) {
1419 m_solver->m_extractedCells->a[ ncell ].m_pointIds[ p ] = -1;
1420 }
1421 }
1422 }
1423 for ( MInt p = 0; p < noNodes; p++ ) {
1424 m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ p ] = -1;
1425 }
1426 m_solver->m_extractedCells->a[ c ].m_pointIds[ i ] = -1;
1427 m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells = 0;
1428 for ( MInt j = 0; j < nDim; j++ ) {
1429 m_solver->m_gridPoints->a[ gridPointId ].m_coordinates[ j ] = std::numeric_limits<float>::max();
1430 }
1431 */
1432 } else {
1433 conn[polyId].push_back(gridPointId);
1434 /*MBool exist = false;
1435 for ( MInt k = 0; k < m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells; k++ ){
1436 if ( m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[k] == cellId ) exist = true;
1437 }
1438 if ( !exist && m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells < noNodes ) {
1439 m_solver->m_gridPoints->a[ gridPointId ].m_cellIds[ m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells
1440 ] = cellId; m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells++;
1441 }*/
1442 }
1443 }
1444 }
1445
1446
1447 const MInt oldNoPoints = m_solver->m_gridPoints->size();
1448 MIntScratchSpace pointRefs(oldNoPoints, FUN_, "pointRefs");
1449 pointRefs.fill(0);
1450 for(MInt gridPointId = 0; gridPointId < oldNoPoints; gridPointId++) {
1451 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
1452 for(MInt p = 0; p < noNodes; p++) {
1453 m_solver->m_gridPoints->a[gridPointId].m_cellIds[p] = -1;
1454 }
1455 }
1456 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1457 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1458 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
1459 ? m_solver->getAssociatedInternalCell(cellId)
1460 : cellId;
1461 // MInt bndryId = a_bndryId( cellId );
1462 MInt pc = polyMap[c];
1463 // MInt noPoints = ( bndryId > -1 && pc > -1 ) ? conn[pc].size() : noNodes;
1464 MInt noPoints = (pc > -1) ? conn[pc].size() : noNodes;
1465 for(MInt q = 0; q < noPoints; q++) {
1466 // for ( MInt q = 0; q < noNodes; q++ ) {
1467 // MInt gridPointId = ( bndryId > -1 && pc > -1 ) ? conn[pc][q] : m_solver->m_extractedCells->a[c].m_pointIds[q];
1468 MInt gridPointId = (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[c].m_pointIds[q];
1469 if(gridPointId >= oldNoPoints) cerr << "point out of range" << endl;
1470 // MInt gridPointId = m_solver->m_extractedCells->a[c].m_pointIds[q];
1471 if(gridPointId > -1) {
1472 MBool found = false;
1473 for(MInt n = 0; n < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; n++) {
1474 if(m_solver->m_gridPoints->a[gridPointId].m_cellIds[n] == rootId) found = true;
1475 }
1476 if(!found && (m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes)) {
1477 m_solver->m_gridPoints->a[gridPointId].m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] =
1478 cellId;
1479 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
1480 }
1481 pointRefs[gridPointId]++;
1482 }
1483 }
1484 }
1485 for(MInt gridPointId = 0; gridPointId < oldNoPoints; gridPointId++) {
1486 if(pointRefs[gridPointId] > noNodes) cerr << "too many point refs " << pointRefs[gridPointId] << endl;
1487 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells != pointRefs[gridPointId]) {
1488 cerr << "strange0 " << m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells << " " << pointRefs[gridPointId]
1489 << endl;
1490 }
1491 // if ( pointRefs[gridPointId] == 0 ) {
1492 // m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells = 0;
1493 //}
1494 }
1495 /*
1496 MInt inact = 0;
1497 for(MInt gridPointId = 0; gridPointId < oldNoPoints; gridPointId++ ) {
1498 if ( pointRefs[gridPointId] == 0 ) {
1499 if ( m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells > 0 ) {
1500 cerr <<"strange " << m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells <<endl;
1501 inact++;
1502 }
1503
1504 }
1505 }
1506 cerr << "inact " << inact << endl;*/
1507 MInt lastId = m_solver->m_gridPoints->size() - 1;
1508 for(MInt gridPointId = 0; gridPointId < lastId; gridPointId++) {
1509 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells > 0) continue;
1510 // if ( m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells > 0 && pointRefs[gridPointId] > 0 ) continue;
1511 // if ( pointRefs[gridPointId] > 0 ) continue;
1512 lastId = m_solver->m_gridPoints->size() - 1;
1513 while(m_solver->m_gridPoints->a[lastId].m_noAdjacentCells == 0) {
1514 // while ( pointRefs[gridPointId] == 0 ) {
1515 m_solver->m_gridPoints->resetSize(lastId);
1516 lastId = m_solver->m_gridPoints->size() - 1;
1517 }
1518 if(gridPointId >= m_solver->m_gridPoints->size()) {
1519 break;
1520 } else {
1521 lastId = m_solver->m_gridPoints->size() - 1;
1522 // if ( m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells > 0 ) {
1523 // cerr << "Warning: grid point still active." << endl;
1524 // continue;
1525 //}
1526 if(gridPointId < lastId) {
1527 for(MInt k = 0; k < m_solver->m_gridPoints->a[lastId].m_noAdjacentCells; k++) {
1528 MInt cellId = m_solver->m_gridPoints->a[lastId].m_cellIds[k];
1529 if(cellId < 0) continue;
1530 if(cellMap[cellId] > -1) {
1531 for(MInt i = 0; i < noNodes; i++) {
1532 if(m_solver->m_extractedCells->a[cellMap[cellId]].m_pointIds[i] == lastId) {
1533 m_solver->m_extractedCells->a[cellMap[cellId]].m_pointIds[i] = gridPointId;
1534 }
1535 }
1536 // MInt bndryId = a_bndryId( cellId );
1537 // if ( bndryId > -1 ) {
1538 MInt polyId = polyMap[cellMap[cellId]];
1539 if(polyId > -1) {
1540 replace(conn[polyId].begin(), conn[polyId].end(), lastId, gridPointId);
1541 }
1542 //}
1543 }
1544 if(m_solver->a_hasProperty(cellId, SolverCell::IsSplitCell)) {
1545 auto it = find(m_solver->m_splitCells.begin(), m_solver->m_splitCells.end(), cellId);
1546 if(it == m_solver->m_splitCells.end()) mTerm(1, AT_, "split cells inconsistency.");
1547 const MInt pos = distance(m_solver->m_splitCells.begin(), it);
1548 ASSERT(m_solver->m_splitCells[pos] == cellId, "");
1549 for(MUint c = 0; c < m_solver->m_splitChilds[pos].size(); c++) {
1550 MInt splitChilId = m_solver->m_splitChilds[pos][c];
1551 if(cellMap[splitChilId] > -1) {
1552 for(MInt i = 0; i < noNodes; i++) {
1553 if(m_solver->m_extractedCells->a[cellMap[splitChilId]].m_pointIds[i] == lastId) {
1554 m_solver->m_extractedCells->a[cellMap[splitChilId]].m_pointIds[i] = gridPointId;
1555 }
1556 }
1557 MInt polyId = polyMap[cellMap[splitChilId]];
1558 if(polyId > -1) {
1559 replace(conn[polyId].begin(), conn[polyId].end(), lastId, gridPointId);
1560 }
1561 }
1562 }
1563 }
1564 }
1565 for(MInt k = 0; k < m_solver->m_gridPoints->a[lastId].m_noAdjacentCells; k++) {
1566 m_solver->m_gridPoints->a[gridPointId].m_cellIds[k] = m_solver->m_gridPoints->a[lastId].m_cellIds[k];
1567 m_solver->m_gridPoints->a[lastId].m_cellIds[k] = -1;
1568 }
1569 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = m_solver->m_gridPoints->a[lastId].m_noAdjacentCells;
1570 m_solver->m_gridPoints->a[lastId].m_noAdjacentCells = 0;
1571 pointRefs[gridPointId] = pointRefs[lastId];
1572 pointRefs[lastId] = 0;
1573 for(MInt j = 0; j < nDim; j++) {
1574 m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] = m_solver->m_gridPoints->a[lastId].m_coordinates[j];
1575 m_solver->m_gridPoints->a[lastId].m_coordinates[j] = std::numeric_limits<float>::max();
1576 }
1577 }
1578 m_solver->m_gridPoints->resetSize(lastId);
1579 }
1580 }
1581
1582
1583 multimap<MFloat, MInt> sortedCPs;
1584 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1585 for(MInt i = 0; i < noNodes; i++) {
1586 if(m_solver->m_extractedCells->a[c].m_pointIds[i] >= m_solver->m_gridPoints->size()) {
1587 cerr << domainId() << ": Error invalid grid point " << c << " " << m_solver->m_extractedCells->a[c].m_cellId
1588 << " " << m_solver->m_extractedCells->a[c].m_pointIds[i] << " " << m_solver->m_gridPoints->size() << endl;
1589 m_solver->m_extractedCells->a[c].m_pointIds[i] = 0;
1590 }
1591 }
1592 MInt bndryId = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
1593 if(bndryId < 0) continue;
1594 MInt polyId = polyMap[c];
1595 for(MUint i = 0; i < conn[polyId].size(); i++) {
1596 if(conn[polyId][i] >= (uint64_t)m_solver->m_gridPoints->size()) {
1597 cerr << domainId() << ": Error invalid boundary grid point " << c << " "
1598 << m_solver->m_extractedCells->a[c].m_cellId << " " << conn[polyId][i] << " "
1599 << m_solver->m_gridPoints->size() << endl;
1600 conn[polyId][i] = 0;
1601 }
1602 }
1603 noInternalNodes[bndryId] = 0;
1604 for(MInt i = 0; i < noNodes; i++) {
1605 if(m_solver->m_extractedCells->a[c].m_pointIds[i] > -1) {
1606 noInternalNodes[bndryId]++;
1607 }
1608 }
1609 }
1610 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1611 MInt bndryId = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
1612 if(bndryId < 0) continue;
1613 MInt polyId = polyMap[c];
1614 MBool outside[8];
1615 for(MInt i = 0; i < noNodes; i++) {
1616 outside[i] = true;
1617 if(m_solver->m_extractedCells->a[c].m_pointIds[i] > -1) {
1618 outside[i] = false;
1619 }
1620 }
1621
1622 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
1623 for(MInt cp = 0; cp < m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints; cp++) {
1624 MInt cutEdge = m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutEdge[cp];
1625 MInt gridPointId = -1;
1626 // check if cut point has been created by a neighbor cell prevoiusly
1627 if(cutEdge > -1) {
1628 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1629 MInt f0 = faceStencil[0][cutEdge];
1630 MInt f1 = faceStencil[1][cutEdge];
1631 MInt nghbrs[4] = {-1, -1, -1, -1};
1632 // MInt nDirs[ 4 ][ 2 ] = {{revDir[f0],f1},{f0,revDir[f1]},{revDir[f0],revDir[f1]},{revDir[f0],revDir[f1]}};
1633 if(!m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild) && m_solver->checkNeighborActive(cellId, f0)
1634 && m_solver->a_hasNeighbor(cellId, f0) > 0) {
1635 nghbrs[0] = m_solver->a_bndryId(m_solver->c_neighborId(cellId, f0));
1636 if(m_solver->checkNeighborActive(m_solver->c_neighborId(cellId, f0), f1)
1637 && m_solver->a_hasNeighbor(m_solver->c_neighborId(cellId, f0), f1) > 0) {
1638 nghbrs[3] = m_solver->a_bndryId(m_solver->c_neighborId(m_solver->c_neighborId(cellId, f0), f1));
1639 }
1640 }
1641 if(!m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild) && m_solver->checkNeighborActive(cellId, f1)
1642 && m_solver->a_hasNeighbor(cellId, f1) > 0) {
1643 nghbrs[1] = m_solver->a_bndryId(m_solver->c_neighborId(cellId, f1));
1644 if(m_solver->checkNeighborActive(m_solver->c_neighborId(cellId, f1), f0)
1645 && m_solver->a_hasNeighbor(m_solver->c_neighborId(cellId, f1), f0) > 0) {
1646 nghbrs[2] = m_solver->a_bndryId(m_solver->c_neighborId(m_solver->c_neighborId(cellId, f1), f0));
1647 }
1648 }
1649 for(MInt nBndryId : nghbrs) {
1650 if(gridPointId > -1) continue;
1651 if(nBndryId > -1) {
1652 if(noInternalNodes[nBndryId] == 0) continue;
1653 if(cellMap[m_solver->m_bndryCells->a[nBndryId].m_cellId] < 0) continue;
1654 MInt npc = polyMap[cellMap[m_solver->m_bndryCells->a[nBndryId].m_cellId]];
1655 if(npc < 0) continue;
1656 MInt noPointsNghbr = conn[npc].size();
1657 for(MInt q = noInternalNodes[nBndryId]; q < noPointsNghbr; q++) {
1658 MInt nghbrGridPointId = conn[npc][q];
1659 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
1660 MFloat delta =
1661 sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0]
1662 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][0])
1663 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1]
1664 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][1])
1665 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2]
1666 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][2]));
1667 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
1668 // cerr << "found another match 0" << endl;
1669 gridPointId = nghbrGridPointId;
1670 break;
1671 }
1672 }
1673 /*
1674 MInt ncnt = 0;
1675 for( MInt nSrfc = 0; nSrfc < m_solver->m_bndryCells->a[ nBndryId ].m_noSrfcs; nSrfc++ ){
1676 for( MInt ncp = 0; ncp < m_solver->m_bndryCells->a[ nBndryId ].m_srfcs[nSrfc]->m_noCutPoints; ncp++ )
1677 { MInt nCutEdge = m_solver->m_bndryCells->a[ nBndryId ].m_srfcs[nSrfc]->m_cutEdge[ ncp ]; if (
1678 m_solver->m_bndryCells->a[ bndryId ].m_srfcs[srfc]->m_bodyId[0] == m_solver->m_bndryCells->a[ nBndryId
1679 ].m_srfcs[nSrfc]->m_bodyId[0] ) { if( nCutEdge > -1 ){ MInt nf0 = faceStencil[0][nCutEdge]; MInt nf1
1680 = faceStencil[1][nCutEdge]; if (((nf0 == nDirs[n][0]) && (nf1 == nDirs[n][1])) ||
1681 ((nf0 == nDirs[n][1]) && (nf1 == nDirs[n][0]))) {
1682 MFloat delta = 0;
1683 for( MInt i = 0; i < nDim; i++ ){
1684 delta += POW2( m_solver->m_bndryCells->a[ nBndryId
1685 ].m_srfcs[nSrfc]->m_cutCoordinates[ncp][i] - m_solver->m_bndryCells->a[ bndryId
1686 ].m_srfcs[srfc]->m_cutCoordinates[cp][i] );
1687 }
1688 if( sqrt(delta) < eps*c_cellLengthAtCell(cellId) ) {
1689 MInt nbCell = m_solver->m_bndryCells->a[ nBndryId ].m_cellId;
1690 if ( cellMap[nbCell] > -1 && polyMap[cellMap[nbCell]] > -1 ) {
1691 if ( noInternalNodes[nBndryId]+ncnt < (signed)conn[ polyMap[cellMap[nbCell]] ].size() ) {
1692 gridPointId = conn[ polyMap[cellMap[nbCell]] ][noInternalNodes[nBndryId]+ncnt];
1693 }
1694 }
1695 }
1696 }
1697 }
1698 }
1699 ncnt++;
1700 }
1701 }*/
1702 }
1703 }
1704 }
1705 if(gridPointId < 0) {
1706 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1707 // const MInt counter = getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
1708 const MInt counter = m_solver->getAdjacentLeafCells_d2(cellId, 1, nghbrList, layerId);
1709 for(MInt n = 0; n < counter; n++) {
1710 MInt nghbrCell = nghbrList[n];
1711 if(nghbrCell < 0) continue;
1712 if(m_solver->a_isHalo(nghbrCell)) continue;
1713 MInt ncell = cellMap[nghbrCell];
1714 if(ncell < 0) continue;
1715 MInt pc = polyMap[ncell];
1716 if(m_solver->a_bndryId(nghbrCell) < 0) continue;
1717 if(pc < 0) continue;
1718 MInt nBndryId = m_solver->a_bndryId(nghbrCell);
1719 MInt noPointsNghbr = conn[pc].size();
1720 for(MInt q = noInternalNodes[nBndryId]; q < noPointsNghbr; q++) {
1721 MInt nghbrGridPointId = conn[pc][q];
1722 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
1723 MFloat delta = sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0]
1724 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][0])
1725 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1]
1726 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][1])
1727 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2]
1728 - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][2]));
1729 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
1730 // cerr << "found another match 0" << endl;
1731 gridPointId = nghbrGridPointId;
1732 break;
1733 }
1734 }
1735 }
1736 }
1737 // if not, create it
1738 if(gridPointId < 0) {
1739 gridPointId = m_solver->m_gridPoints->size();
1740 m_solver->m_gridPoints->append();
1741 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
1742 for(MInt i = 0; i < nDim; i++) {
1743 m_solver->m_gridPoints->a[gridPointId].m_coordinates[i] =
1744 m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][i];
1745 }
1746 }
1747 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells >= noNodes) {
1748 cerr << "nodes " << gridPointId << " " << m_solver->m_extractedCells->a[c].m_cellId << " "
1749 << m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells << "; ";
1750 for(MInt k = 0; k < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; k++)
1751 cerr << m_solver->m_gridPoints->a[gridPointId].m_cellIds[k] << " ";
1752 cerr << endl;
1753 }
1754 ASSERT(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes, "");
1755 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes) {
1756 MBool exist = false;
1757 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1758 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
1759 ? m_solver->getAssociatedInternalCell(cellId)
1760 : cellId;
1761 for(MInt k = 0; k < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; k++) {
1762 if(m_solver->m_gridPoints->a[gridPointId].m_cellIds[k] == rootId) exist = true;
1763 }
1764 if(!exist) {
1765 m_solver->m_gridPoints->a[gridPointId].m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] =
1766 rootId;
1767 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
1768 }
1769 }
1770 conn[polyId].push_back(gridPointId);
1771 }
1772 }
1773
1774 if(m_solver->m_bndryCells->a[bndryId].m_faceVertices.empty()) {
1775 facestream[polyId].push_back(0); // number of faces
1776 MInt pointCntOffset = 1;
1777
1778 MInt cnt = 0;
1779 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
1780 sortedCPs.clear();
1781 MFloat pCoords[3]{};
1782 MFloat vec_a[3]{};
1783 MFloat vec_b[3]{};
1784 MFloat normal[3]{};
1785 for(MInt i = 0; i < nDim; i++) {
1786 normal[i] = m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_normalVector[i];
1787 }
1788 MInt spaceId = 0;
1789 MFloat maxC = fabs(normal[0]);
1790 for(MInt i = 1; i < nDim; i++) {
1791 if(fabs(normal[i]) > maxC) {
1792 maxC = fabs(normal[i]);
1793 spaceId = i;
1794 }
1795 }
1796 MInt spaceId1 = (spaceId + 1) % nDim;
1797 MInt spaceId2 = (spaceId1 + 1) % nDim;
1798 vec_a[spaceId1] = F1;
1799 vec_a[spaceId2] = F1;
1800 vec_a[spaceId] = -(vec_a[spaceId1] * normal[spaceId1] + vec_a[spaceId2] * normal[spaceId2]) / normal[spaceId];
1801 MFloat vecsum = sqrt(POW2(vec_a[0]) + POW2(vec_a[1]) + POW2(vec_a[2]));
1802 for(MInt i = 0; i < nDim; i++) {
1803 vec_a[i] *= 1000.0 / vecsum;
1804 }
1805 vec_b[spaceId] = normal[spaceId1] * vec_a[spaceId2] - normal[spaceId2] * vec_a[spaceId1];
1806 vec_b[spaceId1] = normal[spaceId2] * vec_a[spaceId] - normal[spaceId] * vec_a[spaceId2];
1807 vec_b[spaceId2] = normal[spaceId] * vec_a[spaceId1] - normal[spaceId1] * vec_a[spaceId];
1808 vecsum = sqrt(POW2(vec_b[0]) + POW2(vec_b[1]) + POW2(vec_b[2]));
1809 for(MInt i = 0; i < nDim; i++) {
1810 vec_b[i] *= 1000.0 * vecsum;
1811 }
1812 for(MInt cp = 0; cp < m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints; cp++) {
1813 for(MInt i = 0; i < nDim; i++) {
1814 pCoords[i] = m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutCoordinates[cp][i];
1815 }
1816 MFloat dx = F0;
1817 MFloat dy = F0;
1818 for(MInt i = 0; i < nDim; i++) {
1819 dx += vec_a[i] * (pCoords[i] - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_coordinates[i]);
1820 dy += vec_b[i] * (pCoords[i] - m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_coordinates[i]);
1821 }
1822 sortedCPs.insert(make_pair(atan2(dy, dx), cp));
1823 }
1824 ASSERT((signed)sortedCPs.size() == m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints, "");
1825 facestream[polyId].push_back(0); // number of points for current face
1826 for(auto& sortedCP : sortedCPs) {
1827 facestream[polyId].push_back(conn[polyId][noInternalNodes[bndryId] + cnt + sortedCP.second]);
1828 facestream[polyId][pointCntOffset]++;
1829 }
1830 facestream[polyId][0]++;
1831 pointCntOffset = facestream[polyId].size();
1832 cnt += m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints;
1833 }
1834
1835
1836 for(MInt face = 0; face < m_noDirs; face++) {
1837 // if ( m_solver->m_extractedCells->a[ c ].m_pointIds[ pointCode[face][0] ] < 0 &&
1838 // m_solver->m_extractedCells->a[ c ].m_pointIds[ pointCode[face][1] ] < 0
1839 // && m_solver->m_extractedCells->a[ c ].m_pointIds[ pointCode[face][2] ] < 0 &&
1840 // m_solver->m_extractedCells->a[ c ].m_pointIds[ pointCode[face][3] ] < 0 ) continue;
1841 sortedCPs.clear();
1842 MFloat pCoords[MAX_SPACE_DIMENSIONS] = {F0, F0, F0};
1843 MFloat vec_a[MAX_SPACE_DIMENSIONS] = {F0, F0, F0};
1844 MFloat vec_b[MAX_SPACE_DIMENSIONS] = {F0, F0, F0};
1845 MInt spaceId = face / 2;
1846 MInt spaceId1 = (spaceId + 1) % nDim;
1847 MInt spaceId2 = (spaceId1 + 1) % nDim;
1848 vec_a[spaceId1] = F1;
1849 vec_b[spaceId2] = F1;
1850 MFloat sum = F0;
1851 for(MInt p = 0; p < 4; p++) {
1852 if(!outside[pointCode[face][p]]) {
1853 for(MInt i = 0; i < nDim; i++) {
1854 pCoords[i] += m_solver->m_gridPoints->a[m_solver->m_extractedCells->a[c].m_pointIds[pointCode[face][p]]]
1855 .m_coordinates[i];
1856 }
1857 sum += F1;
1858 }
1859 MInt ccnt = 0;
1860 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
1861 for(MInt cp = 0; cp < m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints; cp++) {
1862 if(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutEdge[cp] == edgeCode[face][p]) {
1863 for(MInt i = 0; i < nDim; i++) {
1864 pCoords[i] +=
1865 m_solver->m_gridPoints->a[conn[polyId][noInternalNodes[bndryId] + ccnt]].m_coordinates[i];
1866 }
1867 sum += F1;
1868 }
1869 ccnt++;
1870 }
1871 }
1872 }
1873 for(MInt i = 0; i < nDim; i++) {
1874 pCoords[i] /= sum;
1875 }
1876
1877 for(MInt p = 0; p < 4; p++) {
1878 if(!outside[pointCode[face][p]]) {
1879 MFloat dx = F0;
1880 MFloat dy = F0;
1881 for(MInt i = 0; i < nDim; i++) {
1882 dx += vec_a[i]
1883 * (m_solver->m_gridPoints->a[m_solver->m_extractedCells->a[c].m_pointIds[pointCode[face][p]]]
1884 .m_coordinates[i]
1885 - pCoords[i]);
1886 dy += vec_b[i]
1887 * (m_solver->m_gridPoints->a[m_solver->m_extractedCells->a[c].m_pointIds[pointCode[face][p]]]
1888 .m_coordinates[i]
1889 - pCoords[i]);
1890 }
1891 sortedCPs.insert(make_pair(atan2(dy, dx), m_solver->m_extractedCells->a[c].m_pointIds[pointCode[face][p]]));
1892 }
1893 MInt ccnt = 0;
1894 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
1895 for(MInt cp = 0; cp < m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_noCutPoints; cp++) {
1896 if(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_cutEdge[cp] == edgeCode[face][p]) {
1897 MFloat dx = F0;
1898 MFloat dy = F0;
1899 for(MInt i = 0; i < nDim; i++) {
1900 dx += vec_a[i]
1901 * (m_solver->m_gridPoints->a[conn[polyId][noInternalNodes[bndryId] + ccnt]].m_coordinates[i]
1902 - pCoords[i]);
1903 dy += vec_b[i]
1904 * (m_solver->m_gridPoints->a[conn[polyId][noInternalNodes[bndryId] + ccnt]].m_coordinates[i]
1905 - pCoords[i]);
1906 }
1907 sortedCPs.insert(
1908 make_pair(atan2(dy, dx), static_cast<MInt>(conn[polyId][noInternalNodes[bndryId] + ccnt])));
1909 }
1910 ccnt++;
1911 }
1912 }
1913 }
1914
1915 if(!sortedCPs.empty()) {
1916 facestream[polyId].push_back(sortedCPs.size()); // number of points for current face
1917 for(auto& sortedCP : sortedCPs) {
1918 facestream[polyId].push_back(sortedCP.second);
1919 }
1920 } else
1921 continue;
1922
1923 /*
1924 facestream[polyId].push_back(0); //number of points for current face
1925 for ( MInt p = 0; p < 4; p++ ) {
1926 if ( !outside[ pointCode[face][p] ] ) {
1927 facestream[polyId].push_back( m_solver->m_extractedCells->a[ c ].m_pointIds[ pointCode[face][p] ] );
1928 facestream[polyId][pointCntOffset]++;
1929 }
1930 MInt ccnt = 0;
1931 for( MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++ ) {
1932 for( MInt cp = 0; cp < m_solver->m_bndryCells->a[ bndryId ].m_srfcs[srfc]->m_noCutPoints; cp++ ) {
1933 if ( m_solver->m_bndryCells->a[ bndryId ].m_srfcs[srfc]->m_cutEdge[cp] == edgeCode[face][p] ) {
1934 facestream[polyId].push_back( conn[polyId][noInternalNodes[bndryId]+ccnt] );
1935 facestream[polyId][pointCntOffset]++;
1936 }
1937 ccnt++;
1938 }
1939 }
1940 }*/
1941
1942
1943 facestream[polyId][0]++;
1944 pointCntOffset = facestream[polyId].size();
1945 }
1946
1947 sort(conn[polyId].begin(), conn[polyId].end(), less<T>());
1948 conn[polyId].erase(unique(conn[polyId].begin(), conn[polyId].end()), conn[polyId].end());
1949 connSize += (uint64_t)conn[polyId].size();
1950 facesSize += (uint64_t)facestream[polyId].size();
1951 }
1952 }
1953
1954 // cells with split faces and other special multi cuts
1955 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
1956 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
1957 MInt bndryId = m_solver->a_bndryId(cellId);
1958 if(bndryId < 0) continue;
1959 MInt polyId = polyMap[c];
1960
1961 if(!m_solver->m_bndryCells->a[bndryId].m_faceVertices.empty()) {
1962 // const MInt gridPointOffset = m_solver->m_gridPoints->size();
1963 facestream[polyId].push_back(0); // number of faces
1964 // econn[polyId].clear();
1965 vector<MInt> tmpPointMap;
1966 tmpPointMap.resize(m_solver->m_bndryCells->a[bndryId].m_faceVertices.size() / 3);
1967 for(MUint v = 0; v < m_solver->m_bndryCells->a[bndryId].m_faceVertices.size() / 3; v++) {
1968 tmpPointMap[v] = -1;
1969 }
1970 for(MUint v = 0; v < m_solver->m_bndryCells->a[bndryId].m_faceVertices.size() / 3; v++) {
1971 MFloat ccoords[3]{};
1972 for(MInt i = 0; i < nDim; i++) {
1973 ccoords[i] = m_solver->m_bndryCells->a[bndryId].m_faceVertices[3 * v + i];
1974 }
1975 MInt gridPointId = -1;
1976 MInt noPoints = conn[polyId].size();
1977 for(MInt q = 0; q < noPoints; q++) {
1978 MInt tmpGridPointId = conn[polyId][q];
1979 if(tmpGridPointId < 0) cerr << "negative grid point id " << endl;
1980 MFloat delta = sqrt(POW2(m_solver->m_gridPoints->a[tmpGridPointId].m_coordinates[0] - ccoords[0])
1981 + POW2(m_solver->m_gridPoints->a[tmpGridPointId].m_coordinates[1] - ccoords[1])
1982 + POW2(m_solver->m_gridPoints->a[tmpGridPointId].m_coordinates[2] - ccoords[2]));
1983 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
1984 gridPointId = tmpGridPointId;
1985 }
1986 }
1987 if(gridPointId < 0) {
1988 // const MInt counter = m_solver->getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
1989 const MInt counter = m_solver->getAdjacentLeafCells_d2(cellId, 1, nghbrList, layerId);
1990 for(MInt n = 0; n < counter; n++) {
1991 MInt nghbrCell = nghbrList[n];
1992 if(nghbrCell < 0) continue;
1993 if(m_solver->a_isHalo(nghbrCell)) continue;
1994 MInt ncell = cellMap[nghbrCell];
1995 if(ncell < 0) continue;
1996 MInt pc = polyMap[ncell];
1997 // if (m_solver->a_bndryId( nghbrCell )> -1 ) continue;
1998 MInt noPointsNghbr = (pc > -1) ? conn[pc].size() : noNodes;
1999 for(MInt q = 0; q < noPointsNghbr; q++) {
2000 MInt nghbrGridPointId = (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[ncell].m_pointIds[q];
2001 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
2002 MFloat delta = sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0] - ccoords[0])
2003 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1] - ccoords[1])
2004 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2] - ccoords[2]));
2005 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
2006 // cerr << "found another match 1" << endl;
2007 gridPointId = nghbrGridPointId;
2008 conn[polyId].push_back(gridPointId);
2009 break;
2010 }
2011 }
2012 }
2013 }
2014 if(gridPointId < 0) {
2015 gridPointId = m_solver->m_gridPoints->size();
2016 m_solver->m_gridPoints->append();
2017 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
2018 for(MInt i = 0; i < nDim; i++) {
2019 m_solver->m_gridPoints->a[gridPointId].m_coordinates[i] =
2020 m_solver->m_bndryCells->a[bndryId].m_faceVertices[3 * v + i];
2021 }
2022 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
2023 ? m_solver->getAssociatedInternalCell(cellId)
2024 : cellId;
2025 m_solver->m_gridPoints->a[gridPointId].m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] =
2026 rootId;
2027 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
2028 conn[polyId].push_back(gridPointId);
2029 }
2030 tmpPointMap[v] = gridPointId;
2031 }
2032 for(MUint f = 0; f < m_solver->m_bndryCells->a[bndryId].m_faceStream.size(); f++) {
2033 facestream[polyId][0]++;
2034 facestream[polyId].push_back(m_solver->m_bndryCells->a[bndryId].m_faceStream[f].size());
2035 for(MUint v = 0; v < m_solver->m_bndryCells->a[bndryId].m_faceStream[f].size(); v++) {
2036 if(tmpPointMap[m_solver->m_bndryCells->a[bndryId].m_faceStream[f][v]] < 0)
2037 cerr << "invalid split point" << endl;
2038 // facestream[polyId].push_back(gridPointOffset+m_solver->m_bndryCells->a[bndryId].m_faceStream[f][v]);
2039 facestream[polyId].push_back(tmpPointMap[m_solver->m_bndryCells->a[bndryId].m_faceStream[f][v]]);
2040 }
2041 }
2042 sort(conn[polyId].begin(), conn[polyId].end(), less<T>());
2043 conn[polyId].erase(unique(conn[polyId].begin(), conn[polyId].end()), conn[polyId].end());
2044 connSize += (uint64_t)conn[polyId].size();
2045 facesSize += (uint64_t)facestream[polyId].size();
2046 }
2047 }
2048
2049
2050 if(getInternalPolyhedra) {
2051 const MInt ltable[6][4] = {{0, 2, 6, 4}, {1, 3, 7, 5}, {0, 1, 5, 4}, {2, 3, 7, 6}, {0, 1, 3, 2}, {4, 5, 7, 6}};
2052 const MInt ltable2[6][2] = {{1, 7}, {0, 6}, {2, 7}, {0, 5}, {4, 7}, {0, 3}};
2053 const MInt ltable3[6][4] = {{4, 3, 5, 2}, {4, 3, 5, 2}, {4, 1, 5, 0}, {4, 1, 5, 0}, {2, 1, 3, 0}, {2, 1, 3, 0}};
2054 const MFloat signStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
2055 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
2056 const MInt dirPlusOne[4] = {1, 2, 3, 0};
2057 const MInt dirMinusOne[4] = {3, 0, 1, 2};
2058 const MInt dirPlusTwo[4] = {2, 3, 0, 1};
2059 const MInt reverseDir[6] = {1, 0, 3, 2, 5, 4};
2060
2061
2062 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
2063 const MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
2064 const MInt bndryId = m_solver->a_bndryId(cellId);
2065 if(bndryId > -1) continue;
2066 MInt polyId = polyMap[c];
2067 if(polyId < 0) continue;
2068 if(polyId >= polyCnt) mTerm(1, AT_, "entries out of range");
2069 facestream[polyId].clear();
2070 /*conn[polyId].clear();
2071 for ( MInt i = 0; i < noNodes; i++ ) {
2072 const MInt gridPointId = m_solver->m_extractedCells->a[ c ].m_pointIds[ i ];
2073 if ( gridPointId < 0 ) continue;
2074 conn[polyId].push_back( gridPointId );
2075 }*/
2076 facestream[polyId].push_back(0); // number of faces
2077 MInt pointCntOffset = 1;
2078 MBool polyFaceSide[6];
2079 for(MInt i = 0; i < m_noDirs; i++) {
2080 MBool polySide = false;
2081 if(m_solver->checkNeighborActive(cellId, i) && m_solver->a_hasNeighbor(cellId, i) > 0) {
2082 if(m_solver->c_noChildren(m_solver->c_neighborId(cellId, i)) > 0) {
2083 polySide = true;
2084 }
2085 }
2086 polyFaceSide[i] = polySide;
2087 if(polySide) {
2088 // store four faces (polygon numbering!)
2089 MInt edgeIds[4] = {-1, -1, -1, -1};
2090 MInt centerId = -1;
2091 MInt nghbrId0 = m_solver->c_childId(m_solver->c_neighborId(cellId, i), ltable2[i][0]);
2092 if(nghbrId0 > -1 && m_solver->a_bndryId(nghbrId0) > -1) nghbrId0 = -1;
2093 MInt nghbrId = (nghbrId0 < 0) ? -1 : cellMap[nghbrId0];
2094
2095 MFloat ccoords[3]{};
2096 for(MInt j = 0; j < nDim; j++) {
2097 ccoords[j] = m_solver->a_coordinate(cellId, j);
2098 }
2099 ccoords[i / 2] +=
2100 (i % 2 == 0) ? -F1B2 * m_solver->c_cellLengthAtCell(cellId) : F1B2 * m_solver->c_cellLengthAtCell(cellId);
2101
2102 MInt gridPointId = -1;
2103 for(MInt e = 0; e < (MInt)conn[polyId].size(); e++) {
2104 MFloat delta = F0;
2105 for(MInt j = 0; j < nDim; j++) {
2106 delta += POW2(ccoords[j] - m_solver->m_gridPoints->a[conn[polyId][e]].m_coordinates[j]);
2107 }
2108 delta = sqrt(delta);
2109 if(delta < eps0 * m_solver->c_cellLengthAtCell(cellId)) {
2110 gridPointId = conn[polyId][e];
2111 break;
2112 }
2113 }
2114 if(gridPointId < 0 && nghbrId > -1) {
2115 gridPointId = m_solver->m_extractedCells->a[nghbrId].m_pointIds[ltable2[i][1]];
2116 conn[polyId].push_back(gridPointId);
2117 }
2118 if(gridPointId < 0) {
2119 // const MInt counter = m_solver->getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
2120 const MInt counter = m_solver->getAdjacentLeafCells_d2(cellId, 1, nghbrList, layerId);
2121 for(MInt n = 0; n < counter; n++) {
2122 MInt nghbrCell = nghbrList[n];
2123 if(nghbrCell < 0) continue;
2124 if(m_solver->a_isHalo(nghbrCell)) continue;
2125 MInt ncell = cellMap[nghbrCell];
2126 if(ncell < 0) continue;
2127 MInt pc = polyMap[ncell];
2128 // if (m_solver->a_bndryId( nghbrCell )> -1 ) continue;
2129 MInt noPointsNghbr = (pc > -1) ? conn[pc].size() : noNodes;
2130 for(MInt q = 0; q < noPointsNghbr; q++) {
2131 MInt nghbrGridPointId = (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[ncell].m_pointIds[q];
2132 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
2133 MFloat delta = sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0] - ccoords[0])
2134 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1] - ccoords[1])
2135 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2] - ccoords[2]));
2136 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
2137 // cerr << "found another match 1" << endl;
2138 gridPointId = nghbrGridPointId;
2139 conn[polyId].push_back(gridPointId);
2140 break;
2141 }
2142 }
2143 }
2144 }
2145 if(gridPointId < 0) {
2146 gridPointId = m_solver->m_gridPoints->size();
2147 m_solver->m_gridPoints->append();
2148 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
2149 for(MInt j = 0; j < nDim; j++) {
2150 m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] = ccoords[j];
2151 }
2152 conn[polyId].push_back(gridPointId);
2153 }
2154 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes) {
2155 MBool exist = false;
2156 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
2157 ? m_solver->getAssociatedInternalCell(cellId)
2158 : cellId;
2159 for(MInt k = 0; k < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; k++) {
2160 if(m_solver->m_gridPoints->a[gridPointId].m_cellIds[k] == rootId) exist = true;
2161 }
2162 if(!exist) {
2163 m_solver->m_gridPoints->a[gridPointId]
2164 .m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] = rootId;
2165 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
2166 }
2167 }
2168 centerId = gridPointId;
2169
2170 for(MInt k = 0; k < 4; k++) {
2171 MInt cpt = ltable[reverseDir[i]][k];
2172 MInt npt = ltable[reverseDir[i]][dirPlusOne[k]];
2173 gridPointId = -1;
2174 for(MInt j = 0; j < nDim; j++) {
2175 ccoords[j] = m_solver->a_coordinate(cellId, j);
2176 }
2177 ccoords[i / 2] +=
2178 (i % 2 == 0) ? -m_solver->c_cellLengthAtCell(cellId) : m_solver->c_cellLengthAtCell(cellId);
2179 for(MInt j = 0; j < nDim; j++) {
2180 ccoords[j] += signStencil[cpt][j] * F1B4 * m_solver->c_cellLengthAtCell(cellId);
2181 ccoords[j] += signStencil[npt][j] * F1B4 * m_solver->c_cellLengthAtCell(cellId);
2182 }
2183 for(MUint e = 0; e < conn[polyId].size(); e++) {
2184 MFloat delta = F0;
2185 for(MInt j = 0; j < nDim; j++) {
2186 delta += POW2(ccoords[j] - m_solver->m_gridPoints->a[conn[polyId][e]].m_coordinates[j]);
2187 }
2188 delta = sqrt(delta);
2189 if(delta < eps0 * m_solver->c_cellLengthAtCell(cellId)) {
2190 gridPointId = conn[polyId][e];
2191 e = (MInt)conn[polyId].size();
2192 }
2193 }
2194 if(gridPointId < 0 && nghbrId > -1) {
2195 MInt nghbr0 = m_solver->c_childId(m_solver->c_neighborId(cellId, i), ltable[reverseDir[i]][k]);
2196 if(nghbr0 > -1 && m_solver->a_bndryId(nghbr0) > -1) nghbr0 = -1;
2197 MInt nghbrId1 = (nghbr0 < 0) ? -1 : cellMap[nghbr0];
2198 if(nghbr0 > -1 && nghbrId1 > -1) {
2199 gridPointId = m_solver->m_extractedCells->a[nghbrId1].m_pointIds[ltable[reverseDir[i]][dirPlusOne[k]]];
2200 conn[polyId].push_back(gridPointId);
2201 }
2202 }
2203 if(gridPointId < 0) {
2204 // const MInt counter = m_solver->getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
2205 const MInt counter = m_solver->getAdjacentLeafCells_d2(cellId, 1, nghbrList, layerId);
2206 for(MInt n = 0; n < counter; n++) {
2207 MInt nghbrCell = nghbrList[n];
2208 if(nghbrCell < 0) continue;
2209 if(m_solver->a_isHalo(nghbrCell)) continue;
2210 MInt ncell = cellMap[nghbrCell];
2211 if(ncell < 0) continue;
2212 MInt pc = polyMap[ncell];
2213 // if (a_bndryId( nghbrCell )> -1 ) continue;
2214 MInt noPointsNghbr = (pc > -1) ? conn[pc].size() : noNodes;
2215 for(MInt q = 0; q < noPointsNghbr; q++) {
2216 MInt nghbrGridPointId = (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[ncell].m_pointIds[q];
2217 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
2218 MFloat delta =
2219 sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0] - ccoords[0])
2220 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1] - ccoords[1])
2221 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2] - ccoords[2]));
2222 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
2223 // cerr << "found another match 2" << endl;
2224 gridPointId = nghbrGridPointId;
2225 conn[polyId].push_back(gridPointId);
2226 break;
2227 }
2228 }
2229 }
2230 }
2231 if(gridPointId < 0) {
2232 gridPointId = m_solver->m_gridPoints->size();
2233 m_solver->m_gridPoints->append();
2234 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
2235 for(MInt j = 0; j < nDim; j++) {
2236 m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] = ccoords[j];
2237 }
2238 conn[polyId].push_back(gridPointId);
2239 }
2240 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes) {
2241 MBool exist = false;
2242 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
2243 ? m_solver->getAssociatedInternalCell(cellId)
2244 : cellId;
2245 for(MInt q = 0; q < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; q++) {
2246 if(m_solver->m_gridPoints->a[gridPointId].m_cellIds[q] == rootId) exist = true;
2247 }
2248 if(!exist) {
2249 m_solver->m_gridPoints->a[gridPointId]
2250 .m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] = rootId;
2251 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
2252 }
2253 }
2254 edgeIds[k] = gridPointId;
2255 }
2256
2257 for(MInt k = 0; k < 4; k++) {
2258 if((signed)centerId >= m_solver->m_gridPoints->size()
2259 || edgeIds[dirMinusOne[k]] >= m_solver->m_gridPoints->size()
2260 || m_solver->m_extractedCells->a[c].m_pointIds[ltable[i][k]] >= m_solver->m_gridPoints->size()
2261 || edgeIds[k] >= m_solver->m_gridPoints->size()) {
2262 cerr << "node out of range" << endl;
2263 polyFaceSide[i] = false;
2264 }
2265 }
2266 if(!polyFaceSide[i]) continue;
2267 for(MInt k = 0; k < 4; k++) {
2268 facestream[polyId].push_back(0); // number of points for current face
2269 facestream[polyId].push_back(centerId);
2270 facestream[polyId].push_back(edgeIds[dirMinusOne[k]]);
2271 facestream[polyId].push_back(m_solver->m_extractedCells->a[c].m_pointIds[ltable[i][k]]);
2272 facestream[polyId].push_back(edgeIds[k]);
2273 facestream[polyId][pointCntOffset] += 4;
2274
2275 if(m_solver->m_extractedCells->a[c].m_pointIds[ltable[i][k]] < 0)
2276 cerr << "Regular vertex not found" << endl;
2277 if(centerId < 0 || edgeIds[dirMinusOne[k]] < 0
2278 || m_solver->m_extractedCells->a[c].m_pointIds[ltable[i][k]] < 0 || edgeIds[k] < 0)
2279 cerr << "invalid index" << endl;
2280 const MFloat dxb2 = F1B2 * m_solver->c_cellLengthAtCell(cellId);
2281 const MFloat dxeps = 0.000001 * m_solver->c_cellLengthAtCell(cellId);
2282 MFloat pcoord[3];
2283 for(MInt j = 0; j < nDim; j++) {
2284 pcoord[j] = m_solver->a_coordinate(cellId, j);
2285 }
2286 pcoord[i / 2] += (i % 2 == 0) ? -dxb2 : dxb2;
2287
2288 for(MInt e = pointCntOffset + 1; e < (MInt)facestream[polyId].size(); e++) {
2289 gridPointId = facestream[polyId][e];
2290 if(gridPointId < 0) cerr << "Invalid point in connectivity set." << endl;
2291 if(count(conn[polyId].begin(), conn[polyId].end(), gridPointId) == 0)
2292 cerr << "Point not in connectivity set." << endl;
2293 if(fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[i / 2] - pcoord[i / 2]) > dxeps) {
2294 cerr << "vertex out of plane (n) " << cellId << " " << i << " " << e << endl;
2295 }
2296 for(MInt j = 0; j < nDim; j++) {
2297 if(j == i / 2) continue;
2298 if((fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] - pcoord[j]) > dxeps)
2299 && ((fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] - pcoord[j]) - dxb2) > dxeps)) {
2300 cerr << "vertex out of plane (t) " << cellId << " " << i << " " << e << endl;
2301 }
2302 }
2303 }
2304 if(i % 2 == 1)
2305 reverse(facestream[polyId].begin() + pointCntOffset + 1,
2306 facestream[polyId].begin() + facestream[polyId].size());
2307
2308 pointCntOffset = facestream[polyId].size();
2309 facestream[polyId][0]++; // add face
2310 }
2311 }
2312 }
2313 for(MInt i = 0; i < m_noDirs; i++) {
2314 if(!polyFaceSide[i]) {
2315 facestream[polyId].push_back(0); // number of points for current face
2316 for(MInt k = 0; k < 4; k++) {
2317 facestream[polyId].push_back(m_solver->m_extractedCells->a[c].m_pointIds[ltable[i][k]]);
2318 facestream[polyId][pointCntOffset]++;
2319 MInt nb0 = (m_solver->checkNeighborActive(cellId, i) && m_solver->a_hasNeighbor(cellId, i) > 0)
2320 ? m_solver->c_neighborId(cellId, i)
2321 : cellId;
2322 MInt nb1 = (m_solver->checkNeighborActive(cellId, ltable3[i][k])
2323 && m_solver->a_hasNeighbor(cellId, ltable3[i][k]) > 0)
2324 ? m_solver->c_neighborId(cellId, ltable3[i][k])
2325 : cellId;
2326 MInt nb01 =
2327 (m_solver->checkNeighborActive(nb0, ltable3[i][k]) && m_solver->a_hasNeighbor(nb0, ltable3[i][k]) > 0)
2328 ? m_solver->c_neighborId(nb0, ltable3[i][k])
2329 : cellId;
2330 MInt nb10 = (m_solver->checkNeighborActive(nb1, i) && m_solver->a_hasNeighbor(nb1, i) > 0)
2331 ? m_solver->c_neighborId(nb1, i)
2332 : cellId;
2333 MBool isPolyEdge = (m_solver->c_noChildren(nb0) > 0) || (m_solver->c_noChildren(nb1) > 0)
2334 || (m_solver->c_noChildren(nb01) > 0) || (m_solver->c_noChildren(nb10) > 0);
2335 if(polyFaceSide[ltable3[i][k]] || isPolyEdge) {
2336 MInt nghbr00 = m_solver->c_neighborId(cellId, ltable3[i][k]);
2337 MInt nghbr0 = (nghbr00 < 0 || m_solver->c_noChildren(nghbr00))
2338 ? -1
2339 : m_solver->c_childId(nghbr00, ltable[i][dirPlusTwo[k]]);
2340 if(nghbr0 > -1 && m_solver->a_bndryId(nghbr0) > -1) nghbr0 = -1;
2341 MInt nghbrId = (nghbr0 < 0) ? -1 : cellMap[nghbr0];
2342 if(!polyFaceSide[ltable3[i][k]]) nghbrId = -1;
2343 MInt cpt = ltable[i][dirPlusTwo[k]];
2344 MInt npt = ltable[i][dirMinusOne[k]];
2345 MFloat pcoord[3] = {F0, F0, F0};
2346 for(MInt j = 0; j < nDim; j++) {
2347 pcoord[j] = m_solver->a_coordinate(cellId, j);
2348 }
2349 pcoord[ltable3[i][k] / 2] += (ltable3[i][k] % 2 == 0) ? -m_solver->c_cellLengthAtCell(cellId)
2350 : m_solver->c_cellLengthAtCell(cellId);
2351 for(MInt j = 0; j < nDim; j++) {
2352 pcoord[j] += signStencil[cpt][j] * F1B4 * m_solver->c_cellLengthAtCell(cellId);
2353 pcoord[j] += signStencil[npt][j] * F1B4 * m_solver->c_cellLengthAtCell(cellId);
2354 }
2355 MInt gridPointId = -1;
2356 for(MInt e = 0; e < (MInt)conn[polyId].size(); e++) {
2357 MFloat delta = F0;
2358 for(MInt j = 0; j < nDim; j++) {
2359 delta += POW2(pcoord[j] - m_solver->m_gridPoints->a[conn[polyId][e]].m_coordinates[j]);
2360 }
2361 delta = sqrt(delta);
2362 if(delta < eps0 * m_solver->c_cellLengthAtCell(cellId)) {
2363 gridPointId = conn[polyId][e];
2364 break;
2365 }
2366 }
2367 if(gridPointId < 0 && nghbrId > -1) {
2368 gridPointId = m_solver->m_extractedCells->a[nghbrId].m_pointIds[ltable[i][dirMinusOne[k]]];
2369 MFloat delta = F0;
2370 for(MInt j = 0; j < nDim; j++) {
2371 delta += POW2(pcoord[j] - m_solver->m_gridPoints->a[gridPointId].m_coordinates[j]);
2372 }
2373 delta = sqrt(delta);
2374 if(delta > eps0 * m_solver->c_cellLengthAtCell(cellId)) {
2375 cerr << "unexpected coordinate mismatch" << endl;
2376 gridPointId = -1;
2377 } else {
2378 conn[polyId].push_back(gridPointId);
2379 }
2380 }
2381 if(gridPointId < 0) {
2382 // const MInt counter = m_solver->getAdjacentLeafCells<2>( cellId, 1, nghbrList, layerId );
2383 const MInt counter = m_solver->getAdjacentLeafCells_d2(cellId, 1, nghbrList, layerId);
2384 for(MInt n = 0; n < counter; n++) {
2385 MInt nghbrCell = nghbrList[n];
2386 if(nghbrCell < 0) continue;
2387 if(m_solver->a_isHalo(nghbrCell)) continue;
2388 MInt ncell = cellMap[nghbrCell];
2389 if(ncell < 0) continue;
2390 MInt pc = polyMap[ncell];
2391 // if (a_bndryId( nghbrCell )> -1 ) continue;
2392 MInt noPointsNghbr = (pc > -1) ? conn[pc].size() : noNodes;
2393 for(MInt q = 0; q < noPointsNghbr; q++) {
2394 MInt nghbrGridPointId =
2395 (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[ncell].m_pointIds[q];
2396 if(nghbrGridPointId < 0) cerr << "negative grid point id " << endl;
2397 MFloat delta =
2398 sqrt(POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[0] - pcoord[0])
2399 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[1] - pcoord[1])
2400 + POW2(m_solver->m_gridPoints->a[nghbrGridPointId].m_coordinates[2] - pcoord[2]));
2401 if(delta < eps * m_solver->c_cellLengthAtCell(cellId)) {
2402 // cerr << "found another match 3" << endl;
2403 gridPointId = nghbrGridPointId;
2404 conn[polyId].push_back(gridPointId);
2405 break;
2406 }
2407 }
2408 }
2409 }
2410 if(gridPointId < 0) {
2411 gridPointId = m_solver->m_gridPoints->size();
2412 m_solver->m_gridPoints->append();
2413 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 0;
2414 for(MInt j = 0; j < nDim; j++) {
2415 m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] = pcoord[j];
2416 }
2417 conn[polyId].push_back(gridPointId);
2418 }
2419 if(m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells < noNodes) {
2420 MBool exist = false;
2421 MInt rootId = (m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild))
2422 ? m_solver->getAssociatedInternalCell(cellId)
2423 : cellId;
2424 for(MInt q = 0; q < m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells; q++) {
2425 if(m_solver->m_gridPoints->a[gridPointId].m_cellIds[q] == rootId) exist = true;
2426 }
2427 if(!exist) {
2428 m_solver->m_gridPoints->a[gridPointId]
2429 .m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] = rootId;
2430 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
2431 }
2432 }
2433 facestream[polyId].push_back(gridPointId);
2434 facestream[polyId][pointCntOffset]++;
2435 }
2436 }
2437
2438 const MFloat dxb2 = F1B2 * m_solver->c_cellLengthAtCell(cellId);
2439 const MFloat dxeps = 0.000001 * m_solver->c_cellLengthAtCell(cellId);
2440 MFloat pcoord[3];
2441 for(MInt j = 0; j < nDim; j++) {
2442 pcoord[j] = m_solver->a_coordinate(cellId, j);
2443 }
2444 pcoord[i / 2] += (i % 2 == 0) ? -dxb2 : dxb2;
2445 for(MInt e = pointCntOffset + 1; e < (MInt)facestream[polyId].size(); e++) {
2446 MInt gridPointId = facestream[polyId][e];
2447 if(count(conn[polyId].begin(), conn[polyId].end(), gridPointId) == 0)
2448 cerr << "Point not in connectivity set." << endl;
2449 if(fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[i / 2] - pcoord[i / 2]) > dxeps) {
2450 cerr << "vertex out of plane 2 (n) " << cellId << " " << i << " " << e << endl;
2451 }
2452 for(MInt j = 0; j < nDim; j++) {
2453 if(j == i / 2) continue;
2454 if((fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] - pcoord[j]) > dxeps)
2455 && ((fabs(m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] - pcoord[j]) - dxb2) > dxeps)) {
2456 cerr << "vertex out of plane 2 (t) " << cellId << " " << i << " " << e << endl;
2457 }
2458 }
2459 }
2460 if(i % 2 == 1)
2461 reverse(facestream[polyId].begin() + pointCntOffset + 1,
2462 facestream[polyId].begin() + facestream[polyId].size());
2463
2464 pointCntOffset = facestream[polyId].size();
2465 facestream[polyId][0]++; // add face
2466 }
2467 }
2468 sort(conn[polyId].begin(), conn[polyId].end(), less<T>());
2469 conn[polyId].erase(unique(conn[polyId].begin(), conn[polyId].end()), conn[polyId].end());
2470
2471 for(MInt e = 0; e < (MInt)conn[polyId].size(); e++) {
2472 for(MInt f = e + 1; f < (MInt)conn[polyId].size(); f++) {
2473 MInt gridPointId0 = conn[polyId][e];
2474 MInt gridPointId1 = conn[polyId][f];
2475 MFloat delta = F0;
2476 for(MInt j = 0; j < nDim; j++) {
2477 delta += POW2(m_solver->m_gridPoints->a[gridPointId1].m_coordinates[j]
2478 - m_solver->m_gridPoints->a[gridPointId0].m_coordinates[j]);
2479 }
2480 delta = sqrt(delta);
2481 if(delta < 0.0001 * m_solver->c_cellLengthAtCell(cellId)) {
2482 cerr << setprecision(16) << "duplicate point " << e << " " << f << " " << gridPointId0 << " "
2483 << gridPointId1 << " " << delta << endl;
2484 }
2485 }
2486 }
2487
2488 connSize += (uint64_t)conn[polyId].size();
2489 facesSize += (uint64_t)facestream[polyId].size();
2490 }
2491 }
2492
2493
2494 { // check polyhedra
2495 MIntScratchSpace pointRef(50, FUN_, "pointRef");
2496 MInt errorCnt = 0;
2497 const MInt maxNoVertices = 4 * IPOW2(nDim - 1);
2498 // originally: 8
2499 // for multipleLevelSet: 4*IPOW2(nDim-1) +4
2500 for(MInt polyId = 0; polyId < polyCnt; polyId++) {
2501 const MInt cellId = revMap[polyId];
2502 const MInt bndryId = m_solver->a_bndryId(cellId);
2503 MBool errorFlag = false;
2504 if(conn[polyId].empty() || facestream[polyId].empty())
2505 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2506 << "): empty face stream" << endl;
2507 const MFloat cellLength = m_solver->c_cellLengthAtCell(cellId);
2508 const MFloat deltaH = eps0 * cellLength;
2509 // const MFloat deltaH2 = 0.05*cellLength;
2510 const MFloat deltaH3 = eps * cellLength;
2511 const auto noVerts = (MInt)conn[polyId].size();
2512 for(MInt e = 0; e < noVerts; e++) {
2513 if(conn[polyId][e] >= (unsigned)m_solver->m_gridPoints->size()) {
2514 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2515 << "): vertex out of range: " << conn[polyId][e] << " " << conn[polyId][e] << " "
2516 << m_solver->m_gridPoints->size() << endl;
2517 errorFlag = true;
2518 }
2519 for(MInt f = e + 1; f < noVerts; f++) {
2520 MInt gridPointId0 = conn[polyId][e];
2521 MInt gridPointId1 = conn[polyId][f];
2522 MFloat delta = F0;
2523 for(MInt j = 0; j < nDim; j++) {
2524 delta += POW2(m_solver->m_gridPoints->a[gridPointId1].m_coordinates[j]
2525 - m_solver->m_gridPoints->a[gridPointId0].m_coordinates[j]);
2526 }
2527 delta = sqrt(delta);
2528 if((bndryId > -1 && delta < deltaH3) || (bndryId < 0 && delta < deltaH)) {
2529 cerr << setprecision(16) << domainId() << " (" << cellId << "/" << bndryId << "/"
2530 << m_solver->c_globalId(cellId) << "): duplicate point " << e << " " << f << " " << gridPointId0 << " "
2531 << gridPointId1 << " " << delta << " " << delta / cellLength << " /split "
2532 << m_solver->a_hasProperty(cellId, SolverCell::IsSplitCell) << " "
2533 << m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild) << " "
2534 << m_solver->a_hasProperty(cellId, SolverCell::HasSplitFace) << endl;
2535 errorFlag = true;
2536 }
2537 }
2538 }
2539
2540 MInt cnt = 0;
2541 const MInt noFaces = facestream[polyId][cnt++];
2542 if(noFaces < 4 || noFaces > 24) {
2543 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2544 << "): no faces strange: " << noFaces << endl;
2545 errorFlag = true;
2546 }
2547 MInt noEdges = 0;
2548 pointRef.fill(0);
2549 edgeCheck.fill(0);
2550 map<MUint, MInt> pointMap;
2551 MInt tmpCnt = 0;
2552 for(MInt e = 0; e < noVerts; e++) {
2553 if(pointMap.count(conn[polyId][e]) > 0) {
2554 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2555 << "): Duplicate vertex: " << e << " " << pointMap.count(conn[polyId][e]) << endl;
2556 errorFlag = true;
2557 }
2558 pointMap[conn[polyId][e]] = tmpCnt++;
2559 }
2560 if(tmpCnt != noVerts) {
2561 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2562 << "): temp count mismatch: " << tmpCnt << " " << noVerts << endl;
2563 errorFlag = true;
2564 }
2565
2566 for(MInt f = 0; f < noFaces; f++) {
2567 const MInt noPoints = facestream[polyId][cnt++];
2568 if(noPoints < 3 || noPoints > maxNoVertices) {
2569 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2570 << "): no points strange: " << noPoints << endl;
2571 errorFlag = true;
2572 }
2573 for(MInt p = 0; p < noPoints; p++) {
2574 const MUint pointId = facestream[polyId][cnt + p];
2575 if(pointId >= (unsigned)m_solver->m_gridPoints->size()) {
2576 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2577 << "): vertex out of range: " << pointId << " " << m_solver->m_gridPoints->size() << endl;
2578 errorFlag = true;
2579 }
2580 for(MInt q = p + 1; q < noPoints; q++) {
2581 if(facestream[polyId][cnt + q] == pointId) {
2582 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2583 << "): duplicate vertex face: " << p << " " << q << " " << pointId << " "
2584 << facestream[polyId][cnt + q] << endl;
2585 errorFlag = true;
2586 }
2587 }
2588 }
2589 MInt p0 = facestream[polyId][cnt];
2590 MInt p1 = facestream[polyId][cnt + 1];
2591 MInt p2 = facestream[polyId][cnt + 2];
2592 MFloat test = F0;
2593 MInt p2cnt = 3;
2594 MFloat d1 = F0;
2595 MFloat d2 = F0;
2596 for(MInt j = 0; j < nDim; j++) {
2597 test += (m_solver->m_gridPoints->a[p1].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j])
2598 * (m_solver->m_gridPoints->a[p2].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2599 d1 += POW2(m_solver->m_gridPoints->a[p1].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2600 d2 += POW2(m_solver->m_gridPoints->a[p2].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2601 }
2602 d1 = sqrt(d1);
2603 d2 = sqrt(d2);
2604 while(fabs(fabs(test / mMax(1e-12, d1 * d2)) - F1) < 0.1 && p2cnt < noPoints) {
2605 p2 = facestream[polyId][cnt + p2cnt];
2606 d1 = F0;
2607 d2 = F0;
2608 test = F0;
2609 for(MInt j = 0; j < nDim; j++) {
2610 test += (m_solver->m_gridPoints->a[p1].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j])
2611 * (m_solver->m_gridPoints->a[p2].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2612 d1 += POW2(m_solver->m_gridPoints->a[p1].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2613 d2 += POW2(m_solver->m_gridPoints->a[p2].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j]);
2614 }
2615 d1 = sqrt(d1);
2616 d2 = sqrt(d2);
2617 p2cnt++;
2618 }
2619 MFloat a[3]{}, b[3]{}, c[3]{}, d[3]{};
2620 if(noPoints >= 3) {
2621 for(MInt j = 0; j < nDim; j++) {
2622 a[j] = m_solver->m_gridPoints->a[p1].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j];
2623 b[j] = m_solver->m_gridPoints->a[p2].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j];
2624 }
2625 c[0] = a[1] * b[2] - a[2] * b[1];
2626 c[1] = a[2] * b[0] - a[0] * b[2];
2627 c[2] = a[0] * b[1] - a[1] * b[0];
2628 MFloat cabs = F0;
2629 for(MInt j = 0; j < nDim; j++) {
2630 cabs += POW2(c[j]);
2631 }
2632 cabs = sqrt(cabs);
2633 for(MInt j = 0; j < nDim; j++) {
2634 c[j] /= mMax(1e-12, cabs);
2635 }
2636
2637 for(MInt p = 0; p < noPoints; p++) {
2638 const MInt pointId = facestream[polyId][cnt + p];
2639 for(MInt j = 0; j < nDim; j++) {
2640 d[j] =
2641 m_solver->m_gridPoints->a[pointId].m_coordinates[j] - m_solver->m_gridPoints->a[p0].m_coordinates[j];
2642 }
2643 test = F0;
2644 for(MInt j = 0; j < nDim; j++) {
2645 test += c[j] * d[j];
2646 }
2647 // if ( (bndryId > -1 && (fabs(test) > deltaH2) ) || (bndryId < 0 && (fabs(test) > deltaH) ) ) {
2648 if((bndryId < 0 && (fabs(test) > deltaH))) {
2649 cerr << setprecision(16) << domainId() << " (" << cellId << "/" << bndryId << "/"
2650 << m_solver->c_globalId(cellId) << "): vertex out of plane: " << p << " " << test << " "
2651 << test / cellLength << endl;
2652 errorFlag = true;
2653 }
2654 }
2655 }
2656
2657 MInt firstPoint = facestream[polyId][cnt];
2658 for(MInt p = 0; p < noPoints; p++) {
2659 const MInt pointId = facestream[polyId][cnt++];
2660 if(count(conn[polyId].begin(), conn[polyId].end(), pointId) != 1) {
2661 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2662 << "): vertex not found or duplicate: " << p << " " << pointId << " "
2663 << count(conn[polyId].begin(), conn[polyId].end(), pointId) << endl;
2664 errorFlag = true;
2665 }
2666 MInt nextPoint = (p == noPoints - 1) ? firstPoint : facestream[polyId][cnt];
2667 MInt v0 = pointMap[pointId];
2668 MInt v1 = pointMap[nextPoint];
2669 pointRef[v0]++;
2670 if(v0 == v1) {
2671 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2672 << "): edge invalid: " << v0 << " " << v1 << " " << endl;
2673 errorFlag = true;
2674 }
2675 if(v1 < v0) swap(v0, v1);
2676 edgeCheck(v0, v1)++;
2677 }
2678 }
2679
2680 MInt noVerts0 = 0;
2681 for(MInt e = 0; e < noVerts; e++) {
2682 if(pointRef[e] > 0) noVerts0++;
2683 /*
2684 if ( pointRef[e] == 0 ) {
2685 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << c_globalId(cellId)
2686 << "): unused vertex: " << e << " " << pointRef[e] << " " << noVerts
2687 << " /split " << m_solver->a_hasProperty( cellId , SolverCell::IsSplitCell ) << " " <<
2688 m_solver->a_hasProperty( cellId , SolverCell::IsSplitChild )
2689 << " " << m_solver->a_hasProperty( cellId , SolverCell::HasSplitFace ) << endl;
2690 }
2691 else */
2692 if(pointRef[e] != 0 && (pointRef[e] < 2 || pointRef[e] > 4)) {
2693 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2694 << "): cell points not watertight: " << e << " " << pointRef[e] << " " << noVerts << endl;
2695 errorFlag = true;
2696 }
2697 for(MInt f = 0; f < noVerts; f++) {
2698 if(edgeCheck(e, f) != 0 && edgeCheck(e, f) != 2) {
2699 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2700 << "): cell not watertight: " << e << " " << f << " " << edgeCheck(e, f) << endl;
2701 errorFlag = true;
2702 }
2703 if(edgeCheck(e, f) == 2) noEdges++;
2704 }
2705 }
2706
2707 if(noVerts0 - noEdges + noFaces != 2) {
2708 cerr << domainId() << " (" << cellId << "/" << bndryId << "/" << m_solver->c_globalId(cellId)
2709 << "): Euler's formula not fulfilled: " << noVerts << " " << noEdges << " " << noFaces << " /split "
2710 << m_solver->a_hasProperty(cellId, SolverCell::IsSplitCell) << " "
2711 << m_solver->a_hasProperty(cellId, SolverCell::IsSplitChild) << " "
2712 << m_solver->a_hasProperty(cellId, SolverCell::HasSplitFace) << endl;
2713 errorFlag = true;
2714 }
2715
2716 if(bndryId > -1) {
2717 // if ( !m_solver->m_bndryCells->a[bndryId].m_faceVertices.empty() ) errorFlag = true;
2718 }
2719
2720 if(errorFlag && errorCnt < 10) {
2721 ofstream ofl;
2722 string fileName = m_solver->m_solutionOutput + "polyhedron_" + to_string(m_solver->c_globalId(cellId)) + ".vtu";
2723 ofl.open(fileName.c_str(), ios_base::out | ios_base::trunc);
2724 if(ofl.is_open() && ofl.good()) {
2725 ofl << "<?xml version=\"1.0\"?>" << endl;
2726 ofl << R"(<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian">)" << endl;
2727 ofl << "<UnstructuredGrid>" << endl;
2728 ofl << "<Piece NumberOfPoints=\"" << noVerts << "\" NumberOfCells=\"" << 1 << "\">" << endl;
2729 ofl << "<Points>" << endl;
2730 ofl << R"(<DataArray type="Float32" NumberOfComponents="3" format="ascii">)" << endl;
2731 map<MUint, MInt> pointMap2;
2732 tmpCnt = 0;
2733 for(MInt v = 0; v < noVerts; v++) {
2734 MInt gridPointId = conn[polyId][v];
2735 pointMap2[gridPointId] = tmpCnt++;
2736 for(MInt j = 0; j < nDim; j++) {
2737 ofl << m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] << " ";
2738 }
2739 ofl << endl;
2740 }
2741 ofl << "</DataArray>" << endl;
2742 ofl << "</Points>" << endl;
2743 ofl << "<Cells>" << endl;
2744 ofl << R"(<DataArray type="Int32" Name="connectivity" format="ascii">)" << endl;
2745 for(MInt v = 0; v < noVerts; v++) {
2746 ofl << v << " ";
2747 }
2748 ofl << endl;
2749 ofl << "</DataArray>" << endl;
2750 ofl << R"(<DataArray type="Int32" Name="offsets" format="ascii">)" << endl;
2751 ofl << noVerts << endl;
2752 ofl << "</DataArray>" << endl;
2753 ofl << R"(<DataArray type="Int32" Name="types" format="ascii">)" << endl;
2754 ofl << "42" << endl;
2755 ofl << "</DataArray>" << endl;
2756 ofl << R"(<DataArray type="Int32" Name="faces" format="ascii">)" << endl;
2757 MInt fcnt = 0;
2758 MInt noFace = facestream[polyId][fcnt++];
2759 ofl << noFace << " ";
2760 for(MInt f = 0; f < noFace; f++) {
2761 MInt noPoints = facestream[polyId][fcnt++];
2762 ofl << noPoints << " ";
2763 for(MInt p = 0; p < noPoints; p++) {
2764 ofl << pointMap2[facestream[polyId][fcnt++]] << " ";
2765 }
2766 }
2767 ofl << endl;
2768 ofl << "</DataArray>" << endl;
2769 ofl << R"(<DataArray type="Int32" Name="faceoffsets" format="ascii">)" << endl;
2770 ofl << facestream[polyId].size() << endl;
2771 ofl << "</DataArray>" << endl;
2772 ofl << "</Cells>" << endl;
2773 ofl << "</Piece>" << endl;
2774 ofl << "</UnstructuredGrid>" << endl;
2775 ofl << "</VTKFile>" << endl;
2776 ofl.close();
2777 ofl.clear();
2778 }
2779 }
2780 if(errorFlag) errorCnt++;
2781 }
2782 }
2783
2784 /*
2785 const MInt oldNoPoints0 = m_solver->m_gridPoints->size();
2786 MIntScratchSpace pointRefs2(oldNoPoints0, FUN_, "pointRefs2");
2787 pointRefs2.fill(0);
2788 for ( MInt c = 0; c < m_solver->m_extractedCells->size(); c++ ) {
2789 //MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
2790 //MInt bndryId = a_bndryId( cellId );
2791 MInt pc = polyMap[c];
2792 MInt noPoints = ( pc > -1 ) ? conn[pc].size() : noNodes;
2793 for ( MInt q = 0; q < noPoints; q++ ) {
2794 MInt gridPointId = ( pc > -1 ) ? conn[pc][q] : m_solver->m_extractedCells->a[c].m_pointIds[q];
2795 if ( gridPointId < 0 ) continue;
2796 pointRefs2[gridPointId]++;
2797 }
2798 }
2799 MInt inact = 0;
2800 for(MInt gridPointId = 0; gridPointId < oldNoPoints0; gridPointId++ ) {
2801 if ( pointRefs2[gridPointId] == 0 ) {
2802 if ( m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells > 0 ) {
2803 cerr <<"strange " << m_solver->m_gridPoints->a[ gridPointId ].m_noAdjacentCells <<endl;
2804 inact++;
2805 }
2806
2807 }
2808 }
2809 MPI_Allreduce(MPI_IN_PLACE, &inact, 1, MPI_INT, MPI_SUM, mpiComm(), AT_, "MPI_IN_PLACE", "inact" );
2810 cerr << "inact " << inact << endl;
2811 */
2812
2813 return connSize;
2814}
Collector::resetSize
MInt resetSize(MInt inputSize)
Definition: collector.h:45
Collector::append
void append()
Definition: collector.h:153
Collector::size
MInt size()
Definition: collector.h:29
FvCartesianSolverXD::m_splitChilds
std::vector< std::vector< MInt > > m_splitChilds
Definition: fvcartesiansolverxd.h:1458
FvCartesianSolverXD::c_neighborId
MLong c_neighborId(const MInt cellId, const MInt dir, const MBool assertNeighborState=true) const
Returns the grid neighbor id of the grid cell cellId dir.
Definition: fvcartesiansolverxd.h:570
FvCartesianSolverXD::m_bndryCells
Collector< FvBndryCell< nDim, SysEqn > > * m_bndryCells
Definition: fvcartesiansolverxd.h:2314
FvCartesianSolverXD::c_childId
MLong c_childId(const MInt cellId, const MInt pos) const
Returns the grid child id of the grid cell cellId at position pos.
Definition: fvcartesiansolverxd.h:564
FvCartesianSolverXD::c_cellLengthAtCell
MFloat c_cellLengthAtCell(const MInt cellId) const
Returns the length of the cell for level.
Definition: fvcartesiansolverxd.h:577
FvCartesianSolverXD::a_isHalo
MBool a_isHalo(const MInt cellId) const
Returns IsHalo of the cell cellId.
Definition: fvcartesiansolverxd.h:362
FvCartesianSolverXD::getAdjacentLeafCells_d2
virtual MInt getAdjacentLeafCells_d2(const MInt, const MInt, MIntScratchSpace &, MIntScratchSpace &)
Definition: fvcartesiansolverxd.cpp:13947
FvCartesianSolverXD::m_gridPoints
Collector< CartesianGridPoint< nDim > > * m_gridPoints
Definition: fvcartesiansolverxd.h:1395
FvCartesianSolverXD::a_noCells
MInt a_noCells() const
Returns the number of cells.
Definition: fvcartesiansolverxd.h:336
FvCartesianSolverXD::a_hasProperty
MBool a_hasProperty(const MInt cellId, const Cell p) const
Returns grid cell property p of the cell cellId.
Definition: fvcartesiansolverxd.h:588
FvCartesianSolverXD::c_noChildren
MInt c_noChildren(const MInt cellId) const
Returns the number of children of the cell cellId.
Definition: fvcartesiansolverxd.h:540
FvCartesianSolverXD::a_coordinate
MFloat & a_coordinate(const MInt cellId, const MInt dir)
Returns the coordinate of the cell from the fvcellcollector cellId for dimension dir.
Definition: fvcartesiansolverxd.h:425
FvCartesianSolverXD::c_globalId
MLong c_globalId(const MInt cellId) const
Returns the global grid id of the grid cell cellId.
Definition: fvcartesiansolverxd.h:552
FvCartesianSolverXD::m_extractedCells
Collector< PointBasedCell< nDim > > * m_extractedCells
Definition: fvcartesiansolverxd.h:1394
FvCartesianSolverXD::a_hasNeighbor
MInt a_hasNeighbor(const MInt cellId, const MInt dir, const MBool assertNeighborState=true) const
Returns noNeighborIds of the cell CellId for direction dir.
Definition: fvcartesiansolverxd.h:601
FvCartesianSolverXD::gridPointIsInside
virtual MBool gridPointIsInside(MInt, MInt)
Definition: fvcartesiansolverxd.cpp:8159
FvCartesianSolverXD::a_bndryId
MInt & a_bndryId(const MInt cellId)
Returns the bndryId of the cell cellId.
Definition: fvcartesiansolverxd.h:672
FvCartesianSolverXD::m_splitCells
std::vector< MInt > m_splitCells
Definition: fvcartesiansolverxd.h:1457
FvCartesianSolverXD::checkNeighborActive
MBool checkNeighborActive(const MInt cellId, const MInt dir) const
Cecks wether the cell cellId has a valid neighbor in direction dir.
Definition: fvcartesiansolverxd.h:509
FvCartesianSolverXD::getAssociatedInternalCell
const MInt & getAssociatedInternalCell(const MInt &cellId) const
Returns the Id of the split cell, if cellId is a split child.
Definition: fvcartesiansolverxd.h:1081
ScratchSpace::size
size_type size() const
Definition: scratch.h:302
ScratchSpace::fill
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
Solver::m_solutionOutput
MString m_solutionOutput
Definition: solver.h:82
VtkIo::m_noDirs
static constexpr MInt m_noDirs
Definition: iovtk.h:34
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
IPOW2
constexpr MLong IPOW2(MInt x)
Definition: maiaconstants.h:154
MUint
uint32_t MUint
Definition: maiatypes.h:63
MFloat
double MFloat
Definition: maiatypes.h:52
maia::collector_memory::cellId
void const MInt cellId
Definition: collector.h:239
maia::math::distance
MFloat distance(const MFloat *a, const MFloat *b)
Definition: maiamath.h:249
maia::tensor::swap
void swap(Tensor< TT > &a, Tensor< TT > &b)
Non-member swap exchanges the contents of two Tensors.
Definition: tensor.h:752
a
Definition: contexttypes.h:19

◆ writeVtkXmlOutput()

template<MInt nDim, class SysEqn >
void VtkIo< nDim, SysEqn >::writeVtkXmlOutput ( const MChar *  fileName)

Parameters: m_solver->m_vtuWritePointData: specifies whether the flow field is written as point or cell data m_solver->m_vtuCoordinatesThreshold: optionally specifies a bounding box to which the output domain is truncated m_solver->m_vtuLevelThreshold: optionally specifies the maximum cell level to be saved

Opening the .vtu file requires Paraview>=4.0, resp. VTK>=6 If more than 4.2 billion points shall be written the function will rerun with 64 bit integer output For data compression compile with 'WITH_ZLIB' macro

For file-format reference see http://www.vtk.org/Wiki/VTK_XML_Formats

Do not attempt to change to data types used (e.g. float or uint64_t) to XYZ data types since this will break the output!

Author
Lennart Scheiders
Date
01/2018

write flow variables as XML VTK

Author
Lennart Schneiders

Definition at line 180 of file iovtk.cpp.

180 {
181 IF_CONSTEXPR(nDim == 2) {
182 TRACE();
183
184 if(noDomains() > 1) {
185 mTerm(1, AT_, "2D VTK output not yet in massive parallel mode.");
186 }
187
188 MInt noCells = 1;
189 MInt noPoints = 1;
190
191 MBool writePointData = false;
192
193 const MInt baseCellType = 8;
194 const MInt polyCellType = 7;
195
196 MInt ghostCellId, nghbrId, counter;
197
198 //-------------------------------
199 const string dataType64 = "Float64";
200#ifdef DOUBLE_PRECISION_OUTPUT
201 const string dataType = "Float64";
202#else
203 const string dataType = "Float32";
204#endif
205 const string iDataType = "Int32";
206 const string uIDataType = "UInt32";
207 const string uI8DataType = "UInt8";
208 // int inumber;
209 MUint uinumber;
210 // unsigned char ui8number;
211#ifdef DOUBLE_PRECISION_OUTPUT
212 // MFloat fnumber;
213 typedef MFloat flt;
214#else
215 // float fnumber;
216 typedef float flt;
217#endif
218 //-------------------------------
219
220 m_solver->m_firstUseInitializeVtkXmlOutput =
221 initializeVtkXmlOutput(fileName, m_solver->outputDir(), m_solver->m_solutionInterval, m_solver->m_restart,
222 m_solver->m_firstUseInitializeVtkXmlOutput);
223
224 if(domainId() == 0) {
225 cerr << "extracted";
226 cerr << "(" << m_solver->m_extractedCells->size() << ")...";
227 }
228
229 noPoints = m_solver->m_gridPoints->size();
230 MInt noExtractedCells = m_solver->m_extractedCells->size();
231 noCells = noExtractedCells;
232 MInt scratchSize = m_solver->maxRefinementLevel() + 1;
233 MFloatScratchSpace cellLength(scratchSize, AT_, "cellLength");
234
235 for(MInt level = 0; level < m_solver->maxRefinementLevel() + 1; level++) {
236 cellLength.p[level] = m_solver->c_cellLengthAtLevel(level);
237 }
238
239 MBool** pointInside = new MBool*[m_solver->m_bndryCells->size()];
240 for(MInt i = 0; i < m_solver->m_bndryCells->size(); i++) {
241 pointInside[i] = new MBool[IPOW2(nDim)];
242 for(MInt j = 0; j < IPOW2(nDim); j++) {
243 pointInside[i][j] = false;
244 }
245 }
246 const MFloat signStencil[4][2] = {{-F1, -F1}, {F1, -F1}, {-F1, F1}, {F1, F1}};
247 MFloat tmpPoint[nDim];
248
249 for(MInt bndryId = 0; bndryId < m_solver->m_bndryCells->size(); bndryId++) {
250 MInt cellId = m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[bndryId].m_cellId;
251 for(MInt p = 0; p < IPOW2(nDim); p++) {
252 for(MInt i = 0; i < nDim; i++) {
253 tmpPoint[i] =
254 m_solver->a_coordinate(cellId, i) + signStencil[p][i] * F1B2 * cellLength.p[m_solver->a_level(cellId)];
255 }
256 pointInside[bndryId][p] = m_solver->m_geometry->pointIsInside(tmpPoint);
257 }
258 }
259
260 // 4 node ids for each direction
261 const MInt ltable[4][3] = {{2, 3}, {0, 2}, {1, 0}, {3, 1}};
262
263 const MInt pointOrder[4] = {0, 1, 3, 2};
264 const MInt dirOrder[4] = {2, 1, 3, 0};
265 const MInt reverseDir[4] = {1, 0, 3, 2};
266
267 // MFloat edgeIds[ 4 ];
268 MInt centerId, cellId;
269
270 ScratchSpace<MInt> cellTypes(noExtractedCells, AT_, "cellTypes");
271 ScratchSpace<MInt> polyIds(noExtractedCells, AT_, "polyIds");
272 ScratchSpace<MInt> reverseMapping(m_solver->a_noCells(), AT_, "reverseMapping");
273 ScratchSpace<MInt> cutPointIds(m_noDirs * m_solver->m_fvBndryCnd->m_solver->m_bndryCells->size(), AT_,
274 "cutPointIds");
275
276
277 MInt noPolyCells;
278
279
280 for(MInt i = 0; i < m_noDirs * m_solver->m_fvBndryCnd->m_solver->m_bndryCells->size(); i++) {
281 cutPointIds.p[i] = -1;
282 }
283
284 if(domainId() == 0) {
285 cerr << "prepare...";
286 }
287 for(MInt c = 0; c < m_solver->a_noCells(); c++) {
288 reverseMapping.p[c] = -1;
289 }
290 for(MInt c = 0; c < noExtractedCells; c++) {
291 reverseMapping.p[m_solver->m_extractedCells->a[c].m_cellId] = c;
292 }
293
294
295 // delete inactive grid points
296 for(MInt p = 0; p < noPoints; p++) {
297 MBool pointIsActive = true;
298 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
299 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
300 if((m_solver->a_bndryId(cellId) > -1) && (reverseMapping.p[cellId] > -1)) {
301 for(MInt k = 0; k < IPOW2(nDim); k++) {
302 if(m_solver->m_extractedCells->a[reverseMapping.p[cellId]].m_pointIds[k] == p) {
303 if(pointInside[m_solver->a_bndryId(cellId)][k]) {
304 pointIsActive = false;
305 break;
306 }
307 }
308 }
309 }
310 }
311 if(!pointIsActive) {
312 noPoints--;
313
314 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
315 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
316 if(reverseMapping.p[cellId] > -1) {
317 for(MInt k = 0; k < IPOW2(nDim); k++) {
318 if(m_solver->m_extractedCells->a[reverseMapping.p[cellId]].m_pointIds[k] == p) {
319 m_solver->m_extractedCells->a[reverseMapping.p[cellId]].m_pointIds[k] = -1;
320 }
321 }
322 }
323 }
324
325 for(MInt n = 0; n < m_solver->m_gridPoints->a[noPoints].m_noAdjacentCells; n++) {
326 cellId = m_solver->m_gridPoints->a[noPoints].m_cellIds[n];
327 if(reverseMapping.p[cellId] > -1) {
328 for(MInt k = 0; k < IPOW2(nDim); k++) {
329 if(m_solver->m_extractedCells->a[reverseMapping.p[cellId]].m_pointIds[k] == noPoints) {
330 m_solver->m_extractedCells->a[reverseMapping.p[cellId]].m_pointIds[k] = p;
331 }
332 }
333 }
334 m_solver->m_gridPoints->a[p].m_cellIds[n] = m_solver->m_gridPoints->a[noPoints].m_cellIds[n];
335 m_solver->m_gridPoints->a[p].m_noAdjacentCells = m_solver->m_gridPoints->a[noPoints].m_noAdjacentCells;
336 }
337 for(MInt i = 0; i < nDim; i++) {
338 m_solver->m_gridPoints->a[p].m_coordinates[i] = m_solver->m_gridPoints->a[noPoints].m_coordinates[i];
339 }
340 m_solver->m_gridPoints->resetSize(noPoints);
341 p--;
342 }
343 }
344
345 //---
346
347 MBool isPolyCell;
348 noPolyCells = 0;
349
350 for(MInt c = 0; c < noExtractedCells; c++) {
351 cellId = m_solver->m_extractedCells->a[c].m_cellId;
352 cellTypes.p[c] = baseCellType;
353 polyIds.p[c] = -1;
354 isPolyCell = false;
355 if(m_solver->a_bndryId(cellId) > -1) {
356 isPolyCell = true;
357 } else {
358 for(MInt i = 0; i < m_noDirs; i++) {
359 if(m_solver->a_hasNeighbor(cellId, i) > 0) {
360 if(m_solver->c_noChildren(m_solver->c_neighborId(cellId, i)) > 0) {
361 isPolyCell = true;
362 }
363 }
364 }
365 }
366 if(isPolyCell) {
367 cellTypes.p[c] = polyCellType;
368 polyIds.p[c] = noPolyCells;
369 noPolyCells++;
370 }
371 }
372
373 MBool isPolyFace;
374 // MBool pointExists;
375 MInt eCell;
376 const MInt maxNoExtraPoints = 8;
377 const MInt maxNoPolyCells = noPolyCells;
378 ScratchSpace<MInt> polyCellLink(maxNoPolyCells, AT_, "polyCellLink");
379 ScratchSpace<MInt> extraPoints(maxNoPolyCells * maxNoExtraPoints, AT_, "extraPoints");
380 ScratchSpace<MInt> noExtraPoints(maxNoPolyCells, AT_, "noExtraPoints");
381 ScratchSpace<MInt> noInternalPoints(maxNoPolyCells, AT_, "noInternalPoints");
382 MInt offset, pointCount; // faceCount;
383
384 for(MInt c = 0; c < noExtractedCells; c++) {
385 if(cellTypes.p[c] == polyCellType) {
386 polyCellLink.p[polyIds.p[c]] = c;
387 }
388 }
389
390
391 const MInt noInternalGridPoints = m_solver->m_gridPoints->size();
392
393 for(MInt pc = 0; pc < noPolyCells; pc++) {
394 eCell = polyCellLink.p[pc];
395 cellId = m_solver->m_extractedCells->a[eCell].m_cellId;
396 noInternalPoints.p[pc] = IPOW2(nDim);
397 noExtraPoints.p[pc] = 0;
398
399 // a) cut cells at the boundary
400 if(m_solver->a_bndryId(cellId) > -1) {
401 MInt bndryId = m_solver->a_bndryId(cellId);
402
403 noInternalPoints.p[pc] = 0;
404
405 for(MInt i = 0; i < m_noDirs; i++) {
406 MInt p = pointOrder[i];
407 MInt dir = dirOrder[i];
408
409 // 0. determine internal vertices
410 if(!pointInside[bndryId][p]) {
411 extraPoints.p[pc * maxNoExtraPoints + noExtraPoints.p[pc]] =
412 m_solver->m_extractedCells->a[eCell].m_pointIds[p];
413 noExtraPoints.p[pc]++;
414 }
415
416 // 1. add cut points
417 for(MInt cp = 0; cp < m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[bndryId].m_srfcs[0]->m_noCutPoints;
418 cp++) {
419 if(m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[bndryId].m_srfcs[0]->m_cutEdge[cp] == dir) {
420 MInt gridPointId = -1;
421 if(m_solver->a_hasNeighbor(cellId, dir) > 0) {
422 if(m_solver->a_bndryId(m_solver->c_neighborId(cellId, dir)) > -1) {
423 gridPointId =
424 cutPointIds
425 .p[m_noDirs * m_solver->a_bndryId(m_solver->c_neighborId(cellId, dir)) + reverseDir[dir]];
426 if(gridPointId > -1) {
427 m_solver->m_gridPoints->a[gridPointId]
428 .m_cellIds[m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells] = cellId;
429 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells++;
430 }
431 }
432 }
433 if(gridPointId >= noPoints) {
434 mTerm(1, AT_, "Error 0 in FvMbSolver3D::writeVtkXmlOutput(..). Quit.");
435 }
436 if(gridPointId < 0) {
437 gridPointId = m_solver->m_gridPoints->size();
438 m_solver->m_gridPoints->append();
439 noPoints++;
440 m_solver->m_gridPoints->a[gridPointId].m_noAdjacentCells = 1;
441 m_solver->m_gridPoints->a[gridPointId].m_cellIds[0] = cellId;
442 for(MInt j = 0; j < nDim; j++) {
443 m_solver->m_gridPoints->a[gridPointId].m_coordinates[j] =
444 m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
445 .m_srfcs[0]
446 ->m_cutCoordinates[cp][j];
447 }
448 }
449 cutPointIds.p[m_noDirs * bndryId + dir] = gridPointId;
450
451 extraPoints.p[pc * maxNoExtraPoints + noExtraPoints.p[pc]] = gridPointId;
452 noExtraPoints.p[pc]++;
453 }
454 }
455 }
456 }
457 // b) cells at fine/coarse mesh interfaces
458 else {
459 noInternalPoints.p[pc] = 0;
460
461 for(MInt i = 0; i < m_noDirs; i++) {
462 MInt p = pointOrder[i];
463 MInt dir = dirOrder[i];
464
465 extraPoints.p[pc * maxNoExtraPoints + noExtraPoints.p[pc]] =
466 m_solver->m_extractedCells->a[eCell].m_pointIds[p];
467 noExtraPoints.p[pc]++;
468
469 isPolyFace = false;
470 if(m_solver->a_hasNeighbor(cellId, dir) > 0) {
471 if(m_solver->c_noChildren(m_solver->c_neighborId(cellId, dir)) > 0) {
472 isPolyFace = true;
473 }
474 }
475
476 if(isPolyFace) {
477 // store four faces (polygon numbering!)
478 nghbrId = reverseMapping.p[m_solver->c_childId(m_solver->c_neighborId(cellId, dir), ltable[i][0])];
479 if(nghbrId < 0) {
480 cerr << endl
481 << cellId << " " << m_solver->c_neighborId(cellId, dir) << " "
482 << m_solver->c_childId(m_solver->c_neighborId(cellId, dir), ltable[i][0]) << " / "
483 << m_solver->a_level(cellId) << " " << m_solver->a_level(m_solver->c_neighborId(cellId, dir)) << " "
484 << m_solver->a_level(m_solver->c_childId(m_solver->c_neighborId(cellId, dir), ltable[i][0])) << " / "
485 << m_solver->c_noChildren(cellId) << " "
486 << m_solver->c_noChildren(m_solver->c_neighborId(cellId, dir)) << " "
487 << m_solver->c_noChildren(m_solver->c_childId(m_solver->c_neighborId(cellId, dir), ltable[i][0]))
488 << endl;
489 cerr << m_solver->c_childId(cellId, 0) << " " << m_solver->c_childId(cellId, 1) << " "
490 << m_solver->c_childId(cellId, 2) << " " << m_solver->c_childId(cellId, 3) << endl;
491 mTerm(1, AT_, "Error 1 in FvMbSolver2D::writeVtkXmlOutput(..). Quit.");
492 }
493 centerId = m_solver->m_extractedCells->a[nghbrId].m_pointIds[ltable[i][1]];
494 extraPoints.p[pc * maxNoExtraPoints + noExtraPoints.p[pc]] = centerId;
495 noExtraPoints.p[pc]++;
496
497 if(m_solver->m_gridPoints->a[centerId].m_noAdjacentCells < IPOW2(nDim)) {
498 m_solver->m_gridPoints->a[centerId].m_cellIds[m_solver->m_gridPoints->a[centerId].m_noAdjacentCells] =
499 cellId;
500 m_solver->m_gridPoints->a[centerId].m_noAdjacentCells++;
501 } else {
502 mTerm(1, AT_, "Error 2 in FvMbSolver2D::writeVtkXmlOutput(..). Quit.");
503 }
504 }
505 }
506 }
507 }
508
509 pointCount = 0;
510 for(MInt c = 0; c < noExtractedCells; c++) {
511 MInt pc = polyIds.p[c];
512 if(pc > -1) {
513 pointCount += noInternalPoints.p[pc];
514 pointCount += noExtraPoints.p[pc];
515 } else {
516 pointCount += IPOW2(nDim);
517 }
518 }
519
520
521 MFloat gradU[2][2];
522 for(MInt c = 0; c < m_solver->m_noActiveCells; c++) {
523 cellId = m_solver->m_activeCellIds[c];
524 for(MInt i = 0; i < nDim; i++) {
525 for(MInt j = 0; j < nDim; j++) {
526 gradU[i][j] =
527 m_solver->a_slope(cellId, sysEqn().PV->VV[i], j) * m_solver->m_referenceLength / m_solver->m_UInfinity;
528 }
529 }
530 MFloat omega = F1B2 * (gradU[1][0] - gradU[0][1]);
531 MFloat S = 0;
532 for(MInt i = 0; i < nDim; i++) {
533 for(MInt j = 0; j < nDim; j++) {
534 S += POW2(F1B2 * (gradU[i][j] + gradU[j][i]));
535 }
536 }
537 MFloat Q = POW2(omega) - F1B2 * POW2(S);
538 m_solver->a_slope(cellId, 2, 0) = Q;
539 m_solver->a_slope(cellId, 0, 0) = omega;
540 }
541
542
543 //=============================================
544 //================ GATHER =====================
545 //=============================================
546
547 if(domainId() == 0) {
548 cerr << "gather...";
549 }
550
551 //-----------
552 // points
553 ScratchSpace<flt> points(3 * noPoints, AT_, "points");
554 for(MInt p = 0; p < noPoints; p++) {
555 for(MInt i = 0; i < nDim; i++) {
556 points.p[3 * p + i] = (flt)m_solver->m_gridPoints->a[p].m_coordinates[i];
557 }
558 points.p[3 * p + 2] = (flt)F0;
559 }
560
561 //-----------
562 // connectivity
563 ScratchSpace<MUint> connectivity(pointCount, AT_, "connectivity");
564 counter = 0;
565 for(MInt c = 0; c < noExtractedCells; c++) {
566 MInt pc = polyIds.p[c];
567 if(pc > -1) {
568 for(MInt p = 0; p < noExtraPoints.p[pc]; p++) {
569 connectivity.p[counter] = (MUint)extraPoints.p[pc * maxNoExtraPoints + p];
570 counter++;
571 }
572 } else {
573 for(MInt p = 0; p < IPOW2(nDim); p++) {
574 connectivity.p[counter] = (MUint)m_solver->m_extractedCells->a[c].m_pointIds[p];
575 counter++;
576 }
577 }
578 }
579 if(counter != pointCount) {
580 mTerm(1, AT_, "E1");
581 }
582
583 //-----------
584 // offsets
585 ScratchSpace<MUint> offsets(noCells, AT_, "offsets");
586 counter = 0;
587 for(MInt c = 0; c < noExtractedCells; c++) {
588 MInt pc = polyIds.p[c];
589 if(pc > -1) {
590 counter += noInternalPoints.p[pc] + noExtraPoints.p[pc];
591 } else {
592 counter += IPOW2(nDim);
593 }
594
595 offsets.p[c] = (MUint)counter;
596 }
597
598
599 //-----------
600 // types
601 ScratchSpace<unsigned char> types(noCells, AT_, "types");
602 for(MInt c = 0; c < noExtractedCells; c++) {
603 types.p[c] = (unsigned char)cellTypes.p[c];
604 }
605
606 //-----------
607 // cellIds
608 ScratchSpace<MInt> cellIds(noCells, AT_, "cellIds");
609 for(MInt c = 0; c < noExtractedCells; c++) {
610 cellIds.p[c] = m_solver->m_extractedCells->a[c].m_cellId;
611 }
612
613 //-----------
614 // bndryIds
615 ScratchSpace<MInt> bndryIds(noCells, AT_, "bndryIds");
616 for(MInt c = 0; c < noExtractedCells; c++) {
617 bndryIds.p[c] = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
618 }
619 /*
620 //-----------
621 // drho_dt
622 ScratchSpace<MFloat> drho_dt(noCells, AT_, "drho_dt" );
623 for( MInt c = 0; c < noExtractedCells; c++ ) {
624 drho_dt.p[ c ] = m_rightHandSide[ m_solver->m_extractedCells->a[ c ].m_cellId ][ CV->RHO ];
625 }
626 */
627 //-----------
628 // levelSetFunction
629 ScratchSpace<flt> levelSetFunction((m_solver->m_levelSet) ? noCells : 0, AT_, "levelSetFunction");
630 if(m_solver->m_levelSet) {
631 for(MInt c = 0; c < noExtractedCells; c++) {
632 if(m_solver->m_combustion) {
633 levelSetFunction.p[c] = m_solver->a_levelSetFunction(m_solver->m_extractedCells->a[c].m_cellId, 0);
634 } else if(!m_solver->m_levelSetMb) {
635 levelSetFunction.p[c] = m_solver->a_levelSetFunction(m_solver->m_extractedCells->a[c].m_cellId, 0);
636 }
637 }
638 }
639
640 //-----------
641 // additional cellData/pointData
642 MInt asize = noCells;
643 if(writePointData) {
644 asize = noPoints;
645 }
646 ScratchSpace<flt> pressure(asize, AT_, "pressure");
647 ScratchSpace<flt> velocity(3 * asize, AT_, "velocity");
648 ScratchSpace<flt> density(asize, AT_, "density");
649 ScratchSpace<flt> vorticity(asize, AT_, "vorticity");
650 ScratchSpace<flt> progressVariable((m_solver->m_combustion) ? asize : 0, AT_, "progressVariable");
651 // ScratchSpace<flt> Qcriterion(asize, AT_, "Qcriterion" );
652
653
654 MFloat tmp;
655 MFloat weights[/*IPOW2(nDim)*/ 2 * 2];
656 if(writePointData) {
657 for(MInt p = 0; p < noInternalGridPoints; p++) {
658 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
659 weights[n] = F0;
660 }
661 tmp = F0;
662 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
663 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
664 if(m_solver->a_bndryId(cellId) > -1) {
665 weights[n] = F1
666 / (sqrt(POW2(m_solver->a_coordinate(cellId, 0)
667 + m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
668 .m_coordinates[0]
669 - m_solver->m_gridPoints->a[p].m_coordinates[0])
670 + POW2(m_solver->a_coordinate(cellId, 1)
671 + m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
672 .m_coordinates[1]
673 - m_solver->m_gridPoints->a[p].m_coordinates[1])));
674 tmp += weights[n];
675 } else {
676 weights[n] =
677 F1
678 / (sqrt(POW2(m_solver->a_coordinate(cellId, 0) - m_solver->m_gridPoints->a[p].m_coordinates[0])
679 + POW2(m_solver->a_coordinate(cellId, 1) - m_solver->m_gridPoints->a[p].m_coordinates[1])));
680 tmp += weights[n];
681 }
682 }
683 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
684 weights[n] /= tmp;
685 }
686 //---
687 pressure.p[p] = (flt)F0;
688 density.p[p] = (flt)F0;
689 vorticity.p[p] = (flt)F0;
690 // Qcriterion.p[ p ] = (flt)F0 ;
691 if(m_solver->m_combustion) {
692 progressVariable.p[p] = (flt)F0;
693 }
694 for(MInt i = 0; i < nDim; i++) {
695 velocity.p[3 * p + i] = (flt)F0;
696 }
697 velocity.p[3 * p + 2] = (flt)F0;
698
699 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
700 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
701 pressure.p[p] += (flt)weights[n] * m_solver->a_pvariable(cellId, sysEqn().PV->P);
702 density.p[p] += (flt)weights[n] * m_solver->a_pvariable(cellId, sysEqn().PV->RHO);
703 vorticity.p[p] += (flt)weights[n] * m_solver->a_slope(cellId, 0, 0);
704 IF_CONSTEXPR(hasPV_C<SysEqn>::value) {
705 if(m_solver->m_combustion) {
706 progressVariable.p[p] += (flt)weights[n] * m_solver->a_pvariable(cellId, sysEqn().PV->C);
707 }
708 }
709 for(MInt i = 0; i < nDim; i++) {
710 velocity.p[3 * p + i] += (flt)weights[n] * m_solver->a_pvariable(cellId, sysEqn().PV->VV[i]);
711 }
712 }
713 }
714 for(MInt p = noInternalGridPoints; p < noPoints; p++) {
715 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
716 weights[n] = F0;
717 }
718 tmp = F0;
719 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
720 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
721 if(m_solver->a_bndryId(cellId) < 0) {
722 continue;
723 }
724 weights[n] = F1
725 / (sqrt(POW2(m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
726 .m_srfcs[0]
727 ->m_coordinates[0]
728 - m_solver->m_gridPoints->a[p].m_coordinates[0])
729 + POW2(m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
730 .m_srfcs[0]
731 ->m_coordinates[1]
732 - m_solver->m_gridPoints->a[p].m_coordinates[1])));
733 tmp += weights[n];
734 }
735 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
736 weights[n] /= tmp;
737 }
738 pressure.p[p] = (flt)F0;
739 density.p[p] = (flt)F0;
740 vorticity.p[p] = (flt)F0;
741 // Qcriterion.p[ p ] = (flt)F0 ;
742 if(m_solver->m_combustion) {
743 progressVariable.p[p] = (flt)F0;
744 }
745 for(MInt i = 0; i < nDim; i++) {
746 velocity.p[3 * p + i] = (flt)F0;
747 }
748 velocity.p[3 * p + 2] = F0;
749 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
750 cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
751 if(m_solver->a_bndryId(cellId) < 0) {
752 continue;
753 }
754 ghostCellId = m_solver->m_fvBndryCnd->m_solver->m_bndryCells->a[m_solver->a_bndryId(cellId)]
755 .m_srfcVariables[0]
756 ->m_ghostCellId;
757 pressure.p[p] +=
758 (flt)weights[n] * F1B2
759 * (m_solver->a_pvariable(cellId, sysEqn().PV->P) + m_solver->a_pvariable(ghostCellId, sysEqn().PV->P));
760 density.p[p] += (flt)weights[n] * F1B2
761 * (m_solver->a_pvariable(cellId, sysEqn().PV->RHO)
762 + m_solver->a_pvariable(ghostCellId, sysEqn().PV->RHO));
763 vorticity.p[p] += (flt)weights[n] * m_solver->a_slope(cellId, 0, 0);
764 IF_CONSTEXPR(hasPV_C<SysEqn>::value) {
765 if(m_solver->m_combustion) {
766 progressVariable.p[p] += (flt)weights[n] * F1B2
767 * (m_solver->a_pvariable(cellId, sysEqn().PV->C)
768 + m_solver->a_pvariable(ghostCellId, sysEqn().PV->C));
769 }
770 }
771 for(MInt i = 0; i < nDim; i++) {
772 velocity.p[3 * p + i] += (flt)weights[n] * F1B2
773 * (m_solver->a_pvariable(cellId, sysEqn().PV->VV[i])
774 + m_solver->a_pvariable(ghostCellId, sysEqn().PV->VV[i]));
775 }
776 }
777 }
778 } else {
779 for(MInt c = 0; c < noExtractedCells; c++) {
780 cellId = m_solver->m_extractedCells->a[c].m_cellId;
781 pressure.p[c] = (flt)m_solver->a_pvariable(cellId, sysEqn().PV->P);
782 density.p[c] = (flt)m_solver->a_pvariable(cellId, sysEqn().PV->RHO);
783 vorticity.p[c] = (flt)m_solver->a_slope(cellId, 0, 0);
784 IF_CONSTEXPR(hasPV_C<SysEqn>::value) {
785 if(m_solver->m_combustion) {
786 progressVariable.p[c] = (flt)m_solver->a_pvariable(cellId, sysEqn().PV->C);
787 }
788 }
789 for(MInt i = 0; i < nDim; i++) {
790 velocity.p[3 * c + i] = (flt)m_solver->a_pvariable(cellId, sysEqn().PV->VV[i]);
791 }
792 velocity.p[3 * c + 2] = (flt)F0;
793 }
794 }
795
796
797 if(domainId() == 0) {
798 cerr << "write...";
799 }
800
801 ofstream ofl;
802 ofl.open(fileName);
803
804 if(ofl) {
805 offset = 0;
806
807 //================== VTKFile =================
808 ofl << "<?xml version=\"1.0\"?>" << endl;
809 ofl << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
810 ofl << "<UnstructuredGrid>" << endl;
811
812 //================== FieldData =================
813 ofl << "<DataArray type=\"" << dataType
814 << "\" Name=\"timeStep\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << m_solver->timeStep()
815 << "\" RangeMax=\"" << m_solver->timeStep() << "\">" << endl;
816 ofl << m_solver->timeStep() << endl;
817 ofl << "</DataArray>" << endl;
818 ofl << "<FieldData>" << endl;
819 ofl << "<DataArray type=\"" << dataType << "\" Name=\"refTime\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
820 << m_solver->m_timeRef << "\" RangeMax=\"" << m_solver->m_timeRef << "\">" << endl;
821 ofl << m_solver->m_timeRef << endl;
822 ofl << "</DataArray>" << endl;
823 if(m_solver->m_outputPhysicalTime) {
824 ofl << "<DataArray type=\"" << dataType << "\" Name=\"time\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
825 << m_solver->m_physicalTime << "\" RangeMax=\"" << m_solver->m_physicalTime << "\">" << endl;
826 ofl << m_solver->m_physicalTime << endl;
827 ofl << "</DataArray>" << endl;
828 ofl << "<DataArray type=\"" << dataType
829 << "\" Name=\"internalTime\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << m_solver->m_time
830 << "\" RangeMax=\"" << m_solver->m_time << "\">" << endl;
831 ofl << m_solver->m_time << endl;
832 } else {
833 ofl << "<DataArray type=\"" << dataType << "\" Name=\"time\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
834 << m_solver->m_time << "\" RangeMax=\"" << m_solver->m_time << "\">" << endl;
835 ofl << m_solver->m_time << endl;
836 ofl << "</DataArray>" << endl;
837 ofl << "<DataArray type=\"" << dataType
838 << "\" Name=\"physicalTime\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << m_solver->m_physicalTime
839 << "\" RangeMax=\"" << m_solver->m_physicalTime << "\">" << endl;
840 ofl << m_solver->m_physicalTime << endl;
841 }
842 ofl << "</DataArray>" << endl;
843 ofl << "<DataArray type=\"" << dataType << "\" Name=\"Ma\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
844 << m_solver->m_Ma << "\" RangeMax=\"" << m_solver->m_Ma << "\">" << endl;
845 ofl << m_solver->m_Ma << endl;
846 ofl << "</DataArray>" << endl;
847 ofl << "<DataArray type=\"" << dataType << "\" Name=\"Re\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
848 << m_solver->m_Re << "\" RangeMax=\"" << m_solver->m_Re << "\">" << endl;
849 ofl << m_solver->m_Re << endl;
850 ofl << "</DataArray>" << endl;
851 ofl << "<DataArray type=\"" << dataType << "\" Name=\"Pr\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
852 << m_solver->m_Pr << "\" RangeMax=\"" << m_solver->m_Pr << "\">" << endl;
853 ofl << m_solver->m_Pr << endl;
854 ofl << "</DataArray>" << endl;
855 ofl << "<DataArray type=\"" << dataType << "\" Name=\"gamma\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
856 << m_solver->sysEqn().gamma_Ref() << "\" RangeMax=\"" << m_solver->sysEqn().gamma_Ref() << "\">" << endl;
857 ofl << m_solver->sysEqn().gamma_Ref() << endl;
858 ofl << "</DataArray>" << endl;
859 ofl << "<DataArray type=\"" << dataType << "\" Name=\"CFL\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
860 << m_solver->m_cfl << "\" RangeMax=\"" << m_solver->m_cfl << "\">" << endl;
861 ofl << m_solver->m_cfl << endl;
862 ofl << "</DataArray>" << endl;
863 ofl << "<DataArray type=\"" << dataType << "\" Name=\"uInf\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
864 << m_solver->m_UInfinity << "\" RangeMax=\"" << m_solver->m_UInfinity << "\">" << endl;
865 ofl << m_solver->m_UInfinity << endl;
866 ofl << "</DataArray>" << endl;
867 ofl << "<DataArray type=\"" << dataType << "\" Name=\"vInf\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
868 << m_solver->m_VInfinity << "\" RangeMax=\"" << m_solver->m_VInfinity << "\">" << endl;
869 ofl << m_solver->m_VInfinity << endl;
870 ofl << "</DataArray>" << endl;
871 /* ofl << "<DataArray type=\"" << dataType << "\" Name=\"wInf\" format=\"ascii\"
872 NumberOfTuples=\"1\" RangeMin=\"" << m_solver->m_WInfinity << "\" RangeMax=\"" << m_solver->m_WInfinity << "\">"
873 << endl; ofl << m_solver->m_WInfinity << endl; ofl << "</DataArray>" << endl;*/
874 ofl << "<DataArray type=\"" << dataType << "\" Name=\"pInf\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
875 << m_solver->m_PInfinity << "\" RangeMax=\"" << m_solver->m_PInfinity << "\">" << endl;
876 ofl << m_solver->m_PInfinity << endl;
877 ofl << "</DataArray>" << endl;
878 ofl << "<DataArray type=\"" << dataType << "\" Name=\"rhoInf\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
879 << m_solver->m_rhoInfinity << "\" RangeMax=\"" << m_solver->m_rhoInfinity << "\">" << endl;
880 ofl << m_solver->m_rhoInfinity << endl;
881 ofl << "</DataArray>" << endl;
882 ofl << "</FieldData>" << endl;
883 if(m_solver->m_combustion) {
884 ofl << "<DataArray type=\"" << dataType
885 << "\" Name=\"flSpeed\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << m_solver->m_flameSpeed
886 << "\" RangeMax=\"" << m_solver->m_flameSpeed << "\">" << endl;
887 ofl << m_solver->m_flameSpeed << endl;
888 ofl << "</DataArray>" << endl;
889 ofl << "<DataArray type=\"" << dataType
890 << "\" Name=\"lamflTh\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
891 << m_solver->m_laminarFlameThickness << "\" RangeMax=\"" << m_solver->m_laminarFlameThickness << "\">"
892 << endl;
893 ofl << m_solver->m_laminarFlameThickness << endl;
894 ofl << "</DataArray>" << endl;
895 if(m_solver->m_forcing) {
896 // flame strouhal number = 2*pi*f/Lf
897 MInt flameStr = m_solver->m_flameStrouhal / (F2 * PI);
898 ofl << "<DataArray type=\"" << dataType
899 << "\" Name=\"flStr\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << flameStr << "\" RangeMax=\""
900 << flameStr << "\">" << endl;
901 ofl << flameStr << endl;
902 ofl << "</DataArray>" << endl;
903 // wave number based on Lf or 1 omega = 2*pi/Lf (= m_flameStrouhal)
904 ofl << "<DataArray type=\"" << dataType
905 << "\" Name=\"waveN\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\"" << m_solver->m_flameStrouhal
906 << "\" RangeMax=\"" << m_solver->m_flameStrouhal << "\">" << endl;
907 ofl << m_solver->m_flameStrouhal << endl;
908 ofl << "</DataArray>" << endl;
909 ofl << "<DataArray type=\"" << dataType
910 << "\" Name=\"forcAmpl\" format=\"ascii\" NumberOfTuples=\"1\" RangeMin=\""
911 << m_solver->m_forcingAmplitude << "\" RangeMax=\"" << m_solver->m_forcingAmplitude << "\">" << endl;
912 ofl << m_solver->m_forcingAmplitude << endl;
913 ofl << "</DataArray>" << endl;
914 }
915 }
916
917 //================== /FieldData =================
918
919 ofl << "<Piece NumberOfPoints=\"" << noPoints << "\" NumberOfCells=\"" << noCells << "\">" << endl;
920
921
922 //================== Points =================
923 ofl << "<Points>" << endl;
924 ofl << "<DataArray type=\"" << dataType << "\" NumberOfComponents=\"3\" format=\"appended\" offset=\"" << offset
925 << "\"/>" << endl;
926 offset += points.m_memsize + sizeof(MUint);
927 ofl << "</Points>" << endl;
928 //================== /Points =================
929
930
931 //================== Cells =================
932 ofl << "<Cells>" << endl;
933
934 ofl << "<DataArray type=\"" << uIDataType << "\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
935 << "\"/>" << endl;
936 offset += connectivity.m_memsize + sizeof(MUint);
937
938 ofl << "<DataArray type=\"" << uIDataType << "\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
939 << "\"/>" << endl;
940 offset += offsets.m_memsize + sizeof(MUint);
941
942 ofl << "<DataArray type=\"" << uI8DataType << "\" Name=\"types\" format=\"appended\" offset=\"" << offset
943 << "\"/>" << endl;
944 offset += types.m_memsize + sizeof(MUint);
945
946 ofl << "</Cells>" << endl;
947 //================== /Cells =================
948
949
950 //================== CellData/PointData =================
951 ofl << "<CellData Scalars=\"scalars\">" << endl;
952
953 ofl << "<DataArray type=\"" << iDataType << "\" Name=\"cellIds\" format=\"appended\" offset=\"" << offset
954 << "\"/>" << endl;
955 offset += cellIds.m_memsize + sizeof(MUint);
956
957 ofl << "<DataArray type=\"" << iDataType << "\" Name=\"bndryIds\" format=\"appended\" offset=\"" << offset
958 << "\"/>" << endl;
959 offset += bndryIds.m_memsize + sizeof(MUint);
960 /*
961 ofl << "<DataArray type=\"" << dataType64 << "\" Name=\"drho_dt\" format=\"appended\"
962 offset=\""<< offset <<"\"/>" << endl; offset += drho_dt.m_memsize + sizeof( MUint );
963 */
964 if(m_solver->m_levelSet || m_solver->m_levelSetMb) {
965 ofl << "<DataArray type=\"" << dataType << "\" Name=\"levelSetFunction\" format=\"appended\" offset=\""
966 << offset << "\"/>" << endl;
967 offset += levelSetFunction.m_memsize + sizeof(MUint);
968 }
969 if(writePointData) {
970 ofl << "</CellData>" << endl;
971 ofl << "<PointData Scalars=\"scalars\">" << endl;
972 }
973
974 ofl << "<DataArray type=\"" << dataType << "\" Name=\"pressure\" format=\"appended\" offset=\"" << offset
975 << "\"/>" << endl;
976 offset += pressure.m_memsize + sizeof(MUint);
977
978 ofl << "<DataArray type=\"" << dataType
979 << "\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"appended\" offset=\"" << offset << "\"/>" << endl;
980 offset += velocity.m_memsize + sizeof(MUint);
981
982 ofl << "<DataArray type=\"" << dataType << "\" Name=\"density\" format=\"appended\" offset=\"" << offset << "\"/>"
983 << endl;
984 offset += density.m_memsize + sizeof(MUint);
985
986 ofl << "<DataArray type=\"" << dataType << "\" Name=\"vorticity\" format=\"appended\" offset=\"" << offset
987 << "\"/>" << endl;
988 offset += vorticity.m_memsize + sizeof(MUint);
989
990 // ofl << "<DataArray type=\"" << dataType << "\" Name=\"Qcriterion\" format=\"appended\"
991 // offset=\""<< offset <<"\"/>" << endl;
992 // offset += Qcriterion.m_memsize + sizeof( MUint );
993
994 if(m_solver->m_combustion) {
995 ofl << "<DataArray type=\"" << dataType << "\" Name=\"progressVariable\" format=\"appended\" offset=\""
996 << offset << "\"/>" << endl;
997 offset += progressVariable.m_memsize + sizeof(MUint);
998 }
999 if(writePointData) {
1000 ofl << "</PointData>" << endl;
1001 } else {
1002 ofl << "</CellData>" << endl;
1003 }
1004 //================== /CellData/PointData =================
1005
1006
1007 ofl << "</Piece>" << endl;
1008 ofl << "</UnstructuredGrid>" << endl;
1009
1010
1011 //================== AppendedData =================
1012 ofl << "<AppendedData encoding=\"raw\">" << endl;
1013 ofl << "_";
1014 ofl.close();
1015 ofl.open(fileName, ios_base::out | ios_base::app | ios_base::binary);
1016
1017 //----------
1018 // point coordinates
1019 uinumber = (MUint)points.m_memsize;
1020 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(MUint));
1021 ofl.write(reinterpret_cast<const char*>(points.getPointer()), uinumber);
1022
1023 //----------
1024 // connectivity
1025 uinumber = (MUint)connectivity.m_memsize;
1026 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(MUint));
1027 ofl.write(reinterpret_cast<const char*>(connectivity.getPointer()), uinumber);
1028
1029 //----------
1030 // offsets
1031 uinumber = (MUint)offsets.m_memsize;
1032 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1033 ofl.write(reinterpret_cast<const char*>(offsets.getPointer()), uinumber);
1034
1035
1036 //----------
1037 // cell types
1038 uinumber = (MUint)types.m_memsize;
1039 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1040 ofl.write(reinterpret_cast<const char*>(types.getPointer()), uinumber);
1041
1042 //----------
1043 // cellIds
1044 uinumber = (MUint)cellIds.m_memsize;
1045 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1046 ofl.write(reinterpret_cast<const char*>(cellIds.getPointer()), uinumber);
1047
1048 //----------
1049 // bndryIds
1050 uinumber = (MUint)bndryIds.m_memsize;
1051 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1052 ofl.write(reinterpret_cast<const char*>(bndryIds.getPointer()), uinumber);
1053 /*
1054 //----------
1055 // drho_dt
1056 uinumber = (MUint)drho_dt.m_memsize;
1057 ofl.write(reinterpret_cast<const char *> (&uinumber), sizeof(uinumber));
1058 ofl.write(reinterpret_cast<const char *> (drho_dt.getPointer()), uinumber);
1059 */
1060
1061 //----------
1062 // levelSetFunction
1063 if(m_solver->m_levelSet && (m_solver->m_levelSetMb)) {
1064 uinumber = (MUint)levelSetFunction.m_memsize;
1065 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1066 ofl.write(reinterpret_cast<const char*>(levelSetFunction.getPointer()), uinumber);
1067 }
1068 //----------
1069 // pressure
1070 uinumber = (MUint)pressure.m_memsize;
1071 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1072 ofl.write(reinterpret_cast<const char*>(pressure.getPointer()), uinumber);
1073
1074 //----------
1075 // velocity
1076 uinumber = (MUint)velocity.m_memsize;
1077 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1078 ofl.write(reinterpret_cast<const char*>(velocity.getPointer()), uinumber);
1079
1080 //----------
1081 // density
1082 uinumber = (MUint)density.m_memsize;
1083 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1084 ofl.write(reinterpret_cast<const char*>(density.getPointer()), uinumber);
1085
1086 //----------
1087 // vorticity
1088 uinumber = (MUint)vorticity.m_memsize;
1089 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1090 ofl.write(reinterpret_cast<const char*>(vorticity.getPointer()), uinumber);
1091 /*
1092 //----------
1093 // Qcriterion
1094 uinumber = (MUint)Qcriterion.m_memsize;
1095 ofl.write(reinterpret_cast<const char *> (&uinumber), sizeof(uinumber));
1096 ofl.write(reinterpret_cast<const char *> (Qcriterion.getPointer()), uinumber);
1097 */
1098
1099 //----------
1100 // progress variable
1101 if(m_solver->m_combustion) {
1102 uinumber = (MUint)progressVariable.m_memsize;
1103 ofl.write(reinterpret_cast<const char*>(&uinumber), sizeof(uinumber));
1104 ofl.write(reinterpret_cast<const char*>(progressVariable.getPointer()), uinumber);
1105 }
1106 //----------
1107 ofl.close();
1108 ofl.open(fileName, ios_base::app);
1109 ofl << endl;
1110 ofl << "</AppendedData>" << endl;
1111 //================== /AppendedData =================
1112
1113
1114 ofl << "</VTKFile>" << endl;
1115 ofl.close();
1116 //================== /VTKFile =================
1117 } else {
1118 cerr << "ERROR! COULD NOT OPEN FILE " << fileName << " for writing! " << endl;
1119 }
1120
1121 for(MInt i = 0; i < m_solver->m_bndryCells->size(); i++) {
1122 delete[] pointInside[i];
1123 }
1124 delete[] pointInside;
1125 pointInside = nullptr;
1126 }
1127 else IF_CONSTEXPR(nDim == 3) {
1128 static_cast<void>(fileName); // silence unused-parameter warning for 3D
1129 MString fileName2 = m_solver->m_solutionOutput + "QOUT_" + to_string(globalTimeStep) + ".vtu";
1130 if(m_solver->m_multipleFvSolver) {
1131 fileName2 = m_solver->m_solutionOutput + "QOUT" + to_string(m_solver->m_solverId) + "_"
1132 + to_string(globalTimeStep) + ".vtu";
1133 }
1134 MString geomFileName = m_solver->m_solutionOutput + "GEOM_" + to_string(globalTimeStep) + ".vtp";
1135 writeVtuOutputParallel<>(fileName2, geomFileName);
1136 }
1137}
FvBndryCndXD::m_solver
FvCartesianSolverXD< nDim, SysEqn > * m_solver
Definition: fvcartesianbndrycndxd.h:138
FvCartesianSolverXD::a_levelSetFunction
MFloat & a_levelSetFunction(const MInt cellId, const MInt set)
Returns the levelSet-value for fv-CellId cellId and set.
Definition: fvcartesiansolverxd.h:999
FvCartesianSolverXD::m_fvBndryCnd
FvBndryCndXD< nDim_, SysEqn > * m_fvBndryCnd
Definition: fvcartesiansolverxd.h:313
FvCartesianSolverXD::c_cellLengthAtLevel
MFloat c_cellLengthAtLevel(const MInt level) const
Returns the length of the cell for level.
Definition: fvcartesiansolverxd.h:582
FvCartesianSolverXD::a_slope
MFloat & a_slope(const MInt cellId, MInt const varId, const MInt dir) override
Returns the slope of the cell cellId for the variable varId in direction dir.
Definition: fvcartesiansolverxd.h:765
FvCartesianSolverXD::m_geometry
Geometry< nDim > * m_geometry
Definition: fvcartesiansolverxd.h:205
FvCartesianSolverXD::a_level
MInt & a_level(const MInt cellId)
Returns the level of the cell from the fvcellcollector cellId.
Definition: fvcartesiansolverxd.h:431
FvCartesianSolverXD::timeStep
MFloat timeStep(MBool canSolver=false) noexcept
Definition: fvcartesiansolverxd.h:2503
FvCartesianSolverXD::a_pvariable
MFloat & a_pvariable(const MInt cellId, const MInt varId)
Returns primitive variable v of the cell cellId for variables varId.
Definition: fvcartesiansolverxd.h:616
Solver::outputDir
MString outputDir() const
Return the directory for output files.
Definition: solver.h:407
Solver::m_restart
MBool m_restart
Definition: solver.h:97
Solver::m_solverId
const MInt m_solverId
a unique solver identifier
Definition: solver.h:90
Solver::m_solutionInterval
MInt m_solutionInterval
The number of timesteps before writing the next solution file.
Definition: solver.h:74
Solver::m_Re
MFloat m_Re
the Reynolds number
Definition: solver.h:68
Solver::m_Ma
MFloat m_Ma
the Mach number
Definition: solver.h:71
VtkIo::initializeVtkXmlOutput
MBool initializeVtkXmlOutput(const MChar *, MString, MInt, MBool, MBool)
Definition: iovtk.cpp:29
VtkIo::sysEqn
SysEqn sysEqn() const
Definition: iovtk.h:43
hasPV_C
Checks if the primitive variable C exists.
Definition: fvcartesiansyseqntraits.h:38
maia::CartesianSolver::maxRefinementLevel
MInt maxRefinementLevel() const
Definition: cartesiansolver.h:68

◆ writeVtuArrayParallel()

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::writeVtuArrayParallel ( MPI_File &  file,
void *  base,
MPI_Offset  offset,
MPI_Offset  count,
MPI_Offset  globalCount 
)

Definition at line 134 of file iovtk.cpp.

135 {
136 MPI_File_seek(file, offset, MPI_SEEK_CUR, AT_);
137#ifdef VTU_OUTPUT_NONCOLLECTIVE
138 MPI_Status status;
139 MInt rcode = (count > 0) ? MPI_File_write(file, base, count, MPI_CHAR, &status) : MPI_SUCCESS;
140#else
141#ifndef VTU_OUTPUT_BLOCKING
142 MPI_Status status;
143 MInt rcode = MPI_File_write_all(file, base, count, MPI_CHAR, &status);
144#else
145 MPI_Request writeRequest = MPI_REQUEST_NULL;
146 MInt rcode = MPI_File_iwrite_all(file, base, count, MPI_CHAR, &writeRequest);
147 MPI_Wait(&writeRequest, MPI_STATUS_IGNORE, AT_);
148#endif
149#endif
150 MPI_File_seek(file, globalCount - offset - count, MPI_SEEK_CUR, AT_);
151 return rcode;
152}
MPI_File_seek
int MPI_File_seek(MPI_File mpi_fh, MPI_Offset offset, int whence, const MString &name)
same as MPI_File_seek
Definition: mpioverride.cpp:2239
MPI_Wait
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
Definition: mpioverride.cpp:732

◆ writeVtuOutputParallel()

template<MInt nDim, class SysEqn >
template<typename uint_t >
void VtkIo< nDim, SysEqn >::writeVtuOutputParallel ( const MString  fileName,
const MString  geomFileName,
const MInt  noSolverSpecificVars = 0,
const MFloatScratchSpace &  solverSpecificVars = MFloatScratchSpace(1, "writeVtuOutputParallel",                                                                        "defaultParameter1"),
const MInt  noUserVars = 0,
const MFloatScratchSpace &  userVars = MFloatScratchSpace(1, "writeVtuOutputParallel", "defaultParameter2") 
)

Definition at line 2818 of file iovtk.cpp.

2821 {
2822 TRACE();
2823
2824 // NOTE: currently VTU output is not working on claix
2825
2826 static_assert(std::is_same<uint_t, uint32_t>::value || std::is_same<uint_t, uint64_t>::value,
2827 "Invalid template type specified.");
2828 static_assert(sizeof(float) == 4, "VTU output will fail since float has an unexpected size on this architecture.");
2829 const char* const uIDataType = (std::is_same<uint_t, uint64_t>::value) ? "UInt64" : "UInt32";
2830 const char* const uI8DataType = "UInt8";
2831 const char* const uI64DataType = "UInt64";
2832 const char* const dataType = "Float32";
2833 const uint64_t noNodes = IPOW2(nDim);
2834 const uint8_t baseCellType = (nDim == 3) ? 11 : 8;
2835 const uint8_t polyCellType = 42;
2836
2837 NEW_TIMER_GROUP_STATIC(t_vtuTimer, "writeVtuOutputParallel");
2838 NEW_TIMER_STATIC(t_vtutotal, "VtuOutputParallel", t_vtuTimer);
2839 NEW_SUB_TIMER_STATIC(t_faces, "faces", t_vtutotal);
2840 NEW_SUB_TIMER_STATIC(t_init, "init", t_vtutotal);
2841 NEW_SUB_TIMER_STATIC(t_prepare, "prepare", t_vtutotal);
2842 NEW_SUB_TIMER_STATIC(t_cells, "cells", t_vtutotal);
2843 NEW_SUB_TIMER_STATIC(t_variables, "variables", t_vtutotal);
2844 NEW_SUB_TIMER_STATIC(t_geometry, "geometry", t_vtutotal);
2845 RECORD_TIMER_START(t_vtutotal);
2846 RECORD_TIMER_START(t_faces);
2847
2848 const auto noCells = (uint64_t)m_solver->m_extractedCells->size();
2849 auto noPoints = (uint64_t)m_solver->m_gridPoints->size();
2850 auto connSize = (uint64_t)(noNodes * noCells);
2851 auto facesSize = (uint64_t)0;
2852 vector<uint_t>* facestream = nullptr;
2853 vector<uint_t>* conn = nullptr;
2854 ScratchSpace<MInt> polyhedralCell(m_solver->a_noCells(), AT_, "polyhedralCell");
2855 MInt noPolyhedra = 0;
2856 polyhedralCell.fill(-1);
2857 const MBool getInternalPolyhedra = true;
2858 if(m_solver->m_vtuCutCellOutput) {
2859 connSize = vtuAssembleFaceStream(facestream, conn, facesSize, polyhedralCell, noPolyhedra, getInternalPolyhedra);
2860 if(facestream == nullptr || conn == nullptr) mTerm(1, AT_, "empty stream");
2861 }
2862 noPoints = (uint64_t)m_solver->m_gridPoints->size(); // update after cut cell assembly
2863 RECORD_TIMER_STOP(t_faces);
2864 RECORD_TIMER_START(t_init);
2865
2866
2867 const MBool mergeGridPoints = true;
2868 const uint64_t noPoints0 = noPoints;
2869 ScratchSpace<MInt> pointDomainId(noPoints, AT_, "pointDomainId");
2870 ScratchSpace<MInt> pointRemapId(noPoints, AT_, "pointRemapId");
2871 pointDomainId.fill(domainId());
2872 pointRemapId.fill(-1);
2873 if(mergeGridPoints) {
2874 MIntScratchSpace nghbrList(300, AT_, "nghbrList");
2875 MIntScratchSpace layerId(300, AT_, "layerId");
2876 ScratchSpace<MInt> cellMap(m_solver->a_noCells(), AT_, "cellMap");
2877 cellMap.fill(-1);
2878 MInt noSendDat0 = 0;
2879 MInt noRecvDat0 = 0;
2880 MInt noSendDatTotal = 0;
2881 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2882 noSendDat0 += m_solver->noWindowCells(i);
2883 noRecvDat0 += m_solver->noHaloCells(i);
2884 }
2885 ScratchSpace<MInt> noCellPoints(noSendDat0 + 1, AT_, "noCellPoints");
2886 ScratchSpace<MInt> noCellPointsRecv(noRecvDat0 + 1, AT_, "noCellPointsRecv");
2887 ScratchSpace<MInt> noSendDat(noDomains(), AT_, "noSendDat");
2888 ScratchSpace<MInt> noRecvDat(noDomains(), AT_, "noRecvDat");
2889 noCellPoints.fill(0);
2890 noSendDat.fill(0);
2891 noRecvDat.fill(0);
2892 for(MInt c = 0; c < (signed)noCells; c++) {
2893 cellMap(m_solver->m_extractedCells->a[c].m_cellId) = c;
2894 }
2895 MInt scnt0 = 0;
2896 noSendDatTotal = 0;
2897 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2898 noSendDat[i] = 0;
2899 for(MInt j = 0; j < m_solver->noWindowCells(i); j++) {
2900 MInt c = cellMap(m_solver->windowCellId(i, j));
2901 MInt pc = (c > -1) ? polyhedralCell[c] : -1;
2902 MInt np = 0;
2903 if(c > -1) np = (pc > -1) ? conn[pc].size() : noNodes;
2904 noCellPoints(scnt0) = np;
2905 noSendDatTotal += np;
2906 noSendDat[i] += np;
2907 scnt0++;
2908 }
2909 }
2910 ScratchSpace<MPI_Request> sendReq(m_solver->noNeighborDomains(), AT_, "sendReq");
2911 sendReq.fill(MPI_REQUEST_NULL);
2912 scnt0 = 0;
2913 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2914 MPI_Issend(&(noCellPoints[scnt0]), m_solver->noWindowCells(i), MPI_INT, m_solver->neighborDomain(i), 132,
2915 mpiComm(), &sendReq[i], AT_, "(noCellPoints[scnt0])");
2916 scnt0 += m_solver->noWindowCells(i);
2917 }
2918 MInt rcnt0 = 0;
2919 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2920 MPI_Recv(&(noCellPointsRecv[rcnt0]), m_solver->noHaloCells(i), MPI_INT, m_solver->neighborDomain(i), 132,
2921 mpiComm(), MPI_STATUS_IGNORE, AT_, "(noCellPointsRecv[rcnt0])");
2922 rcnt0 += m_solver->noHaloCells(i);
2923 }
2924 if(m_solver->noNeighborDomains() > 0)
2925 MPI_Waitall(m_solver->noNeighborDomains(), &sendReq[0], MPI_STATUSES_IGNORE, AT_);
2926 MInt noRecvDatTotal = 0;
2927 rcnt0 = 0;
2928 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2929 noRecvDat[i] = 0;
2930 for(MInt j = 0; j < m_solver->noHaloCells(i); j++) {
2931 noRecvDat[i] += noCellPointsRecv[rcnt0];
2932 noRecvDatTotal += noCellPointsRecv[rcnt0];
2933 rcnt0++;
2934 }
2935 }
2936 ScratchSpace<MInt> sendGridPoints(noSendDatTotal + 1, AT_, "sendGridPoints");
2937 ScratchSpace<MInt> recvGridPoints(noRecvDatTotal + 1, AT_, "recvGridPoints");
2938 ScratchSpace<MFloat> sendGridPointsCoord(noSendDatTotal + 1, 3, AT_, "sendGridPointsCoord");
2939 ScratchSpace<MFloat> recvGridPointsCoord(noRecvDatTotal + 1, 3, AT_, "recvGridPointsCoord");
2940 scnt0 = 0;
2941 MInt scnt = 0;
2942 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2943 for(MInt j = 0; j < m_solver->noWindowCells(i); j++) {
2944 MInt c = cellMap(m_solver->windowCellId(i, j));
2945 if(noCellPoints(scnt0) > 0) {
2946 if(c < 0) {
2947 cerr << "case should not occur" << endl;
2948 continue;
2949 }
2950 MInt pc = polyhedralCell[c];
2951 MInt noGridPoints = (pc > -1) ? conn[pc].size() : noNodes;
2952 for(MInt p = 0; p < noGridPoints; p++) {
2953 MInt gridPointId = (pc > -1) ? conn[pc][p] : m_solver->m_extractedCells->a[c].m_pointIds[p];
2954 if(gridPointId < 0) cerr << "negative grid point id " << endl;
2955 sendGridPoints(scnt) = gridPointId;
2956 for(MInt k = 0; k < nDim; k++) {
2957 sendGridPointsCoord(scnt, k) = m_solver->m_gridPoints->a[gridPointId].m_coordinates[k];
2958 }
2959 scnt++;
2960 }
2961 }
2962 scnt0++;
2963 }
2964 }
2965
2966 sendReq.fill(MPI_REQUEST_NULL);
2967 scnt = 0;
2968 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2969 MPI_Issend(&(sendGridPoints[scnt]), noSendDat[i], MPI_INT, m_solver->neighborDomain(i), 133, mpiComm(),
2970 &sendReq[i], AT_, "(sendGridPoints[scnt])");
2971 scnt += noSendDat[i];
2972 }
2973 MInt rcnt = 0;
2974 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2975 MPI_Recv(&(recvGridPoints[rcnt]), noRecvDat[i], MPI_INT, m_solver->neighborDomain(i), 133, mpiComm(),
2976 MPI_STATUS_IGNORE, AT_, "(recvGridPoints[rcnt])");
2977 rcnt += noRecvDat[i];
2978 }
2979 if(m_solver->noNeighborDomains() > 0)
2980 MPI_Waitall(m_solver->noNeighborDomains(), &sendReq[0], MPI_STATUSES_IGNORE, AT_);
2981
2982 sendReq.fill(MPI_REQUEST_NULL);
2983 scnt = 0;
2984 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2985 MPI_Issend(&(sendGridPointsCoord[scnt]), 3 * noSendDat[i], MPI_DOUBLE, m_solver->neighborDomain(i), 134,
2986 mpiComm(), &sendReq[i], AT_, "(sendGridPointsCoord[scnt])");
2987 scnt += 3 * noSendDat[i];
2988 }
2989 rcnt = 0;
2990 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
2991 MPI_Recv(&(recvGridPointsCoord[rcnt]), 3 * noRecvDat[i], MPI_DOUBLE, m_solver->neighborDomain(i), 134, mpiComm(),
2992 MPI_STATUS_IGNORE, AT_, "(recvGridPointsCoord[rcnt])");
2993 rcnt += 3 * noRecvDat[i];
2994 }
2995 if(m_solver->noNeighborDomains() > 0)
2996 MPI_Waitall(m_solver->noNeighborDomains(), &sendReq[0], MPI_STATUSES_IGNORE, AT_);
2997
2998 rcnt0 = 0;
2999 rcnt = 0;
3000 map<MInt, MInt> haloMap;
3001 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
3002 for(MInt j = 0; j < m_solver->noHaloCells(i); j++) {
3003 MInt haloId = m_solver->haloCellId(i, j);
3004 MInt noRecvPoints = noCellPointsRecv[rcnt0];
3005 if(noRecvPoints > 0) {
3006 // const MInt counter = m_solver->getAdjacentLeafCells<2>( haloId, 2, nghbrList, layerId );
3007 const MInt counter = m_solver->getAdjacentLeafCells_d2(haloId, 2, nghbrList, layerId);
3008 for(MInt n = 0; n < counter; n++) {
3009 MInt nghbrId = nghbrList[n];
3010 if(nghbrId < 0) continue;
3011 if(m_solver->a_isHalo(nghbrId)) continue;
3012 MInt ncell = cellMap[nghbrId];
3013 if(ncell < 0) continue;
3014 MInt pc = polyhedralCell[ncell];
3015 MFloat delta0 = (pc > -1) ? 1e-12 * m_solver->c_cellLengthAtCell(nghbrId)
3016 : 0.0001 * m_solver->c_cellLengthAtCell(nghbrId);
3017 MInt noPointsNghbr = (pc > -1) ? conn[pc].size() : noNodes;
3018 for(MInt p = 0; p < noRecvPoints; p++) {
3019 for(MInt q = 0; q < noPointsNghbr; q++) {
3020 MInt gridPointId = (pc > -1) ? conn[pc][q] : m_solver->m_extractedCells->a[ncell].m_pointIds[q];
3021 if(gridPointId < 0) cerr << "negative grid point id " << endl;
3022 MFloat delta = sqrt(
3023 POW2(m_solver->m_gridPoints->a[gridPointId].m_coordinates[0] - recvGridPointsCoord(rcnt + p, 0))
3024 + POW2(m_solver->m_gridPoints->a[gridPointId].m_coordinates[1] - recvGridPointsCoord(rcnt + p, 1))
3025 + POW2(m_solver->m_gridPoints->a[gridPointId].m_coordinates[2] - recvGridPointsCoord(rcnt + p, 2)));
3026 if(delta < delta0) {
3027 if(m_solver->neighborDomain(i) < pointDomainId[gridPointId]) {
3028 pointDomainId[gridPointId] = m_solver->neighborDomain(i);
3029 pointRemapId[gridPointId] = recvGridPoints[rcnt + p];
3030 haloMap[gridPointId] = rcnt + p;
3031 }
3032 }
3033 }
3034 }
3035 }
3036 }
3037 rcnt += noRecvPoints;
3038 rcnt0++;
3039 }
3040 }
3041
3042 // update number of points
3043 noPoints = 0;
3044 for(uint64_t p = 0; p < noPoints0; p++) {
3045 if(pointRemapId[p] < 0) {
3046 pointRemapId[p] = noPoints;
3047 if(pointDomainId[p] != domainId()) cerr << "unexpected case" << endl;
3048 noPoints++;
3049 }
3050 if(pointRemapId[p] < 0) {
3051 cerr << "unexpected case 2" << endl;
3052 }
3053 }
3054
3055
3056 // communicate point shifts
3057 scnt0 = 0;
3058 scnt = 0;
3059 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
3060 for(MInt j = 0; j < m_solver->noWindowCells(i); j++) {
3061 MInt c = cellMap(m_solver->windowCellId(i, j));
3062 if(noCellPoints(scnt0) > 0) {
3063 if(c < 0) {
3064 cerr << "case should not occur" << endl;
3065 continue;
3066 }
3067 MInt pc = polyhedralCell[c];
3068 MInt noGridPoints = (pc > -1) ? conn[pc].size() : noNodes;
3069 for(MInt p = 0; p < noGridPoints; p++) {
3070 MInt gridPointId = (pc > -1) ? conn[pc][p] : m_solver->m_extractedCells->a[c].m_pointIds[p];
3071 if(gridPointId < 0) cerr << "negative grid point id " << endl;
3072 sendGridPoints(scnt) = (pointDomainId[gridPointId] == domainId()) ? pointRemapId[gridPointId] : -1;
3073 scnt++;
3074 }
3075 }
3076 scnt0++;
3077 }
3078 }
3079 sendReq.fill(MPI_REQUEST_NULL);
3080 scnt = 0;
3081 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
3082 MPI_Issend(&(sendGridPoints[scnt]), noSendDat[i], MPI_INT, m_solver->neighborDomain(i), 135, mpiComm(),
3083 &sendReq[i], AT_, "(sendGridPoints[scnt])");
3084 scnt += noSendDat[i];
3085 }
3086 rcnt = 0;
3087 for(MInt i = 0; i < m_solver->noNeighborDomains(); i++) {
3088 MPI_Recv(&(recvGridPoints[rcnt]), noRecvDat[i], MPI_INT, m_solver->neighborDomain(i), 135, mpiComm(),
3089 MPI_STATUS_IGNORE, AT_, "(recvGridPoints[rcnt])");
3090 rcnt += noRecvDat[i];
3091 }
3092 if(m_solver->noNeighborDomains() > 0)
3093 MPI_Waitall(m_solver->noNeighborDomains(), &sendReq[0], MPI_STATUSES_IGNORE, AT_);
3094
3095 for(auto& it : haloMap) {
3096 MInt gridPointId = it.first;
3097 MInt otherId = recvGridPoints[it.second];
3098 if(gridPointId < 0 || gridPointId >= (signed)noPoints0) cerr << "point out of range" << endl;
3099 pointRemapId[gridPointId] = otherId;
3100 if(otherId < 0 || (pointDomainId[gridPointId] >= domainId())) {
3101 cerr << "inconsistent point shifts " << gridPointId << " " << otherId << " " << pointDomainId[gridPointId]
3102 << " " << domainId() << endl;
3103 continue;
3104 }
3105 }
3106
3107
3108 } else {
3109 for(uint64_t p = 0; p < noPoints; p++) {
3110 pointRemapId[p] = p;
3111 }
3112 }
3113
3114 ScratchSpace<uint64_t> noPointsPerDomain(noDomains(), AT_, "noPointsPerDomain");
3115 ScratchSpace<uint64_t> pointOffsetsGlobal(noDomains() + 1, AT_, "pointOffsetsGlobal");
3116 ScratchSpace<uint64_t> noCellsPerDomain(noDomains(), AT_, "noCellsPerDomain");
3117 ScratchSpace<uint64_t> connSizePerDomain(noDomains(), AT_, "connSizePerDomain");
3118 ScratchSpace<uint64_t> facesPerDomain(noDomains(), AT_, "facesPerDomain");
3119 ScratchSpace<uint64_t> dataPerDomain(noDomains(), 4, AT_, "dataPerDomain");
3120 dataPerDomain(domainId(), 0) = noPoints;
3121 dataPerDomain(domainId(), 1) = noCells;
3122 dataPerDomain(domainId(), 2) = connSize;
3123 dataPerDomain(domainId(), 3) = facesSize;
3124 MPI_Allgather(MPI_IN_PLACE, 4, MPI_UINT64_T, &dataPerDomain[0], 4, MPI_UINT64_T, mpiComm(), AT_, "MPI_IN_PLACE",
3125 "dataPerDomain[0]");
3126 pointOffsetsGlobal[0] = 0;
3127 for(MInt d = 0; d < noDomains(); d++) {
3128 pointOffsetsGlobal[d + 1] = pointOffsetsGlobal[d] + dataPerDomain(d, 0);
3129 noPointsPerDomain[d] = dataPerDomain(d, 0);
3130 noCellsPerDomain[d] = dataPerDomain(d, 1);
3131 connSizePerDomain[d] = dataPerDomain(d, 2);
3132 facesPerDomain[d] = dataPerDomain(d, 3);
3133 }
3134
3135 uint64_t globalNoPoints = 0;
3136 uint64_t globalNoCells = 0;
3137 uint64_t globalConnSize = 0;
3138 uint64_t globalFacesSize = 0;
3139 uint64_t globalPointOffset = 0;
3140 uint64_t globalCellOffset = 0;
3141 uint64_t globalConnSizeOffset = 0;
3142 uint64_t globalFacesOffset = 0;
3143 uint64_t offset = 0;
3144 for(MInt d = 0; d < noDomains(); d++) {
3145 globalNoPoints += noPointsPerDomain[d];
3146 globalNoCells += noCellsPerDomain[d];
3147 globalConnSize += connSizePerDomain[d];
3148 globalFacesSize += facesPerDomain[d];
3149 }
3150 for(MInt d = 0; d < domainId(); d++) {
3151 globalPointOffset += noPointsPerDomain[d];
3152 globalCellOffset += noCellsPerDomain[d];
3153 globalConnSizeOffset += connSizePerDomain[d];
3154 globalFacesOffset += facesPerDomain[d];
3155 }
3156
3157 // if output size exceeds 32bit boundary rerun with 64 bit data types
3158 if(std::is_same<uint_t, uint32_t>::value) {
3159 if(globalNoCells > (uint64_t)std::numeric_limits<uint_t>::max()
3160 || globalNoPoints > (uint64_t)std::numeric_limits<uint_t>::max()) {
3161 return writeVtuOutputParallel<uint64_t>(fileName, geomFileName);
3162 }
3163 }
3164
3165 RECORD_TIMER_STOP(t_init);
3166 RECORD_TIMER_START(t_prepare);
3167
3168 if(string2enum(m_solver->m_outputFormat) == VTU) {
3169 //-----------
3170 // points
3171 uint64_t pointsMemsize = 3 * noPoints * sizeof(float);
3172 uint64_t pointsMemsizeGlobal = 3 * globalNoPoints * sizeof(float);
3173 uint64_t pointsOffset = 3 * globalPointOffset * sizeof(float);
3174 const MInt pointsMaxSize = estimateMemorySizeSolverwise(noPoints, noPointsPerDomain, 3 * sizeof(float));
3175 ScratchSpace<float> points(pointsMaxSize, 3, AT_, "points");
3176 uint_t cnt = 0;
3177 for(uint_t p = 0; p < noPoints0; p++) {
3178 if(pointDomainId[p] != domainId()) continue;
3179 for(MInt i = 0; i < nDim; i++) {
3180 points((MInt)cnt, i) = (float)m_solver->m_gridPoints->a[p].m_coordinates[i];
3181 }
3182 IF_CONSTEXPR(nDim == 2) points((MInt)cnt, 2) = (float)0.0;
3183 cnt++;
3184 }
3185 insertDataHeader(reinterpret_cast<char*>(&points(0)), pointsMemsize, pointsMemsizeGlobal, pointsOffset);
3186
3187 //-----------
3188 // connectivity
3189 uint64_t connectivityMemsize = connSize * sizeof(uint_t);
3190 uint64_t connectivityOffset = globalConnSizeOffset * sizeof(uint_t);
3191 uint64_t connectivityMemsizeGlobal = globalConnSize * sizeof(uint_t);
3192 const MInt connectivityMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, noNodes * sizeof(uint_t));
3193 ScratchSpace<uint_t> connectivity(mMax(connectivityMaxSize * noNodes, connSize + 20), AT_, "connectivity");
3194 cnt = 0;
3195 for(MInt c = 0; c < (signed)noCells; c++) {
3196 MInt pc = polyhedralCell[c];
3197 if(pc < 0) {
3198 for(MInt p = 0; p < (signed)noNodes; p++) {
3199 if(m_solver->m_extractedCells->a[c].m_pointIds[p] < 0)
3200 cerr << domainId() << ": warning negative point " << c << " " << m_solver->m_extractedCells->a[c].m_cellId
3201 << " " << m_solver->m_extractedCells->a[c].m_pointIds[p] << endl;
3202 MInt gridPointId = m_solver->m_extractedCells->a[c].m_pointIds[p];
3203 connectivity(cnt) =
3204 (uint_t)pointOffsetsGlobal[pointDomainId[gridPointId]] + (uint_t)pointRemapId[gridPointId];
3205 cnt++;
3206 }
3207 } else {
3208 ASSERT(cnt + conn[pc].size() <= connectivity.size(), "");
3209 for(MUint i = 0; i < conn[pc].size(); i++) {
3210 uint_t gridPointId = conn[pc][i];
3211 connectivity(cnt) =
3212 (uint_t)pointOffsetsGlobal[pointDomainId[gridPointId]] + (uint_t)pointRemapId[gridPointId];
3213 cnt++;
3214 }
3215 }
3216 }
3217 for(MUint i = 0; i < connSize; i++) {
3218 if(connectivity[i] >= globalNoPoints) {
3219 cerr << domainId() << ": invalid grid point " << i << " " << connectivity[i] << " " << globalPointOffset << " "
3220 << noPoints << " " << globalNoPoints << " " << globalTimeStep << endl;
3221 connectivity[i] = 0;
3222 }
3223 }
3224 insertDataHeader(reinterpret_cast<char*>(&connectivity(0)), connectivityMemsize, connectivityMemsizeGlobal,
3225 connectivityOffset);
3226
3227 //-----------
3228 // offsets
3229 uint64_t offsetsMemsize = noCells * sizeof(uint64_t);
3230 uint64_t offsetsOffset = globalCellOffset * sizeof(uint64_t);
3231 uint64_t offsetsMemsizeGlobal = globalNoCells * sizeof(uint64_t);
3232 const MInt offsetsMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, sizeof(uint64_t));
3233 ScratchSpace<uint64_t> offsets(offsetsMaxSize, AT_, "offsets");
3234 uint64_t offs = globalConnSizeOffset;
3235 for(uint64_t c = 0; c < noCells; c++) {
3236 MInt pc = polyhedralCell[c];
3237 if(pc < 0)
3238 offs += (uint64_t)noNodes;
3239 else
3240 offs += (uint64_t)conn[pc].size();
3241 offsets(c) = offs;
3242 if(offsets(c) - globalConnSizeOffset > connSize || offsets(c) > globalConnSize) {
3243 cerr << domainId() << ": invalid offset: " << offsets(c) << " " << connSize << " " << globalConnSizeOffset
3244 << " " << globalConnSize << endl;
3245 }
3246 }
3247 insertDataHeader(reinterpret_cast<char*>(&offsets(0)), offsetsMemsize, offsetsMemsizeGlobal, offsetsOffset);
3248
3249 //-----------
3250 // types
3251 uint64_t typesMemsize = noCells * sizeof(uint8_t);
3252 uint64_t typesOffset = globalCellOffset * sizeof(uint8_t);
3253 uint64_t typesMemsizeGlobal = globalNoCells * sizeof(uint8_t);
3254 const MInt typesMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, sizeof(uint8_t));
3255 ScratchSpace<uint8_t> types(typesMaxSize, AT_, "types");
3256 for(MInt c = 0; c < (signed)noCells; c++) {
3257 MInt pc = polyhedralCell[c];
3258 types(c) = (pc < 0) ? baseCellType : polyCellType;
3259 }
3260 insertDataHeader(reinterpret_cast<char*>(&types(0)), typesMemsize, typesMemsizeGlobal, typesOffset);
3261
3262 //-----------
3263 // faces
3264 uint64_t facesMemsize = facesSize * sizeof(uint_t);
3265 uint64_t facesOffset = globalFacesOffset * sizeof(uint_t);
3266 uint64_t facesMemsizeGlobal = globalFacesSize * sizeof(uint_t);
3267 const MInt facesMaxSize = estimateMemorySizeSolverwise(facesSize, facesPerDomain, sizeof(uint_t));
3268 ScratchSpace<uint_t> faces(facesMaxSize, AT_, "faces");
3269 cnt = 0;
3270 for(MInt c = 0; c < (signed)noCells; c++) {
3271 MInt pc = polyhedralCell[c];
3272 if(pc < 0) {
3273 faces(cnt) = 0; // can skip?
3274 cnt++; // can skip?
3275 } else {
3276 if(facestream[pc].empty()) mTerm(1, AT_, "empty face stream");
3277 if(conn[pc].empty()) mTerm(1, AT_, "empty conn stream");
3278 // copy( facestream[pc].begin(), facestream[pc].end(), &faces(cnt) );
3279 MInt fcnt = 0;
3280 const MInt noFaces = facestream[pc][fcnt];
3281 faces(cnt) = facestream[pc][fcnt];
3282 cnt++;
3283 fcnt++;
3284 for(MInt face = 0; face < noFaces; face++) {
3285 MInt pointCnt = facestream[pc][fcnt];
3286 faces(cnt) = facestream[pc][fcnt];
3287 cnt++;
3288 fcnt++;
3289 for(MInt p = 0; p < pointCnt; p++) {
3290 uint_t gridPointId = facestream[pc][fcnt];
3291 faces(cnt) = (uint_t)pointOffsetsGlobal[pointDomainId[gridPointId]] + (uint_t)pointRemapId[gridPointId];
3292 cnt++;
3293 fcnt++;
3294 }
3295 }
3296 }
3297 }
3298
3299 insertDataHeader(reinterpret_cast<char*>(&faces(0)), facesMemsize, facesMemsizeGlobal, facesOffset);
3300
3301 //-----------
3302 // faceoffsets
3303 uint64_t faceoffsetsMemsize = noCells * sizeof(uint64_t);
3304 uint64_t faceoffsetsOffset = globalCellOffset * sizeof(uint64_t);
3305 uint64_t faceoffsetsMemsizeGlobal = globalNoCells * sizeof(uint64_t);
3306 const MInt faceoffsetsMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, sizeof(uint64_t));
3307 ScratchSpace<uint64_t> faceoffsets(faceoffsetsMaxSize, AT_, "faceoffsets");
3308 offs = globalFacesOffset;
3309 for(MInt c = 0; c < (signed)noCells; c++) {
3310 MInt pc = polyhedralCell[c];
3311 if(pc < 0)
3312 offs++;
3313 else
3314 offs += (uint64_t)facestream[pc].size();
3315 faceoffsets(c) = offs;
3316 }
3317 insertDataHeader(reinterpret_cast<char*>(&faceoffsets(0)), faceoffsetsMemsize, faceoffsetsMemsizeGlobal,
3318 faceoffsetsOffset);
3319
3320 //-----------
3321 // globalId
3322 uint64_t globalIdMemsize = m_solver->m_vtuGlobalIdOutput ? noCells * sizeof(uint_t) : 0;
3323 uint64_t globalIdOffset = globalCellOffset * sizeof(uint_t);
3324 uint64_t globalIdMemsizeGlobal = globalNoCells * sizeof(uint_t);
3325 const MInt globalIdMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, sizeof(uint_t));
3326 ScratchSpace<uint_t> globalId(globalIdMaxSize, AT_, "globalId");
3327 if(m_solver->m_vtuGlobalIdOutput) {
3328 for(uint_t c = 0; c < noCells; c++) {
3329 globalId(c) = (uint_t)m_solver->c_globalId(
3330 m_solver->a_hasProperty(m_solver->m_extractedCells->a[c].m_cellId, SolverCell::IsSplitChild)
3331 ? m_solver->getAssociatedInternalCell(m_solver->m_extractedCells->a[c].m_cellId)
3332 : m_solver->m_extractedCells->a[c].m_cellId);
3333 }
3334 insertDataHeader(reinterpret_cast<char*>(&globalId(0)), globalIdMemsize, globalIdMemsizeGlobal, globalIdOffset);
3335 }
3336
3337 //-----------
3338 // domainIds
3339 uint64_t domainIdsMemsize = m_solver->m_vtuDomainIdOutput ? noCells * sizeof(uint_t) : 0;
3340 uint64_t domainIdsOffset = globalCellOffset * sizeof(uint_t);
3341 uint64_t domainIdsMemsizeGlobal = globalNoCells * sizeof(uint_t);
3342 const MInt domainIdsMaxSize = estimateMemorySizeSolverwise(noCells, noCellsPerDomain, sizeof(uint_t));
3343 ScratchSpace<uint_t> domainIds(domainIdsMaxSize, AT_, "domainIds");
3344 if(m_solver->m_vtuGlobalIdOutput) {
3345 for(uint_t c = 0; c < noCells; c++) {
3346 domainIds(c) = (uint_t)domainId();
3347 }
3348
3349 insertDataHeader(reinterpret_cast<char*>(&domainIds(0)), domainIdsMemsize, domainIdsMemsizeGlobal,
3350 domainIdsOffset);
3351 }
3352
3353 //-----------
3354 // variables
3355 const uint64_t dataSize = m_solver->m_vtuWritePointData ? noPoints : noCells;
3356 const uint64_t globalDataSize = m_solver->m_vtuWritePointData ? globalNoPoints : globalNoCells;
3357 const uint64_t globalDataOffset = m_solver->m_vtuWritePointData ? globalPointOffset : globalCellOffset;
3358 ScratchSpace<uint64_t>* noDataPerDomain = m_solver->m_vtuWritePointData ? &noPointsPerDomain : &noCellsPerDomain;
3359 uint64_t velocityMemsize = 3 * dataSize * sizeof(float);
3360 uint64_t vorticityMemsize = m_solver->m_vtuVorticityOutput ? 3 * dataSize * sizeof(float) : 0;
3361 uint64_t velGradMemsize = m_solver->m_vtuVelocityGradientOutput ? 9 * dataSize * sizeof(float) : 0;
3362 uint64_t pressureMemsize = dataSize * sizeof(float);
3363 uint64_t densityMemsize = m_solver->m_vtuDensityOutput ? dataSize * sizeof(float) : 0;
3364 uint64_t progressMemsize = dataSize * m_solver->m_noSpecies * sizeof(float);
3365 uint64_t levelSetMemsize = m_solver->m_vtuLevelSetOutput ? dataSize * sizeof(float) : 0;
3366 uint64_t QMemsize = m_solver->m_vtuQCriterionOutput ? dataSize * sizeof(float) : 0;
3367 uint64_t L2Memsize = m_solver->m_vtuLambda2Output ? dataSize * sizeof(float) : 0;
3368 uint64_t solverSpecificVarsMemsize = noSolverSpecificVars > 0 ? dataSize * sizeof(float) * noSolverSpecificVars : 0;
3369 uint64_t userVarsMemsize = noUserVars > 0 ? dataSize * sizeof(float) * noUserVars : 0;
3370 uint64_t velocityOffset = 3 * globalDataOffset * sizeof(float);
3371 uint64_t vorticityOffset = m_solver->m_vtuVorticityOutput ? 3 * globalDataOffset * sizeof(float) : 0;
3372 uint64_t velGradOffset = m_solver->m_vtuVelocityGradientOutput ? 9 * globalDataOffset * sizeof(float) : 0;
3373 uint64_t pressureOffset = globalDataOffset * sizeof(float);
3374 uint64_t densityOffset = m_solver->m_vtuDensityOutput ? globalDataOffset * sizeof(float) : 0;
3375 uint64_t progressOffset = globalDataOffset * m_solver->m_noSpecies * sizeof(float);
3376 uint64_t levelSetOffset = m_solver->m_vtuLevelSetOutput ? globalDataOffset * sizeof(float) : 0;
3377 uint64_t QOffset = m_solver->m_vtuQCriterionOutput ? globalDataOffset * sizeof(float) : 0;
3378 uint64_t L2Offset = m_solver->m_vtuLambda2Output ? globalDataOffset * sizeof(float) : 0;
3379 uint64_t solverSpecificVarsOffset =
3380 noSolverSpecificVars > 0 ? globalDataOffset * sizeof(float) * noSolverSpecificVars : 0;
3381 uint64_t userVarsOffset = noUserVars > 0 ? globalDataOffset * sizeof(float) * noUserVars : 0;
3382 uint64_t userVarsMemsizeGlobal = noUserVars > 0 ? globalDataSize * noUserVars * sizeof(float) : 0;
3383 uint64_t solverSpecificVarsMemsizeGlobal =
3384 noSolverSpecificVars > 0 ? globalDataSize * noSolverSpecificVars * sizeof(float) : 0;
3385 uint64_t velocityMemsizeGlobal = 3 * globalDataSize * sizeof(float);
3386 uint64_t vorticityMemsizeGlobal = m_solver->m_vtuVorticityOutput ? 3 * globalDataSize * sizeof(float) : 0;
3387 uint64_t velGradMemsizeGlobal = m_solver->m_vtuVelocityGradientOutput ? 9 * globalDataSize * sizeof(float) : 0;
3388 uint64_t pressureMemsizeGlobal = globalDataSize * sizeof(float);
3389 uint64_t densityMemsizeGlobal = m_solver->m_vtuDensityOutput ? globalDataSize * sizeof(float) : 0;
3390 uint64_t progressMemsizeGlobal = globalDataSize * m_solver->m_noSpecies * sizeof(float);
3391 uint64_t levelSetMemsizeGlobal = m_solver->m_vtuLevelSetOutput ? globalDataSize * sizeof(float) : 0;
3392 uint64_t QMemsizeGlobal = m_solver->m_vtuQCriterionOutput ? globalDataSize * sizeof(float) : 0;
3393 uint64_t L2MemsizeGlobal = m_solver->m_vtuLambda2Output ? globalDataSize * sizeof(float) : 0;
3394 const MInt vectorMaxSize = estimateMemorySizeSolverwise(dataSize, *noDataPerDomain, 3 * sizeof(float));
3395 const MInt scalarMaxSize = estimateMemorySizeSolverwise(dataSize, *noDataPerDomain, sizeof(float));
3396 ScratchSpace<float> velocity(vectorMaxSize, 3, AT_, "velocity");
3397 ScratchSpace<float> vorticity(m_solver->m_vtuVorticityOutput ? vectorMaxSize : 0, 3, AT_, "vorticity");
3398 ScratchSpace<float> velGrad(m_solver->m_vtuVelocityGradientOutput ? vectorMaxSize : 0, 9, AT_, "velGrad");
3399 ScratchSpace<float> pressure(scalarMaxSize, AT_, "pressure");
3400 ScratchSpace<float> density(m_solver->m_vtuDensityOutput ? scalarMaxSize : 0, AT_, "density");
3401 ScratchSpace<float> progress(scalarMaxSize * m_solver->m_noSpecies, AT_, "progress");
3402 ScratchSpace<float> levelSet(m_solver->m_vtuLevelSetOutput ? scalarMaxSize : 0, AT_, "levelSet");
3403 ScratchSpace<float> Q(m_solver->m_vtuQCriterionOutput ? scalarMaxSize : 0, AT_, "Q");
3404 ScratchSpace<float> lambda2(m_solver->m_vtuLambda2Output ? scalarMaxSize : 0, AT_, "lambda2");
3405 ScratchSpace<float> solverSpecificVarsOut(noSolverSpecificVars > 0 ? scalarMaxSize * noSolverSpecificVars : 0, AT_,
3406 "solverSpecificVars");
3407 ScratchSpace<float> userVarsOut(noUserVars > 0 ? scalarMaxSize * noUserVars : 0, AT_, "userVars");
3408 if(m_solver->m_vtuWritePointData) {
3409 for(MInt p = 0; p < (signed)noPoints; p++) {
3410 for(MInt i = 0; i < 3; i++)
3411 velocity(p, i) = (float)0.0;
3412 for(MInt i = 0; i < noUserVars; i++)
3413 userVarsOut[p * noUserVars + i] = (float)0.0;
3414 for(MInt i = 0; i < noSolverSpecificVars; i++)
3415 solverSpecificVarsOut[p * noSolverSpecificVars + i] = (float)0.0;
3416 if(m_solver->m_vtuVorticityOutput)
3417 for(MInt i = 0; i < nDim; i++)
3418 vorticity(p, i) = (float)0.0;
3419 if(m_solver->m_vtuVelocityGradientOutput)
3420 for(MInt i = 0; i < nDim * nDim; i++)
3421 velGrad(p, i) = (float)0.0;
3422 pressure(p) = (float)0.0;
3423 if(m_solver->m_vtuDensityOutput) density(p) = (float)0.0;
3424 if(m_solver->m_noSpecies > 0) progress(p) = (float)0.0;
3425 if(m_solver->m_vtuLevelSetOutput) levelSet(p) = (float)0.0;
3426 if(m_solver->m_vtuQCriterionOutput) Q(p) = (float)0.0;
3427 if(m_solver->m_vtuLambda2Output) lambda2(p) = (float)0.0;
3428 MFloat sum = F0;
3429 for(MInt n = 0; n < m_solver->m_gridPoints->a[p].m_noAdjacentCells; n++) {
3430 MInt cellId = m_solver->m_gridPoints->a[p].m_cellIds[n];
3431 MFloat dx = F0;
3432 for(MInt i = 0; i < nDim; i++) {
3433 dx += POW2(m_solver->a_coordinate(cellId, i) - m_solver->m_gridPoints->a[p].m_coordinates[i]);
3434 }
3435 dx = mMax(m_solver->m_eps, sqrt(dx));
3436 if(m_solver->m_vtuQCriterionOutput) {
3437 float q = F0;
3438 IF_CONSTEXPR(nDim == 3) {
3439 for(MInt i = 0; i < nDim; i++)
3440 for(MInt j = 0; j < nDim; j++)
3441 q += POW2(F1B2
3442 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3443 - m_solver->a_slope(cellId, sysEqn().PV->VV[j], i)))
3444 - POW2(F1B2
3445 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3446 + m_solver->a_slope(cellId, sysEqn().PV->VV[j], i)));
3447 q *= F1B2;
3448 }
3449 else {
3450 q += dx * F1B2
3451 * (m_solver->a_slope(cellId, sysEqn().PV->VV[1], 0)
3452 - m_solver->a_slope(cellId, sysEqn().PV->VV[0], 1));
3453 }
3454 Q(p) += dx * q;
3455 }
3456 for(MInt i = 0; i < nDim; i++)
3457 velocity(p, i) += dx * m_solver->a_pvariable(cellId, sysEqn().PV->VV[i]);
3458 if(m_solver->m_vtuVorticityOutput) {
3459 ASSERT(nDim == 3, "");
3460 vorticity(p, 0) +=
3461 dx * F1B2
3462 * (m_solver->a_slope(cellId, sysEqn().PV->VV[2], 1) - m_solver->a_slope(cellId, sysEqn().PV->VV[1], 2));
3463 vorticity(p, 1) +=
3464 dx * F1B2
3465 * (m_solver->a_slope(cellId, sysEqn().PV->VV[0], 2) - m_solver->a_slope(cellId, sysEqn().PV->VV[2], 0));
3466 vorticity(p, 2) +=
3467 dx * F1B2
3468 * (m_solver->a_slope(cellId, sysEqn().PV->VV[1], 0) - m_solver->a_slope(cellId, sysEqn().PV->VV[0], 1));
3469 }
3470 if(m_solver->m_vtuVelocityGradientOutput) {
3471 for(MInt i = 0; i < nDim; i++) {
3472 for(MInt j = 0; j < nDim; j++) {
3473 velGrad(p, i * nDim + j) += dx * m_solver->a_slope(cellId, sysEqn().PV->VV[i], j);
3474 }
3475 }
3476 }
3477 if(m_solver->m_vtuLambda2Output) {
3478 MFloat vGrad[3][3]{};
3479 MFloat O[3][3]{};
3480 MFloat S[3][3]{};
3481 MFloat eigenVal[3]{};
3482 for(MInt i = 0; i < nDim; i++) {
3483 for(MInt j = 0; j < nDim; j++) {
3484 O[i][j] = F1B2
3485 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3486 - m_solver->a_slope(cellId, sysEqn().PV->VV[j], i))
3487 / m_solver->m_UInfinity;
3488 S[i][j] = F1B2
3489 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3490 + m_solver->a_slope(cellId, sysEqn().PV->VV[j], i))
3491 / m_solver->m_UInfinity;
3492 }
3493 }
3494 for(MInt i = 0; i < nDim; i++) {
3495 for(MInt j = 0; j < nDim; j++) {
3496 vGrad[i][j] = F0;
3497 for(MInt k = 0; k < nDim; k++) {
3498 vGrad[i][j] += O[i][k] * O[k][j];
3499 vGrad[i][j] += S[i][k] * S[k][j];
3500 }
3501 }
3502 }
3503 maia::math::calcEigenValues(vGrad, eigenVal);
3504 sort(eigenVal, eigenVal + 3);
3505 lambda2(p) += dx * eigenVal[1];
3506 }
3507 pressure(p) += dx * m_solver->a_pvariable(cellId, sysEqn().PV->P);
3508 if(m_solver->m_vtuDensityOutput) density(p) += dx * m_solver->a_pvariable(cellId, sysEqn().PV->RHO);
3509 IF_CONSTEXPR(hasPV_C<SysEqn>::value)
3510 if(m_solver->m_noSpecies > 0) progress(p) += dx * m_solver->a_pvariable(cellId, sysEqn().PV->C);
3511 if(m_solver->m_vtuLevelSetOutput) {
3512 ASSERT(m_solver->m_levelSet, "");
3513 levelSet(p) += dx * m_solver->a_levelSetFunction(cellId, 0);
3514 }
3515 sum += dx;
3516 }
3517 for(MInt i = 0; i < nDim; i++)
3518 velocity(p, i) /= sum;
3519 if(m_solver->m_vtuVorticityOutput)
3520 for(MInt i = 0; i < nDim; i++)
3521 vorticity(p, i) /= sum;
3522 if(m_solver->m_vtuVelocityGradientOutput)
3523 for(MInt i = 0; i < nDim * nDim; i++)
3524 velGrad(p, i) /= sum;
3525 pressure(p) /= sum;
3526 if(m_solver->m_vtuDensityOutput) density(p) /= sum;
3527 if(m_solver->m_noSpecies > 0) progress(p) /= sum;
3528 if(m_solver->m_vtuQCriterionOutput) Q(p) /= sum;
3529 if(m_solver->m_vtuLambda2Output) lambda2(p) /= sum;
3530 }
3531 } else {
3532 for(MInt c = 0; c < (signed)noCells; c++) {
3533 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
3534 for(MInt i = 0; i < nDim; i++)
3535 velocity(c, i) = m_solver->a_pvariable(cellId, sysEqn().PV->VV[i]);
3536 for(MInt i = 0; i < noUserVars; i++)
3537 userVarsOut[c * noUserVars + i] = userVars[noUserVars * cellId + i];
3538 for(MInt i = 0; i < noSolverSpecificVars; i++)
3539 solverSpecificVarsOut[c * noSolverSpecificVars + i] = solverSpecificVars[noSolverSpecificVars * cellId + i];
3540 pressure(c) = m_solver->a_pvariable(cellId, sysEqn().PV->P);
3541 if(m_solver->m_vtuDensityOutput) density(c) = m_solver->a_pvariable(cellId, sysEqn().PV->RHO);
3542 IF_CONSTEXPR(hasPV_C<SysEqn>::value)
3543 if(m_solver->m_noSpecies > 0) progress(c) = m_solver->a_pvariable(cellId, sysEqn().PV->C);
3544 if(m_solver->m_vtuLevelSetOutput) {
3545 ASSERT(!(m_solver->m_levelSetMb && !m_solver->m_constructGField), "");
3546 if(m_solver->m_combustion) {
3547 levelSet(c) = m_solver->a_levelSetFunction(cellId, 0);
3548 } else if(m_solver->m_levelSetMb && m_solver->m_constructGField) {
3549 if(noSolverSpecificVars > 0) {
3550 levelSet(c) = solverSpecificVars[noSolverSpecificVars * cellId];
3551 } else {
3552 ASSERT(false, "levelSet output failed for levelSetMb with constructGField.");
3553 levelSet(c) = F0; // Set 0 as default without message or termination to avoid, if no solverSpecificVars
3554 // not provided.
3555 }
3556 } else {
3557 ASSERT(m_solver->m_levelSet, "");
3558 levelSet(c) = m_solver->a_levelSetFunction(cellId, 0);
3559 }
3560 }
3561 if(m_solver->m_vtuQCriterionOutput) {
3562 Q(c) = F0;
3563 IF_CONSTEXPR(nDim == 3) {
3564 for(MInt i = 0; i < nDim; i++)
3565 for(MInt j = 0; j < nDim; j++)
3566 Q(c) += POW2(F1B2
3567 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3568 - m_solver->a_slope(cellId, sysEqn().PV->VV[j], i)))
3569 - POW2(F1B2
3570 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3571 + m_solver->a_slope(cellId, sysEqn().PV->VV[j], i)));
3572 Q(c) *= F1B2;
3573 }
3574 else {
3575 Q(c) =
3576 F1B2
3577 * (m_solver->a_slope(cellId, sysEqn().PV->VV[1], 0) - m_solver->a_slope(cellId, sysEqn().PV->VV[0], 1));
3578 }
3579 }
3580 if(m_solver->m_vtuLambda2Output) {
3581 MFloat vGrad[3][3];
3582 MFloat O[3][3];
3583 MFloat S[3][3];
3584 MFloat eigenVal[3];
3585 for(MInt i = 0; i < nDim; i++) {
3586 for(MInt j = 0; j < nDim; j++) {
3587 O[i][j] = F1B2
3588 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3589 - m_solver->a_slope(cellId, sysEqn().PV->VV[j], i))
3590 / m_solver->m_UInfinity;
3591 S[i][j] = F1B2
3592 * (m_solver->a_slope(cellId, sysEqn().PV->VV[i], j)
3593 + m_solver->a_slope(cellId, sysEqn().PV->VV[j], i))
3594 / m_solver->m_UInfinity;
3595 }
3596 }
3597 for(MInt i = 0; i < nDim; i++) {
3598 for(MInt j = 0; j < nDim; j++) {
3599 vGrad[i][j] = F0;
3600 for(MInt k = 0; k < nDim; k++) {
3601 vGrad[i][j] += O[i][k] * O[k][j];
3602 vGrad[i][j] += S[i][k] * S[k][j];
3603 }
3604 }
3605 }
3606 maia::math::calcEigenValues(vGrad, eigenVal);
3607 sort(eigenVal, eigenVal + 3);
3608 lambda2(c) = eigenVal[1];
3609 }
3610 if(m_solver->m_vtuVorticityOutput) {
3611 ASSERT(nDim == 3, "");
3612 vorticity(c, 0) =
3613 F1B2
3614 * (m_solver->a_slope(cellId, sysEqn().PV->VV[2], 1) - m_solver->a_slope(cellId, sysEqn().PV->VV[1], 2));
3615 vorticity(c, 1) =
3616 F1B2
3617 * (m_solver->a_slope(cellId, sysEqn().PV->VV[0], 2) - m_solver->a_slope(cellId, sysEqn().PV->VV[2], 0));
3618 vorticity(c, 2) =
3619 F1B2
3620 * (m_solver->a_slope(cellId, sysEqn().PV->VV[1], 0) - m_solver->a_slope(cellId, sysEqn().PV->VV[0], 1));
3621 }
3622 if(m_solver->m_vtuVelocityGradientOutput) {
3623 for(MInt i = 0; i < nDim; i++) {
3624 for(MInt j = 0; j < nDim; j++) {
3625 velGrad(c, i * nDim + j) = m_solver->a_slope(cellId, sysEqn().PV->VV[i], j);
3626 }
3627 }
3628 }
3629 }
3630 }
3631 IF_CONSTEXPR(nDim == 2)
3632 for(MInt c = 0; c < (signed)noCells; c++)
3633 velocity(c, 2) = (float)0.0;
3634
3635 insertDataHeader(reinterpret_cast<char*>(&velocity(0)), velocityMemsize, velocityMemsizeGlobal, velocityOffset);
3636 if(m_solver->m_vtuVorticityOutput)
3637 insertDataHeader(reinterpret_cast<char*>(&vorticity(0)), vorticityMemsize, vorticityMemsizeGlobal,
3638 vorticityOffset);
3639 if(m_solver->m_vtuVelocityGradientOutput)
3640 insertDataHeader(reinterpret_cast<char*>(&velGrad(0)), velGradMemsize, velGradMemsizeGlobal, velGradOffset);
3641 insertDataHeader(reinterpret_cast<char*>(&pressure(0)), pressureMemsize, pressureMemsizeGlobal, pressureOffset);
3642 if(m_solver->m_vtuDensityOutput)
3643 insertDataHeader(reinterpret_cast<char*>(&density(0)), densityMemsize, densityMemsizeGlobal, densityOffset);
3644 if(m_solver->m_noSpecies > 0)
3645 insertDataHeader(reinterpret_cast<char*>(&progress(0)), progressMemsize, progressMemsizeGlobal, progressOffset);
3646 if(m_solver->m_vtuLevelSetOutput)
3647 insertDataHeader(reinterpret_cast<char*>(&levelSet(0)), levelSetMemsize, levelSetMemsizeGlobal, levelSetOffset);
3648 if(m_solver->m_vtuQCriterionOutput)
3649 insertDataHeader(reinterpret_cast<char*>(&Q(0)), QMemsize, QMemsizeGlobal, QOffset);
3650 if(m_solver->m_vtuLambda2Output)
3651 insertDataHeader(reinterpret_cast<char*>(&lambda2(0)), L2Memsize, L2MemsizeGlobal, L2Offset);
3652 if(noUserVars > 0)
3653 insertDataHeader(reinterpret_cast<char*>(&userVarsOut(0)), userVarsMemsize, userVarsMemsizeGlobal,
3654 userVarsOffset);
3655 if(noSolverSpecificVars > 0)
3656 insertDataHeader(reinterpret_cast<char*>(&solverSpecificVarsOut(0)), solverSpecificVarsMemsize,
3657 solverSpecificVarsMemsizeGlobal, solverSpecificVarsOffset);
3658
3659
3660 RECORD_TIMER_STOP(t_prepare);
3661 RECORD_TIMER_START(t_cells);
3662
3663 if(domainId() == 0) {
3664 ofstream ofl;
3665 ofl.open(fileName.c_str(), ios_base::out | ios_base::trunc);
3666 if(ofl.is_open() && ofl.good()) {
3667 //================== VTKFile =================
3668 ofl << "<?xml version=\"1.0\"?>" << endl;
3669 ofl << R"(<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">)"
3670 << endl;
3671 ofl << "<UnstructuredGrid>" << endl;
3672
3673 //================== Dimensions =================
3674 ofl << "<Piece NumberOfPoints=\"" << globalNoPoints << "\" NumberOfCells=\"" << globalNoCells << "\">" << endl;
3675 //================== /Dimensions =================
3676
3677 //================== Points =================
3678 ofl << "<Points>" << endl;
3679 ofl << "<DataArray type=\"" << dataType << R"(" NumberOfComponents="3" format="appended" offset=")" << offset
3680 << "\"/>" << endl;
3681 offset += pointsMemsizeGlobal;
3682 ofl << "</Points>" << endl;
3683 //================== /Points =================
3684
3685 //================== Cells =================
3686 ofl << "<Cells>" << endl;
3687 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="connectivity" format="appended" offset=")" << offset
3688 << "\"/>" << endl;
3689 offset += connectivityMemsizeGlobal;
3690 ofl << "<DataArray type=\"" << uI64DataType << R"(" Name="offsets" format="appended" offset=")" << offset
3691 << "\"/>" << endl;
3692 offset += offsetsMemsizeGlobal;
3693 ofl << "<DataArray type=\"" << uI8DataType << R"(" Name="types" format="appended" offset=")" << offset << "\"/>"
3694 << endl;
3695 offset += typesMemsizeGlobal;
3696 if(m_solver->m_vtuCutCellOutput) {
3697 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="faces" format="appended" offset=")" << offset
3698 << "\"/>" << endl;
3699 offset += facesMemsizeGlobal;
3700 ofl << "<DataArray type=\"" << uI64DataType << R"(" Name="faceoffsets" format="appended" offset=")" << offset
3701 << "\"/>" << endl;
3702 offset += faceoffsetsMemsizeGlobal;
3703 }
3704 ofl << "</Cells>" << endl;
3705 //================== /Cells =================
3706
3707 //========== CellData - PointData ===========
3708 ofl << "<CellData Scalars=\"scalars\">" << endl;
3709 if(m_solver->m_vtuGlobalIdOutput) {
3710 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="globalId" format="appended" offset=")" << offset
3711 << "\"/>" << endl;
3712 offset += globalIdMemsizeGlobal;
3713 }
3714 if(m_solver->m_vtuDomainIdOutput) {
3715 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="domainId" format="appended" offset=")" << offset
3716 << "\"/>" << endl;
3717 offset += domainIdsMemsizeGlobal;
3718 }
3719 if(m_solver->m_vtuWritePointData) ofl << "</CellData>" << endl;
3720 if(m_solver->m_vtuWritePointData) ofl << "<PointData Scalars=\"scalars\">" << endl;
3721 ofl << "<DataArray type=\"" << dataType
3722 << R"(" Name="velocity" NumberOfComponents="3" format="appended" offset=")" << offset << "\"/>" << endl;
3723 offset += velocityMemsizeGlobal;
3724 if(m_solver->m_vtuVorticityOutput)
3725 ofl << "<DataArray type=\"" << dataType
3726 << R"(" Name="vorticity" NumberOfComponents="3" format="appended" offset=")" << offset << "\"/>" << endl;
3727 if(m_solver->m_vtuVorticityOutput) offset += vorticityMemsizeGlobal;
3728 if(m_solver->m_vtuVelocityGradientOutput)
3729 ofl << "<DataArray type=\"" << dataType
3730 << R"(" Name="velocityGradient" NumberOfComponents="9" format="appended" offset=")" << offset << "\"/>"
3731 << endl;
3732 if(m_solver->m_vtuVelocityGradientOutput) offset += velGradMemsizeGlobal;
3733 ofl << "<DataArray type=\"" << dataType << R"(" Name="pressure" format="appended" offset=")" << offset << "\"/>"
3734 << endl;
3735 offset += pressureMemsizeGlobal;
3736 if(m_solver->m_vtuDensityOutput) {
3737 ofl << "<DataArray type=\"" << dataType << R"(" Name="density" format="appended" offset=")" << offset
3738 << "\"/>" << endl;
3739 offset += densityMemsizeGlobal;
3740 }
3741 if(m_solver->m_noSpecies > 0) {
3742 ofl << "<DataArray type=\"" << dataType << R"(" Name="C" format="appended" offset=")" << offset << "\"/>"
3743 << endl;
3744 offset += progressMemsizeGlobal;
3745 }
3746 if(m_solver->m_vtuLevelSetOutput) {
3747 ofl << "<DataArray type=\"" << dataType << R"(" Name="levelSet" format="appended" offset=")" << offset
3748 << "\"/>" << endl;
3749 offset += levelSetMemsizeGlobal;
3750 }
3751 if(m_solver->m_vtuQCriterionOutput) {
3752 IF_CONSTEXPR(nDim == 3)
3753 ofl << "<DataArray type=\"" << dataType << R"(" Name="Q" format="appended" offset=")" << offset << "\"/>"
3754 << endl;
3755 IF_CONSTEXPR(nDim == 2)
3756 ofl << "<DataArray type=\"" << dataType << R"(" Name="vorticity" format="appended" offset=")" << offset
3757 << "\"/>" << endl;
3758 offset += QMemsizeGlobal;
3759 }
3760 if(m_solver->m_vtuLambda2Output) {
3761 ofl << "<DataArray type=\"" << dataType << R"(" Name="lambda2" format="appended" offset=")" << offset
3762 << "\"/>" << endl;
3763 offset += L2MemsizeGlobal;
3764 }
3765 if(m_solver->m_vtuWritePointData)
3766 ofl << "</PointData>" << endl;
3767 else
3768 ofl << "</CellData>" << endl;
3769
3770 if(m_solver->m_levelSetMb && m_solver->m_constructGField ? noSolverSpecificVars > 1
3771 : noSolverSpecificVars > 0) {
3772 ofl << "<DataArray type=\"" << dataType << "\" Name=\"userVars\" NumberOfComponents=\""
3773 << noSolverSpecificVars << "\" format=\"appended\" offset=\"" << offset << "\"/>" << endl;
3774 offset += solverSpecificVarsMemsizeGlobal;
3775 }
3776
3777 if(noUserVars > 0) {
3778 ofl << "<DataArray type=\"" << dataType << "\" Name=\"userVars\" NumberOfComponents=\"" << noUserVars
3779 << "\" format=\"appended\" offset=\"" << offset << "\"/>" << endl;
3780 offset += userVarsMemsizeGlobal;
3781 }
3782 //============ /CellData - /PointData ========
3783
3784 ofl << "</Piece>" << endl;
3785
3786 //================== FieldData =================
3787 ofl << "<FieldData>" << endl;
3788
3789 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="globalTimeStep" format="ascii" NumberOfTuples="1" > )"
3790 << globalTimeStep << " </DataArray>" << endl;
3791 if(m_solver->m_outputPhysicalTime) {
3792 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
3793 << m_solver->m_physicalTime << " </DataArray>" << endl;
3794 ofl << "<DataArray type=\"" << dataType << R"(" Name="internalTime" format="ascii" NumberOfTuples="1" > )"
3795 << m_solver->m_time << " </DataArray>" << endl;
3796 } else {
3797 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
3798 << m_solver->m_time << " </DataArray>" << endl;
3799 ofl << "<DataArray type=\"" << dataType << R"(" Name="physicalTime" format="ascii" NumberOfTuples="1" > )"
3800 << m_solver->m_physicalTime << " </DataArray>" << endl;
3801 }
3802 ofl << "<DataArray type=\"" << dataType << R"(" Name="timeStep" format="ascii" NumberOfTuples="1" > )"
3803 << m_solver->timeStep() << " </DataArray>" << endl;
3804 ofl << "<DataArray type=\"" << dataType << R"(" Name="timeRef" format="ascii" NumberOfTuples="1" > )"
3805 << m_solver->m_timeRef << " </DataArray>" << endl;
3806 ofl << "<DataArray type=\"" << dataType << R"(" Name="Ma" format="ascii" NumberOfTuples="1" > )"
3807 << m_solver->m_Ma << " </DataArray>" << endl;
3808 ofl << "<DataArray type=\"" << dataType << R"(" Name="Re" format="ascii" NumberOfTuples="1" > )"
3809 << m_solver->m_Re << " </DataArray>" << endl;
3810 ofl << "<DataArray type=\"" << dataType << R"(" Name="Re0" format="ascii" NumberOfTuples="1" > )"
3811 << sysEqn().m_Re0 << " </DataArray>" << endl;
3812 ofl << "<DataArray type=\"" << dataType << R"(" Name="Pr" format="ascii" NumberOfTuples="1" > )"
3813 << m_solver->m_Pr << " </DataArray>" << endl;
3814 ofl << "<DataArray type=\"" << dataType << R"(" Name="gamma" format="ascii" NumberOfTuples="1" > )"
3815 << m_solver->sysEqn().gamma_Ref() << " </DataArray>" << endl;
3816 ofl << "<DataArray type=\"" << dataType << R"(" Name="Tref" format="ascii" NumberOfTuples="1" > )"
3817 << m_solver->m_referenceTemperature << " </DataArray>" << endl;
3818 ofl << "<DataArray type=\"" << dataType << R"(" Name="Lref" format="ascii" NumberOfTuples="1" > )"
3819 << m_solver->m_referenceLength << " </DataArray>" << endl;
3820 ofl << "<DataArray type=\"" << dataType << R"(" Name="LrefPhys" format="ascii" NumberOfTuples="1" > )" << F1
3821 << " </DataArray>" << endl;
3822 ofl << "<DataArray type=\"" << dataType << R"(" Name="CFL" format="ascii" NumberOfTuples="1" > )"
3823 << m_solver->m_cfl << " </DataArray>" << endl;
3824 ofl << "<DataArray type=\"" << dataType << R"(" Name="uInf" format="ascii" NumberOfTuples="1" > )"
3825 << m_solver->m_UInfinity << " </DataArray>" << endl;
3826 ofl << "<DataArray type=\"" << dataType << R"(" Name="vInf" format="ascii" NumberOfTuples="1" > )"
3827 << m_solver->m_VInfinity << " </DataArray>" << endl;
3828 IF_CONSTEXPR(nDim == 3)
3829 ofl << "<DataArray type=\"" << dataType << R"(" Name="wInf" format="ascii" NumberOfTuples="1" > )"
3830 << m_solver->m_WInfinity << " </DataArray>" << endl;
3831 ofl << "<DataArray type=\"" << dataType << R"(" Name="pInf" format="ascii" NumberOfTuples="1" > )"
3832 << m_solver->m_PInfinity << " </DataArray>" << endl;
3833 ofl << "<DataArray type=\"" << dataType << R"(" Name="rhoInf" format="ascii" NumberOfTuples="1" > )"
3834 << m_solver->m_rhoInfinity << " </DataArray>" << endl;
3835 ofl << "<DataArray type=\"" << dataType << R"(" Name="TInf" format="ascii" NumberOfTuples="1" > )"
3836 << m_solver->m_TInfinity << " </DataArray>" << endl;
3837 ofl << "<DataArray type=\"" << dataType << R"(" Name="muInf" format="ascii" NumberOfTuples="1" > )"
3838 << sysEqn().m_muInfinity << " </DataArray>" << endl;
3839 ofl << "<DataArray type=\"" << dataType << R"(" Name="length0" format="ascii" NumberOfTuples="1" > )"
3840 << m_solver->c_cellLengthAtLevel(0) << " </DataArray>" << endl;
3841 ofl << "<DataArray type=\"" << dataType << R"(" Name="lengthMax" format="ascii" NumberOfTuples="1" > )"
3842 << m_solver->c_cellLengthAtLevel(m_solver->minLevel()) << " </DataArray>" << endl;
3843 ofl << "<DataArray type=\"" << dataType << R"(" Name="lengthMin" format="ascii" NumberOfTuples="1" > )"
3844 << m_solver->c_cellLengthAtLevel(m_solver->maxLevel()) << " </DataArray>" << endl;
3845 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="minLevel" format="ascii" NumberOfTuples="1" > )"
3846 << m_solver->minLevel() << " </DataArray>" << endl;
3847 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="maxLevel" format="ascii" NumberOfTuples="1" > )"
3848 << m_solver->maxLevel() << " </DataArray>" << endl;
3849 if(noDomains() > 1) {
3850 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="domainOffsets" format="ascii" NumberOfTuples=")"
3851 << noDomains() << "\" > ";
3852 for(MInt i = 0; i < noDomains(); i++) {
3853 ofl << m_solver->domainOffset(i) << " ";
3854 }
3855 ofl << " </DataArray>" << endl;
3856 }
3857 ofl << "</FieldData>" << endl;
3858 //================== /FieldData =================
3859
3860 ofl << "</UnstructuredGrid>" << endl;
3861
3862 //================== AppendedData =================
3863 ofl << "<AppendedData encoding=\"raw\">" << endl;
3864 ofl << "_";
3865 ofl.close();
3866 ofl.clear();
3867
3868 if(m_solver->m_vtuSaveHeaderTesting) {
3869 cout << system(("cp " + fileName + " " + fileName + "_header_testing").c_str());
3870 }
3871 } else {
3872 cerr << "ERROR! COULD NOT OPEN FILE " << fileName << " for writing! (1)" << endl;
3873 }
3874 }
3875
3876 //###################################################
3877 //################### parallel IO ###################
3878 MPI_File file = nullptr;
3879 MInt rcode = MPI_File_open(mpiComm(), const_cast<char*>(fileName.c_str()), MPI_MODE_APPEND | MPI_MODE_WRONLY,
3880 MPI_INFO_NULL, &file, AT_);
3881 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error opening file " + fileName);
3882
3883 //=================== Points =======================
3884 rcode = writeVtuArrayParallel(file, &points(0), pointsOffset, pointsMemsize, pointsMemsizeGlobal);
3885 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (1) writing to file " + fileName);
3886
3887 //=================== Connectivity =======================
3888 rcode = writeVtuArrayParallel(file, &connectivity(0), connectivityOffset, connectivityMemsize,
3889 connectivityMemsizeGlobal);
3890 if(rcode != MPI_SUCCESS) {
3891 mTerm(1, AT_, "Error (2) writing to file " + fileName);
3892 }
3893
3894 //=================== Offsets =======================
3895 rcode = writeVtuArrayParallel(file, &offsets(0), offsetsOffset, offsetsMemsize, offsetsMemsizeGlobal);
3896 if(rcode != MPI_SUCCESS) {
3897 mTerm(1, AT_, "Error (3) writing to file " + fileName);
3898 }
3899
3900 //=================== Types =======================
3901 rcode = writeVtuArrayParallel(file, &types(0), typesOffset, typesMemsize, typesMemsizeGlobal);
3902 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (4) writing to file " + fileName);
3903
3904 //=================== Faces =======================
3905 if(m_solver->m_vtuCutCellOutput)
3906 rcode = writeVtuArrayParallel(file, &faces(0), facesOffset, facesMemsize, facesMemsizeGlobal);
3907 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (5) writing to file " + fileName);
3908
3909 //=================== FaceOffsets =======================
3910 if(m_solver->m_vtuCutCellOutput)
3911 rcode =
3912 writeVtuArrayParallel(file, &faceoffsets(0), faceoffsetsOffset, faceoffsetsMemsize, faceoffsetsMemsizeGlobal);
3913 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (5) writing to file " + fileName);
3914
3915 RECORD_TIMER_STOP(t_cells);
3916 RECORD_TIMER_START(t_variables);
3917
3918 //=================== GlobalId =======================
3919 if(m_solver->m_vtuGlobalIdOutput)
3920 rcode = writeVtuArrayParallel(file, &globalId(0), globalIdOffset, globalIdMemsize, globalIdMemsizeGlobal);
3921 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (5) writing to file " + fileName);
3922 //=================== GlobalId =======================
3923
3924 if(m_solver->m_vtuDomainIdOutput)
3925 rcode = writeVtuArrayParallel(file, &domainIds(0), domainIdsOffset, domainIdsMemsize, domainIdsMemsizeGlobal);
3926 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (5) writing to file " + fileName);
3927
3928 //=================== velocity =======================
3929 rcode = writeVtuArrayParallel(file, &velocity(0), velocityOffset, velocityMemsize, velocityMemsizeGlobal);
3930 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (6) writing to file " + fileName);
3931
3932 //=================== vorticity =======================
3933 if(m_solver->m_vtuVorticityOutput)
3934 rcode = writeVtuArrayParallel(file, &vorticity(0), vorticityOffset, vorticityMemsize, vorticityMemsizeGlobal);
3935 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (7) writing to file " + fileName);
3936
3937 //=================== velGrad =======================
3938 if(m_solver->m_vtuVelocityGradientOutput)
3939 rcode = writeVtuArrayParallel(file, &velGrad(0), velGradOffset, velGradMemsize, velGradMemsizeGlobal);
3940 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (7) writing to file " + fileName);
3941
3942 //=================== pressure =======================
3943 rcode = writeVtuArrayParallel(file, &pressure(0), pressureOffset, pressureMemsize, pressureMemsizeGlobal);
3944 if(rcode != MPI_SUCCESS) {
3945 mTerm(1, AT_, "Error (8) writing to file " + fileName);
3946 }
3947
3948 //=================== density =======================
3949 if(m_solver->m_vtuDensityOutput)
3950 rcode = writeVtuArrayParallel(file, &density(0), densityOffset, densityMemsize, densityMemsizeGlobal);
3951 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (9) writing to file " + fileName);
3952
3953 if(m_solver->m_noSpecies > 0) {
3954 //=================== progress =======================
3955 rcode = writeVtuArrayParallel(file, &progress(0), progressOffset, progressMemsize, progressMemsizeGlobal);
3956 if(rcode != MPI_SUCCESS) {
3957 mTerm(1, AT_, "Error (9) writing to file " + fileName);
3958 }
3959 }
3960 //=================== levelSet =======================
3961 if(m_solver->m_vtuLevelSetOutput) {
3962 rcode = writeVtuArrayParallel(file, &levelSet(0), levelSetOffset, levelSetMemsize, levelSetMemsizeGlobal);
3963 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (9) writing to file " + fileName);
3964 }
3965 //=================== Q =======================
3966 if(m_solver->m_vtuQCriterionOutput) rcode = writeVtuArrayParallel(file, &Q(0), QOffset, QMemsize, QMemsizeGlobal);
3967 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (10) writing to file " + fileName);
3968
3969 //=================== lambda2 =======================
3970 if(m_solver->m_vtuLambda2Output)
3971 rcode = writeVtuArrayParallel(file, &lambda2(0), L2Offset, L2Memsize, L2MemsizeGlobal);
3972 if(rcode != MPI_SUCCESS) {
3973 mTerm(1, AT_, "Error (11) writing to file " + fileName);
3974 }
3975
3976 //=================== solverSpecificVars =======================
3977 if(m_solver->m_levelSetMb && m_solver->m_constructGField ? noSolverSpecificVars > 1 : noSolverSpecificVars > 0)
3978 rcode = writeVtuArrayParallel(file, &solverSpecificVarsOut(0), solverSpecificVarsOffset,
3979 solverSpecificVarsMemsize, solverSpecificVarsMemsizeGlobal);
3980 if(rcode != MPI_SUCCESS) {
3981 mTerm(1, AT_, "Error (12) writing to file " + fileName);
3982 }
3983
3984 //=================== userVars =======================
3985 if(noUserVars > 0)
3986 rcode = writeVtuArrayParallel(file, &userVarsOut(0), userVarsOffset, userVarsMemsize, userVarsMemsizeGlobal);
3987 if(rcode != MPI_SUCCESS) {
3988 mTerm(1, AT_, "Error (13) writing to file " + fileName);
3989 }
3990
3991 MPI_File_close(&file, AT_);
3992 //###################################################
3993 //###################################################
3994
3995 if(domainId() == 0) {
3996 ofstream ofl;
3997 ofl.open(fileName.c_str(), ios_base::out | ios_base::app);
3998 if(ofl.is_open() && ofl.good()) {
3999 ofl << endl;
4000 ofl << "</AppendedData>" << endl;
4001 ofl << "</VTKFile>" << endl;
4002 ofl.close();
4003 ofl.clear();
4004 //================== /AppendedData =================
4005 //================== /VTKFile ======================
4006 } else {
4007 cerr << "ERROR! COULD NOT OPEN FILE " << fileName << " for writing! (2)" << endl;
4008 }
4009 }
4010
4011 //-----------------------
4012 //------ QOUT.pvd -------
4013 //-----------------------
4014
4015 if(domainId() == 0) {
4016 if(m_solver->m_firstUseWriteVtuOutputParallelQout) {
4017 m_solver->m_firstUseWriteVtuOutputParallelQout = !m_solver->m_restart;
4018 }
4019 ofstream ofile;
4020 ofstream ofile2;
4021 if(m_solver->m_firstUseWriteVtuOutputParallelQout)
4022 ofile.open("out/QOUT.pvd.tmp", ios_base::out | ios_base::trunc);
4023 else
4024 ofile.open("out/QOUT.pvd.tmp", ios_base::out | ios_base::app);
4025 if(ofile.is_open() && ofile.good()) {
4026 if(m_solver->m_firstUseWriteVtuOutputParallelQout) {
4027 ofile << R"(<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">)" << endl;
4028 ofile << "<Collection>" << endl;
4029 }
4030 const size_t pos = fileName.find(m_solver->m_solutionOutput);
4031 MString fname =
4032 (pos == string::npos) ? fileName : fileName.substr(pos + m_solver->m_solutionOutput.size(), string::npos);
4033 ofile << "<DataSet part=\"" << 0 << "\" timestep=\"" << globalTimeStep << "\" file=\""
4034 << "./" << fname << "\"/>" << endl;
4035 ofile.close();
4036 ofile.clear();
4037 } else {
4038 cerr << "Error opening file out/QOUT.pvd.tmp" << endl;
4039 }
4040 m_log << "Command executed with return value: " << system("cp out/QOUT.pvd.tmp out/QOUT.pvd") << " at " << AT_
4041 << endl;
4042
4043 ofile2.open("out/QOUT.pvd", ios_base::out | ios_base::app);
4044 if(ofile2.is_open() && ofile2.good()) {
4045 ofile2 << "</Collection>" << endl;
4046 ofile2 << "</VTKFile>" << endl;
4047 ofile2.close();
4048 ofile2.clear();
4049 } else {
4050 cerr << "Error opening file out/QOUT.pvd" << endl;
4051 }
4052 m_solver->m_firstUseWriteVtuOutputParallelQout = false;
4053 }
4054
4055 RECORD_TIMER_STOP(t_variables);
4056 RECORD_TIMER_START(t_geometry);
4057 }
4058
4059
4060 if(m_solver->m_vtuCutCellOutput
4061 && (m_solver->m_vtuGeometryOutput.size() > 0 || string2enum(m_solver->m_outputFormat) == VTP)) {
4062 if(m_solver->m_vtuLevelThreshold < m_solver->maxRefinementLevel()) {
4063 if(m_solver->m_extractedCells) {
4064 delete m_solver->m_extractedCells;
4065 m_solver->m_extractedCells = nullptr;
4066 }
4067 if(m_solver->m_gridPoints) {
4068 delete m_solver->m_gridPoints;
4069 m_solver->m_gridPoints = nullptr;
4070 }
4071 m_solver->extractPointIdsFromGrid(m_solver->m_extractedCells, m_solver->m_gridPoints, false,
4072 m_solver->m_splitChildToSplitCell, m_solver->maxRefinementLevel(),
4073 m_solver->m_vtuCoordinatesThreshold);
4074 facesSize = 0;
4075 noPolyhedra = 0;
4076 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
4077 MInt cellId = m_solver->m_extractedCells->a[c].m_cellId;
4078 MInt bndryId = m_solver->a_bndryId(cellId);
4079 if(bndryId > -1) {
4080 polyhedralCell[c] = noPolyhedra;
4081 noPolyhedra++;
4082 }
4083 }
4084 connSize = vtuAssembleFaceStream(facestream, conn, facesSize, polyhedralCell, noPolyhedra, false);
4085 }
4086
4087 if(!(facestream == nullptr || conn == nullptr)) {
4088 uint64_t noPolygons = 0;
4089 uint64_t geomConnSize = 0;
4090 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
4091 const MInt bndryId = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
4092 const MInt pc = polyhedralCell[c];
4093 if(bndryId < 0) continue;
4094 uint_t cnt = 1;
4095 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
4096 if(m_solver->m_vtuGeometryOutput.count(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_bndryCndId) == 0)
4097 continue;
4098 noPolygons++;
4099 geomConnSize += (uint64_t)facestream[pc][cnt];
4100 cnt += (uint_t)(facestream[pc][cnt] + 1);
4101 }
4102 }
4103
4104 ScratchSpace<float> geomPoints(geomConnSize + 4, 3, AT_, "geomPoints");
4105 ScratchSpace<uint_t> geomConn(geomConnSize + 4, AT_, "geomConn");
4106 ScratchSpace<uint_t> geomOffsets(noPolygons + 4, AT_, "geomOffsets");
4107 ScratchSpace<uint_t> cutPoints(m_solver->m_gridPoints->size(), AT_, "cutPoints");
4108 cutPoints.fill(std::numeric_limits<uint_t>::max());
4109 uint64_t noGeomPoints = 0;
4110 geomConnSize = 0;
4111 noPolygons = 0;
4112 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
4113 const MInt bndryId = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
4114 const MInt pc = polyhedralCell[c];
4115 if(pc < 0) continue;
4116 if(bndryId < 0) continue;
4117 if(facestream[pc].empty()) continue;
4118 ASSERT((signed)facestream[pc][0] >= m_solver->m_bndryCells->a[bndryId].m_noSrfcs, "");
4119 uint_t cnt = 1;
4120 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
4121 if(m_solver->m_vtuGeometryOutput.count(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_bndryCndId) == 0)
4122 continue;
4123 ASSERT(cnt < facestream[pc].size(), "");
4124 MInt noCutPoints = (signed)facestream[pc][cnt];
4125 cnt++;
4126 for(MInt cp = 0; cp < noCutPoints; cp++) {
4127 ASSERT(cnt < facestream[pc].size(), "");
4128 uint_t gridPointId = (uint_t)facestream[pc][cnt];
4129 uint_t pId = cutPoints[gridPointId];
4130 if(pId >= (unsigned)m_solver->m_gridPoints->size()) {
4131 pId = noGeomPoints;
4132 cutPoints[gridPointId] = pId;
4133 for(MInt i = 0; i < nDim; i++) {
4134 geomPoints((signed)pId, i) = m_solver->m_gridPoints->a[gridPointId].m_coordinates[i];
4135 }
4136 noGeomPoints++;
4137 }
4138 geomConn[geomConnSize] = pId;
4139 cnt++;
4140 geomConnSize++;
4141 }
4142 geomOffsets[noPolygons] = geomConnSize;
4143 noPolygons++;
4144 }
4145 }
4146
4147 ScratchSpace<uint64_t> noGeomPointsPerDomain(noDomains(), AT_, "noGeomPointsPerDomain");
4148 ScratchSpace<uint64_t> noPolygonsPerDomain(noDomains(), AT_, "noPolygonsPerDomain");
4149 ScratchSpace<uint64_t> geomConnSizePerDomain(noDomains(), AT_, "geomConnSizePerDomain");
4150 ScratchSpace<uint64_t> geomDataPerDomain(noDomains(), 3, AT_, "geomDataPerDomain");
4151 geomDataPerDomain(domainId(), 0) = noGeomPoints;
4152 geomDataPerDomain(domainId(), 1) = noPolygons;
4153 geomDataPerDomain(domainId(), 2) = geomConnSize;
4154 MPI_Allgather(MPI_IN_PLACE, 3, MPI_UINT64_T, &geomDataPerDomain[0], 3, MPI_UINT64_T, mpiComm(), AT_,
4155 "MPI_IN_PLACE", "geomDataPerDomain[0]");
4156 for(MInt d = 0; d < noDomains(); d++) {
4157 noGeomPointsPerDomain[d] = geomDataPerDomain(d, 0);
4158 noPolygonsPerDomain[d] = geomDataPerDomain(d, 1);
4159 geomConnSizePerDomain[d] = geomDataPerDomain(d, 2);
4160 }
4161 uint64_t globalNoGeomPoints = 0;
4162 uint64_t globalGeomPointOffset = 0;
4163 uint64_t globalNoPolygons = 0;
4164 uint64_t globalPolygonsOffset = 0;
4165 uint64_t globalGeomConnSize = 0;
4166 uint64_t globalGeomConnSizeOffset = 0;
4167 for(MInt d = 0; d < noDomains(); d++) {
4168 globalNoGeomPoints += noGeomPointsPerDomain[d];
4169 globalNoPolygons += noPolygonsPerDomain[d];
4170 globalGeomConnSize += geomConnSizePerDomain[d];
4171 }
4172 for(MInt d = 0; d < domainId(); d++) {
4173 globalGeomPointOffset += noGeomPointsPerDomain[d];
4174 globalPolygonsOffset += noPolygonsPerDomain[d];
4175 globalGeomConnSizeOffset += geomConnSizePerDomain[d];
4176 }
4177
4178 for(MUint p = 0; p < noPolygons; p++) {
4179 geomOffsets[p] += globalGeomConnSizeOffset;
4180 }
4181 for(MUint c = 0; c < geomConnSize; c++) {
4182 geomConn[c] += globalGeomPointOffset;
4183 }
4184
4185 if(globalNoPolygons == 0) {
4186 RECORD_TIMER_STOP(t_geometry);
4187 RECORD_TIMER_STOP(t_vtutotal);
4188 DISPLAY_TIMER_INTERM(t_vtutotal);
4189 return;
4190 }
4191
4192
4193 //-----------
4194 // points
4195 uint64_t geomPointsMemsize = 3 * noGeomPoints * sizeof(float);
4196 uint64_t geomPointsMemsizeGlobal = 3 * globalNoGeomPoints * sizeof(float);
4197 uint64_t geomPointsOffset = 3 * globalGeomPointOffset * sizeof(float);
4198 insertDataHeader(reinterpret_cast<char*>(&geomPoints(0)), geomPointsMemsize, geomPointsMemsizeGlobal,
4199 geomPointsOffset);
4200
4201 //-----------
4202 // connectivity
4203 uint64_t geomConnMemsize = geomConnSize * sizeof(uint_t);
4204 uint64_t geomConnOffset = globalGeomConnSizeOffset * sizeof(uint_t);
4205 uint64_t geomConnMemsizeGlobal = globalGeomConnSize * sizeof(uint_t);
4206 insertDataHeader(reinterpret_cast<char*>(&geomConn(0)), geomConnMemsize, geomConnMemsizeGlobal, geomConnOffset);
4207
4208 //-----------
4209 // offsets
4210 uint64_t geomOffsetsMemsize = noPolygons * sizeof(uint_t);
4211 uint64_t geomOffsetsOffset = globalPolygonsOffset * sizeof(uint_t);
4212 uint64_t geomOffsetsMemsizeGlobal = globalNoPolygons * sizeof(uint_t);
4213 insertDataHeader(reinterpret_cast<char*>(&geomOffsets(0)), geomOffsetsMemsize, geomOffsetsMemsizeGlobal,
4214 geomOffsetsOffset);
4215
4216 //-----------
4217 // variables
4218 uint64_t bodyIdMemsize = noPolygons * sizeof(uint_t);
4219 uint64_t velocityMemsize = 3 * noPolygons * sizeof(float);
4220 uint64_t domainIdsMemsize = noPolygons * sizeof(uint_t);
4221 uint64_t pressureMemsize = noPolygons * sizeof(float);
4222 uint64_t densityMemsize = noPolygons * sizeof(float);
4223 uint64_t normalsMemsize = 3 * noPolygons * sizeof(float);
4224 uint64_t shearMemsize = 3 * noPolygons * sizeof(float);
4225 uint64_t bodyIdOffset = globalPolygonsOffset * sizeof(uint_t);
4226 uint64_t velocityOffset = 3 * globalPolygonsOffset * sizeof(float);
4227 uint64_t domainIdsOffset = globalPolygonsOffset * sizeof(uint_t);
4228 uint64_t pressureOffset = globalPolygonsOffset * sizeof(float);
4229 uint64_t densityOffset = globalPolygonsOffset * sizeof(float);
4230 uint64_t normalsOffset = 3 * globalPolygonsOffset * sizeof(float);
4231 uint64_t shearOffset = 3 * globalPolygonsOffset * sizeof(float);
4232 uint64_t bodyIdMemsizeGlobal = globalNoPolygons * sizeof(uint_t);
4233 uint64_t velocityMemsizeGlobal = 3 * globalNoPolygons * sizeof(float);
4234 uint64_t domainIdsMemsizeGlobal = globalNoPolygons * sizeof(uint_t);
4235 uint64_t pressureMemsizeGlobal = globalNoPolygons * sizeof(float);
4236 uint64_t densityMemsizeGlobal = globalNoPolygons * sizeof(float);
4237 uint64_t normalsMemsizeGlobal = 3 * globalNoPolygons * sizeof(float);
4238 uint64_t shearMemsizeGlobal = 3 * globalNoPolygons * sizeof(float);
4239 const MInt bodyIdMaxSize = estimateMemorySizeSolverwise(noPolygons, noPolygonsPerDomain, sizeof(uint_t));
4240 const MInt vectorMaxSize = estimateMemorySizeSolverwise(noPolygons, noPolygonsPerDomain, 3 * sizeof(float));
4241 const MInt scalarMaxSize = estimateMemorySizeSolverwise(noPolygons, noPolygonsPerDomain, sizeof(float));
4242 ScratchSpace<uint_t> bodyId(bodyIdMaxSize, AT_, "bodyId");
4243 ScratchSpace<float> velocity(vectorMaxSize, 3, AT_, "velocity");
4244 ScratchSpace<uint_t> domainIds(m_solver->m_vtuGeometryOutputExtended ? bodyIdMaxSize : 1, AT_, "domainIds");
4245 ScratchSpace<float> pressure(m_solver->m_vtuGeometryOutputExtended ? scalarMaxSize : 1, AT_, "pressure");
4246 ScratchSpace<float> density(m_solver->m_vtuGeometryOutputExtended ? scalarMaxSize : 1, AT_, "density");
4247 ScratchSpace<float> normals(m_solver->m_vtuGeometryOutputExtended ? vectorMaxSize : 1, 3, AT_, "normals");
4248 ScratchSpace<float> shear(m_solver->m_vtuGeometryOutputExtended ? vectorMaxSize : 1, 3, AT_, "shear");
4249 const MFloat F1BRe = F1 / (m_solver->m_referenceLength * sysEqn().m_Re0);
4250 MInt cnt = 0;
4251 for(MInt c = 0; c < m_solver->m_extractedCells->size(); c++) {
4252 MInt bndryId = m_solver->a_bndryId(m_solver->m_extractedCells->a[c].m_cellId);
4253 MInt pc = polyhedralCell[c];
4254 if(pc < 0) continue;
4255 if(bndryId < 0) continue;
4256 for(MInt srfc = 0; srfc < m_solver->m_bndryCells->a[bndryId].m_noSrfcs; srfc++) {
4257 if(m_solver->m_vtuGeometryOutput.count(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_bndryCndId) == 0)
4258 continue;
4259 bodyId(cnt) = (uint_t)m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_bodyId[0];
4260 if(m_solver->m_internalBodyId != nullptr && (signed)bodyId(cnt) > -1
4261 && (signed)bodyId(cnt) >= m_solver->m_noEmbeddedBodies)
4262 bodyId(cnt) = (uint_t)(m_solver->m_internalBodyId[bodyId(cnt)]);
4263 if(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_bodyId[0] < 0)
4264 bodyId(cnt) = std::numeric_limits<uint_t>::max();
4265 for(MInt i = 0; i < nDim; i++) {
4266 velocity(cnt, i) =
4267 (float)m_solver->m_bndryCells->a[bndryId].m_srfcVariables[srfc]->m_primVars[sysEqn().PV->VV[i]];
4268 }
4269
4270 if(m_solver->m_vtuGeometryOutputExtended) {
4271 domainIds(cnt) = (uint_t)domainId();
4272 pressure(cnt) = (float)m_solver->m_bndryCells->a[bndryId].m_srfcVariables[srfc]->m_primVars[sysEqn().PV->P];
4273 density(cnt) =
4274 (float)m_solver->m_bndryCells->a[bndryId].m_srfcVariables[srfc]->m_primVars[sysEqn().PV->RHO];
4275 MFloat mue = m_solver->sysEqn().sutherlandLaw(m_solver->sysEqn().temperature_ES(
4276 m_solver->m_bndryCells->a[bndryId].m_srfcVariables[srfc]->m_primVars[sysEqn().PV->RHO],
4277 m_solver->m_bndryCells->a[bndryId].m_srfcVariables[srfc]->m_primVars[sysEqn().PV->P]));
4278 for(MInt i = 0; i < nDim; i++) {
4279 normals(cnt, i) = (float)(m_solver->m_bndryCells->a[bndryId].m_srfcs[srfc]->m_normalVector[i]);
4280 shear(cnt, i) = (float)(F1BRe * mue
4281 * m_solver->m_bndryCells->a[bndryId]
4282 .m_srfcVariables[srfc]
4283 ->m_normalDeriv[sysEqn().PV->VV[i]]);
4284 }
4285 }
4286 cnt++;
4287 }
4288 }
4289 insertDataHeader(reinterpret_cast<char*>(&bodyId(0)), bodyIdMemsize, bodyIdMemsizeGlobal, bodyIdOffset);
4290 insertDataHeader(reinterpret_cast<char*>(&velocity(0)), velocityMemsize, velocityMemsizeGlobal, velocityOffset);
4291 if(m_solver->m_vtuGeometryOutputExtended) {
4292 insertDataHeader(reinterpret_cast<char*>(&domainIds(0)), domainIdsMemsize, domainIdsMemsizeGlobal,
4293 domainIdsOffset);
4294 insertDataHeader(reinterpret_cast<char*>(&pressure(0)), pressureMemsize, pressureMemsizeGlobal, pressureOffset);
4295 insertDataHeader(reinterpret_cast<char*>(&density(0)), densityMemsize, densityMemsizeGlobal, densityOffset);
4296 insertDataHeader(reinterpret_cast<char*>(&normals(0)), normalsMemsize, normalsMemsizeGlobal, normalsOffset);
4297 insertDataHeader(reinterpret_cast<char*>(&shear(0)), shearMemsize, shearMemsizeGlobal, shearOffset);
4298 }
4299
4300 if(m_solver->m_vtuWriteGeometryFile) {
4301 if(domainId() == 0) {
4302 offset = 0;
4303 ofstream ofl;
4304 ofl.open(geomFileName.c_str(), ios_base::out | ios_base::trunc);
4305 if(ofl.is_open() && ofl.good()) {
4306 //================== VTKFile =================
4307 ofl << "<?xml version=\"1.0\"?>" << endl;
4308 ofl << R"(<VTKFile type="PolyData" version="1.0" byte_order="LittleEndian" header_type="UInt64">)" << endl;
4309 ofl << "<PolyData>" << endl;
4310
4311 //================== FieldData =================
4312 ofl << "<FieldData>" << endl;
4313
4314 ofl << "<DataArray type=\"" << uIDataType
4315 << R"(" Name="globalTimeStep" format="ascii" NumberOfTuples="1" > )" << globalTimeStep
4316 << " </DataArray>" << endl;
4317 if(m_solver->m_outputPhysicalTime) {
4318 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
4319 << m_solver->m_physicalTime << " </DataArray>" << endl;
4320 ofl << "<DataArray type=\"" << dataType << R"(" Name="internalTime" format="ascii" NumberOfTuples="1" > )"
4321 << m_solver->m_time << " </DataArray>" << endl;
4322 } else {
4323 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
4324 << m_solver->m_time << " </DataArray>" << endl;
4325 ofl << "<DataArray type=\"" << dataType << R"(" Name="physicalTime" format="ascii" NumberOfTuples="1" > )"
4326 << m_solver->m_physicalTime << " </DataArray>" << endl;
4327 }
4328 ofl << "<DataArray type=\"" << dataType << R"(" Name="timeStep" format="ascii" NumberOfTuples="1" > )"
4329 << m_solver->timeStep() << " </DataArray>" << endl;
4330 ofl << "<DataArray type=\"" << dataType << R"(" Name="timeRef" format="ascii" NumberOfTuples="1" > )"
4331 << m_solver->m_timeRef << " </DataArray>" << endl;
4332 ofl << "<DataArray type=\"" << dataType << R"(" Name="Ma" format="ascii" NumberOfTuples="1" > )"
4333 << m_solver->m_Ma << " </DataArray>" << endl;
4334 ofl << "<DataArray type=\"" << dataType << R"(" Name="Re" format="ascii" NumberOfTuples="1" > )"
4335 << m_solver->m_Re << " </DataArray>" << endl;
4336 ofl << "<DataArray type=\"" << dataType << R"(" Name="Re0" format="ascii" NumberOfTuples="1" > )"
4337 << sysEqn().m_Re0 << " </DataArray>" << endl;
4338 ofl << "<DataArray type=\"" << dataType << R"(" Name="Pr" format="ascii" NumberOfTuples="1" > )"
4339 << m_solver->m_Pr << " </DataArray>" << endl;
4340 ofl << "<DataArray type=\"" << dataType << R"(" Name="gamma" format="ascii" NumberOfTuples="1" > )"
4341 << m_solver->sysEqn().gamma_Ref() << " </DataArray>" << endl;
4342 ofl << "<DataArray type=\"" << dataType << R"(" Name="Tref" format="ascii" NumberOfTuples="1" > )"
4343 << m_solver->m_referenceTemperature << " </DataArray>" << endl;
4344 ofl << "<DataArray type=\"" << dataType << R"(" Name="Lref" format="ascii" NumberOfTuples="1" > )"
4345 << m_solver->m_referenceLength << " </DataArray>" << endl;
4346 ofl << "<DataArray type=\"" << dataType << R"(" Name="LrefPhys" format="ascii" NumberOfTuples="1" > )" << F1
4347 << " </DataArray>" << endl;
4348 ofl << "<DataArray type=\"" << dataType << R"(" Name="CFL" format="ascii" NumberOfTuples="1" > )"
4349 << m_solver->m_cfl << " </DataArray>" << endl;
4350 ofl << "<DataArray type=\"" << dataType << R"(" Name="uInf" format="ascii" NumberOfTuples="1" > )"
4351 << m_solver->m_UInfinity << " </DataArray>" << endl;
4352 ofl << "<DataArray type=\"" << dataType << R"(" Name="vInf" format="ascii" NumberOfTuples="1" > )"
4353 << m_solver->m_VInfinity << " </DataArray>" << endl;
4354 IF_CONSTEXPR(nDim == 3)
4355 ofl << "<DataArray type=\"" << dataType << R"(" Name="wInf" format="ascii" NumberOfTuples="1" > )"
4356 << m_solver->m_WInfinity << " </DataArray>" << endl;
4357 ofl << "<DataArray type=\"" << dataType << R"(" Name="pInf" format="ascii" NumberOfTuples="1" > )"
4358 << m_solver->m_PInfinity << " </DataArray>" << endl;
4359 ofl << "<DataArray type=\"" << dataType << R"(" Name="rhoInf" format="ascii" NumberOfTuples="1" > )"
4360 << m_solver->m_rhoInfinity << " </DataArray>" << endl;
4361 ofl << "<DataArray type=\"" << dataType << R"(" Name="TInf" format="ascii" NumberOfTuples="1" > )"
4362 << m_solver->m_TInfinity << " </DataArray>" << endl;
4363 ofl << "<DataArray type=\"" << dataType << R"(" Name="muInf" format="ascii" NumberOfTuples="1" > )"
4364 << sysEqn().m_muInfinity << " </DataArray>" << endl;
4365 ofl << "<DataArray type=\"" << dataType << R"(" Name="length0" format="ascii" NumberOfTuples="1" > )"
4366 << m_solver->c_cellLengthAtLevel(0) << " </DataArray>" << endl;
4367 ofl << "<DataArray type=\"" << dataType << R"(" Name="lengthMax" format="ascii" NumberOfTuples="1" > )"
4368 << m_solver->c_cellLengthAtLevel(m_solver->minLevel()) << " </DataArray>" << endl;
4369 ofl << "<DataArray type=\"" << dataType << R"(" Name="lengthMin" format="ascii" NumberOfTuples="1" > )"
4370 << m_solver->c_cellLengthAtLevel(m_solver->maxLevel()) << " </DataArray>" << endl;
4371 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="minLevel" format="ascii" NumberOfTuples="1" > )"
4372 << m_solver->minLevel() << " </DataArray>" << endl;
4373 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="maxLevel" format="ascii" NumberOfTuples="1" > )"
4374 << m_solver->maxLevel() << " </DataArray>" << endl;
4375 if(noDomains() > 1) {
4376 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="domainOffsets" format="ascii" NumberOfTuples=")"
4377 << noDomains() << "\" > ";
4378 for(MInt i = 0; i < noDomains(); i++) {
4379 ofl << m_solver->domainOffset(i) << " ";
4380 }
4381 ofl << " </DataArray>" << endl;
4382 }
4383
4384 ofl << "</FieldData>" << endl;
4385 //================== /FieldData =================
4386
4387
4388 //================== Dimensions =================
4389 ofl << "<Piece NumberOfPoints=\"" << globalNoGeomPoints << "\" NumberOfPolys=\"" << globalNoPolygons
4390 << "\">" << endl;
4391 //================== /Dimensions ================
4392
4393
4394 //================== Points =================
4395 ofl << "<Points>" << endl;
4396 ofl << "<DataArray type=\"" << dataType << R"(" NumberOfComponents="3" format="appended" offset=")"
4397 << offset << "\"/>" << endl;
4398 offset += geomPointsMemsizeGlobal;
4399 ofl << "</Points>" << endl;
4400 //================== /Points =================
4401
4402
4403 //================== Polys =================
4404 ofl << "<Polys>" << endl;
4405
4406 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="connectivity" format="appended" offset=")" << offset
4407 << "\"/>" << endl;
4408 offset += geomConnMemsizeGlobal;
4409
4410 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="offsets" format="appended" offset=")" << offset
4411 << "\"/>" << endl;
4412 offset += geomOffsetsMemsizeGlobal;
4413
4414 ofl << "</Polys>" << endl;
4415 //================== /Polys =================
4416
4417
4418 //================== PolyData =================
4419 ofl << "<CellData Scalars=\"scalars\">" << endl;
4420
4421 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="bodyId" format="appended" offset=")" << offset
4422 << "\"/>" << endl;
4423 offset += bodyIdMemsizeGlobal;
4424
4425 ofl << "<DataArray type=\"" << dataType
4426 << R"(" Name="velocity" NumberOfComponents="3" format="appended" offset=")" << offset << "\"/>" << endl;
4427 offset += velocityMemsizeGlobal;
4428
4429 if(m_solver->m_vtuGeometryOutputExtended) {
4430 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="domainId" format="appended" offset=")" << offset
4431 << "\"/>" << endl;
4432 offset += domainIdsMemsizeGlobal;
4433
4434 ofl << "<DataArray type=\"" << dataType << R"(" Name="pressure" format="appended" offset=")" << offset
4435 << "\"/>" << endl;
4436 offset += pressureMemsizeGlobal;
4437
4438 ofl << "<DataArray type=\"" << dataType << R"(" Name="density" format="appended" offset=")" << offset
4439 << "\"/>" << endl;
4440 offset += densityMemsizeGlobal;
4441
4442 ofl << "<DataArray type=\"" << dataType
4443 << R"(" Name="cell_normals" NumberOfComponents="3" format="appended" offset=")" << offset << "\"/>"
4444 << endl;
4445 offset += normalsMemsizeGlobal;
4446
4447 ofl << "<DataArray type=\"" << dataType
4448 << R"(" Name="shear" NumberOfComponents="3" format="appended" offset=")" << offset << "\"/>" << endl;
4449 offset += shearMemsizeGlobal;
4450 }
4451
4452 ofl << "</CellData>" << endl;
4453 //================== /PolyData =================
4454
4455
4456 ofl << "</Piece>" << endl;
4457 ofl << "</PolyData>" << endl;
4458
4459
4460 //================== AppendedData =================
4461 ofl << "<AppendedData encoding=\"raw\">" << endl;
4462 ofl << "_";
4463
4464 ofl.close();
4465 ofl.clear();
4466
4467 } else {
4468 cerr << "ERROR! COULD NOT OPEN FILE " << geomFileName << " for writing! (1)" << endl;
4469 }
4470 }
4471
4472
4473 //###################################################
4474 //################### parallel IO ###################
4475 MPI_File gfile = nullptr;
4476 MInt rcode = MPI_File_open(mpiComm(), const_cast<char*>(geomFileName.c_str()),
4477 MPI_MODE_APPEND | MPI_MODE_WRONLY, MPI_INFO_NULL, &gfile, AT_);
4478 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error opening file " + geomFileName);
4479
4480 //=================== Points =======================
4481 rcode =
4482 writeVtuArrayParallel(gfile, &geomPoints(0), geomPointsOffset, geomPointsMemsize, geomPointsMemsizeGlobal);
4483 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (1) writing to file " + geomFileName);
4484
4485 //=================== Connectivity =======================
4486 rcode = writeVtuArrayParallel(gfile, &geomConn(0), geomConnOffset, geomConnMemsize, geomConnMemsizeGlobal);
4487 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (2) writing to file " + geomFileName);
4488
4489 //=================== Offsets =======================
4490 rcode = writeVtuArrayParallel(gfile, &geomOffsets(0), geomOffsetsOffset, geomOffsetsMemsize,
4491 geomOffsetsMemsizeGlobal);
4492 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (3) writing to file " + geomFileName);
4493
4494 //=================== BodyId =======================
4495 rcode = writeVtuArrayParallel(gfile, &bodyId(0), bodyIdOffset, bodyIdMemsize, bodyIdMemsizeGlobal);
4496 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (4) writing to file " + geomFileName);
4497
4498 //=================== Velocity =======================
4499 rcode = writeVtuArrayParallel(gfile, &velocity(0), velocityOffset, velocityMemsize, velocityMemsizeGlobal);
4500 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (9) writing to file " + geomFileName);
4501
4502 if(m_solver->m_vtuGeometryOutputExtended) {
4503 //=================== DomainId =======================
4504 rcode =
4505 writeVtuArrayParallel(gfile, &domainIds(0), domainIdsOffset, domainIdsMemsize, domainIdsMemsizeGlobal);
4506 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (4) writing to file " + geomFileName);
4507
4508 //=================== Pressure =======================
4509 rcode = writeVtuArrayParallel(gfile, &pressure(0), pressureOffset, pressureMemsize, pressureMemsizeGlobal);
4510 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (5) writing to file " + geomFileName);
4511
4512 //=================== Density =======================
4513 rcode = writeVtuArrayParallel(gfile, &density(0), densityOffset, densityMemsize, densityMemsizeGlobal);
4514 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (6) writing to file " + geomFileName);
4515
4516 //=================== Normals =======================
4517 rcode = writeVtuArrayParallel(gfile, &normals(0), normalsOffset, normalsMemsize, normalsMemsizeGlobal);
4518 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (7) writing to file " + geomFileName);
4519
4520 //=================== Shear =======================
4521 rcode = writeVtuArrayParallel(gfile, &shear(0), shearOffset, shearMemsize, shearMemsizeGlobal);
4522 if(rcode != MPI_SUCCESS) mTerm(1, AT_, "Error (8) writing to file " + geomFileName);
4523 }
4524
4525 MPI_File_close(&gfile, AT_);
4526 //###################################################
4527 //###################################################
4528
4529
4530 if(domainId() == 0) {
4531 ofstream ofl;
4532 ofl.open(geomFileName.c_str(), ios_base::out | ios_base::app);
4533 if(ofl.is_open() && ofl.good()) {
4534 ofl << endl;
4535 ofl << "</AppendedData>" << endl;
4536 //================== /AppendedData =================
4537
4538 ofl << "</VTKFile>" << endl;
4539 ofl.close();
4540 ofl.clear();
4541 //================== /VTKFile =================
4542 } else {
4543 cerr << "ERROR! COULD NOT OPEN FILE " << geomFileName << " for writing! (3)" << endl;
4544 }
4545 }
4546
4547
4548 //-----------------------
4549 //------ GEOM.pvd -------
4550 //-----------------------
4551
4552 if(domainId() == 0) {
4553 if(m_solver->m_firstUseWriteVtuOutputParallelGeom) {
4554 m_solver->m_firstUseWriteVtuOutputParallelGeom = !m_solver->m_restart;
4555 }
4556 ofstream ofile;
4557 ofstream ofile2;
4558 if(m_solver->m_firstUseWriteVtuOutputParallelGeom)
4559 ofile.open("out/GEOM.pvd.tmp", ios_base::out | ios_base::trunc);
4560 else
4561 ofile.open("out/GEOM.pvd.tmp", ios_base::out | ios_base::app);
4562 if(ofile.is_open() && ofile.good()) {
4563 if(m_solver->m_firstUseWriteVtuOutputParallelGeom) {
4564 ofile << R"(<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">)" << endl;
4565 ofile << "<Collection>" << endl;
4566 }
4567 const size_t pos = geomFileName.find(m_solver->m_solutionOutput);
4568 MString fname = (pos == string::npos)
4569 ? geomFileName
4570 : geomFileName.substr(pos + m_solver->m_solutionOutput.size(), string::npos);
4571 ofile << "<DataSet part=\"" << 0 << "\" timestep=\"" << globalTimeStep << "\" file=\""
4572 << "./" << fname << "\"/>" << endl;
4573 ofile.close();
4574 ofile.clear();
4575 } else {
4576 cerr << "Error opening file out/GEOM.pvd.tmp" << endl;
4577 }
4578 m_log << "Command executed with return value: " << system("cp out/GEOM.pvd.tmp out/GEOM.pvd") << " at " << AT_
4579 << endl;
4580
4581 ofile2.open("out/GEOM.pvd", ios_base::out | ios_base::app);
4582 if(ofile2.is_open() && ofile2.good()) {
4583 ofile2 << "</Collection>" << endl;
4584 ofile2 << "</VTKFile>" << endl;
4585 ofile2.close();
4586 ofile2.clear();
4587 } else {
4588 cerr << "Error opening file out/GEOM.pvd" << endl;
4589 }
4590 m_solver->m_firstUseWriteVtuOutputParallelGeom = false;
4591 }
4592 }
4593
4594 if(domainId() == 0 && m_solver->m_levelSetMb && m_solver->m_vtuWriteParticleFile) {
4595 offset = 0;
4596 ofstream ofl;
4597 const MString partFileName = m_solver->m_solutionOutput + "POUT_00" + to_string(globalTimeStep) + ".vtp";
4598 ofl.open(partFileName.c_str(), ios_base::out | ios_base::trunc);
4599 if(ofl.is_open() && ofl.good()) {
4600 //================== VTKFile =================
4601 ofl << "<?xml version=\"1.0\"?>" << endl;
4602 ofl << R"(<VTKFile type="PolyData" version="1.0" byte_order="LittleEndian" header_type="UInt64">)" << endl;
4603 ofl << "<PolyData>" << endl;
4604
4605 //================== FieldData =================
4606 ofl << "<FieldData>" << endl;
4607 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="globalTimeStep" format="ascii" NumberOfTuples="1" > )"
4608 << globalTimeStep << " </DataArray>" << endl;
4609 if(m_solver->m_outputPhysicalTime) {
4610 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
4611 << m_solver->m_physicalTime << " </DataArray>" << endl;
4612 ofl << "<DataArray type=\"" << dataType << R"(" Name="internalTime" format="ascii" NumberOfTuples="1" > )"
4613 << m_solver->m_time << " </DataArray>" << endl;
4614 } else {
4615 ofl << "<DataArray type=\"" << dataType << R"(" Name="time" format="ascii" NumberOfTuples="1" > )"
4616 << m_solver->m_time << " </DataArray>" << endl;
4617 ofl << "<DataArray type=\"" << dataType << R"(" Name="physicalTime" format="ascii" NumberOfTuples="1" > )"
4618 << m_solver->m_physicalTime << " </DataArray>" << endl;
4619 }
4620 ofl << "</FieldData>" << endl;
4621 //================== /FieldData =================
4622
4623 //================== Dimensions =================
4624 ofl << "<Piece NumberOfPoints=\"" << m_solver->m_noEmbeddedBodies << "\" NumberOfVerts=\""
4625 << m_solver->m_noEmbeddedBodies << "\">" << endl;
4626 //================== /Dimensions ================
4627
4628 //================== Points =================
4629 ofl << "<Points>" << endl;
4630
4631 ofl << setprecision(12);
4632 ofl << "<DataArray type=\"" << dataType << R"(" NumberOfComponents="3" format="ascii">)" << endl;
4633 for(MInt k = 0; k < m_solver->m_noEmbeddedBodies; k++) {
4634 for(MInt i = 0; i < nDim; i++) {
4635 ofl << m_solver->m_bodyCenter[k * nDim + i] << " ";
4636 }
4637 ofl << endl;
4638 }
4639 ofl << "</DataArray>" << endl;
4640 ofl << "</Points>" << endl;
4641 //================== /Points =================
4642
4643
4644 //================== Verts =================
4645 ofl << "<Verts>" << endl;
4646 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="connectivity" format="ascii">)" << endl;
4647 for(MInt k = 0; k < m_solver->m_noEmbeddedBodies; k++) {
4648 ofl << k << " ";
4649 }
4650 ofl << endl;
4651 ofl << "</DataArray>" << endl;
4652 ofl << "<DataArray type=\"" << uIDataType << R"(" Name="offsets" format="ascii">)" << endl;
4653 for(MInt k = 0; k < m_solver->m_noEmbeddedBodies; k++) {
4654 ofl << k + 1 << " ";
4655 }
4656 ofl << endl;
4657 ofl << "</DataArray>" << endl;
4658 ofl << "</Verts>" << endl;
4659 //================== /Verts =================
4660
4661 //================== PointData =================
4662 ofl << "<PointData Scalars=\"scalars\">" << endl;
4663 ofl << "<DataArray type=\"" << dataType << R"(" Name="velocity" NumberOfComponents="3" format="ascii">)"
4664 << endl;
4665 for(MInt k = 0; k < m_solver->m_noEmbeddedBodies; k++) {
4666 for(MInt i = 0; i < nDim; i++) {
4667 ofl << m_solver->m_bodyVelocity[k * nDim + i] << " ";
4668 }
4669 ofl << endl;
4670 }
4671 ofl << "</DataArray>" << endl;
4672 ofl << "</PointData>" << endl;
4673 //================== /PointData =================
4674
4675 ofl << "</Piece>" << endl;
4676 ofl << "</PolyData>" << endl;
4677 ofl << "</VTKFile>" << endl;
4678 ofl.close();
4679 ofl.clear();
4680 //================== /VTKFile =================
4681 } else {
4682 cerr << "ERROR! COULD NOT OPEN FILE " << partFileName << " for writing! (1)" << endl;
4683 }
4684 }
4685 }
4686 }
4687
4688 if(facestream != nullptr) {
4689 delete[] facestream;
4690 facestream = nullptr;
4691 }
4692 if(conn != nullptr) {
4693 delete[] conn;
4694 conn = nullptr;
4695 }
4696
4697 RECORD_TIMER_STOP(t_geometry);
4698 RECORD_TIMER_STOP(t_vtutotal);
4699 DISPLAY_TIMER_INTERM(t_vtutotal);
4700}
FvCartesianSolverXD::m_splitChildToSplitCell
std::map< MInt, MInt > m_splitChildToSplitCell
Definition: fvcartesiansolverxd.h:2316
FvCartesianSolverXD::m_vtuGeometryOutput
std::set< MInt > m_vtuGeometryOutput
Definition: fvcartesiansolverxd.h:1675
InfoOutFile::open
void open(const MString &filename, const MString &projectName, MInt fileType=0, MPI_Comm mpiComm=MPI_COMM_WORLD, MBool rootOnlyHardwired=false)
Opens a file by passing the parameters to InfoOut_<xyz>FileBuffer::open(...).
Definition: infoout.cpp:975
VtkIo::insertDataHeader
void insertDataHeader(char *data, uint64_t &memsize, uint64_t &memsizeGlobal, uint64_t &offset)
The zeroth domain stores the header for an uncompressed vtu file appended data array specifying its n...
Definition: iovtk.cpp:1157
VtkIo::writeVtuArrayParallel
MInt writeVtuArrayParallel(MPI_File &, void *, MPI_Offset, MPI_Offset, MPI_Offset)
Definition: iovtk.cpp:134
VtkIo::vtuAssembleFaceStream
uint64_t vtuAssembleFaceStream(std::vector< T > *&, std::vector< T > *&, uint64_t &, ScratchSpace< MInt > &, MInt &, const MBool)
Definition: iovtk.cpp:1195
VtkIo::estimateMemorySizeSolverwise
MInt estimateMemorySizeSolverwise(uint64_t, ScratchSpace< uint64_t > &, uint64_t)
Definition: iovtk.cpp:1142
maia::CartesianSolver::noWindowCells
MInt noWindowCells(const MInt domainId) const
Definition: cartesiansolver.h:76
maia::CartesianSolver::domainOffset
MLong domainOffset(const MInt id) const
Definition: cartesiansolver.h:72
maia::CartesianSolver::extractPointIdsFromGrid
void extractPointIdsFromGrid(Collector< PointBasedCell< nDim > > *&, Collector< CartesianGridPoint< nDim > > *&, const MBool, const std::map< MInt, MInt > &, MInt levelThreshold=999999, MFloat *bBox=nullptr, MBool levelSetMb=false) const
Creates a list of unique corner points for all cells using a hash map levelThreshold optionally speci...
Definition: cartesiansolver.h:3524
maia::CartesianSolver::noHaloCells
MInt noHaloCells(const MInt domainId) const
Definition: cartesiansolver.h:74
maia::CartesianSolver::windowCellId
MInt windowCellId(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:77
maia::CartesianSolver::noNeighborDomains
MInt noNeighborDomains() const
Definition: cartesiansolver.h:70
maia::CartesianSolver::haloCellId
MInt haloCellId(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:75
maia::CartesianSolver::neighborDomain
MInt neighborDomain(const MInt id) const
Definition: cartesiansolver.h:71
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
VTP
@ VTP
Definition: enums.h:18
VTU
@ VTU
Definition: enums.h:18
MPI_Recv
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status, const MString &name, const MString &varname)
same as MPI_Recv
Definition: mpioverride.cpp:565
MPI_Issend
int MPI_Issend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Issend
Definition: mpioverride.cpp:523
MPI_File_open
int MPI_File_open(MPI_Comm comm, const char *filename, int amode, MPI_Info info, MPI_File *mpi_fh, const MString &name)
same as MPI_File_open
Definition: mpioverride.cpp:2203
MPI_Waitall
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
Definition: mpioverride.cpp:766
MPI_File_close
int MPI_File_close(MPI_File *mpi_fh, const MString &name)
same as MPI_File_close
Definition: mpioverride.cpp:2273
maia::math::calcEigenValues
void calcEigenValues(MFloat A[3][3], MFloat w[3])
Definition: maiamath.cpp:475

Member Data Documentation

◆ m_bndryCellIds

template<MInt nDim, class SysEqn >
List<MInt>* VtkIo< nDim, SysEqn >::m_bndryCellIds = nullptr
private

Definition at line 30 of file iovtk.h.

◆ m_bndryCells

template<MInt nDim, class SysEqn >
Collector<FvBndryCell<nDim, SysEqn> >* VtkIo< nDim, SysEqn >::m_bndryCells = nullptr
private

Definition at line 27 of file iovtk.h.

◆ m_extractedCells

template<MInt nDim, class SysEqn >
Collector<PointBasedCell<nDim> >* VtkIo< nDim, SysEqn >::m_extractedCells = nullptr
private

Definition at line 28 of file iovtk.h.

◆ m_gridPoints

template<MInt nDim, class SysEqn >
Collector<CartesianGridPoint<nDim> >* VtkIo< nDim, SysEqn >::m_gridPoints = nullptr
private

Definition at line 29 of file iovtk.h.

◆ m_noDirs

template<MInt nDim, class SysEqn >
constexpr MInt VtkIo< nDim, SysEqn >::m_noDirs = 2 * nDim
staticconstexprprivate

Definition at line 34 of file iovtk.h.

◆ m_noSolver

template<MInt nDim, class SysEqn >
MInt VtkIo< nDim, SysEqn >::m_noSolver
private

Definition at line 32 of file iovtk.h.

◆ m_solver

template<MInt nDim, class SysEqn >
FvCartesianSolverXD<nDim, SysEqn>* VtkIo< nDim, SysEqn >::m_solver = nullptr
private

Definition at line 26 of file iovtk.h.

The documentation for this class was generated from the following files: