surface.hpp 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395
  1. #ifndef _SURFACE_HPP_
  2. #define _SURFACE_HPP_
  3. #include <sctl.hpp>
  4. namespace biest {
  5. enum class SurfType {
  6. AxisymCircleWide, // 125 x 50
  7. AxisymCircleNarrow, // 250 x 50
  8. AxisymWide, // 125 x 50
  9. AxisymNarrow, // 250 x 50
  10. RotatingEllipseWide, // 125 x 50
  11. RotatingEllipseNarrow, // 250 x 50
  12. Quas3, // 250 x 50
  13. LHD, // 250 x 50
  14. W7X, // 250 x 45
  15. Stell, // 250 x 50
  16. W7X_ // 250 x 50
  17. };
  18. template <class Real> class Surface {
  19. static constexpr sctl::Integer COORD_DIM = 3;
  20. public:
  21. Surface() : Nt0(0), Np0(0) {}
  22. Surface(sctl::Long Nt, sctl::Long Np, SurfType type = SurfType::AxisymCircleWide);
  23. sctl::Long NTor() const {return Nt0;}
  24. sctl::Long NPol() const {return Np0;}
  25. sctl::Vector<Real>& Coord() {return X0_;}
  26. const sctl::Vector<Real>& Coord() const {return X0_;}
  27. Real& Coord(sctl::Long i, sctl::Long j, sctl::Long k) {return X0_[(k*Nt0+i)*Np0+j];}
  28. const Real& Coord(sctl::Long i, sctl::Long j, sctl::Long k) const {return X0_[(k*Nt0+i)*Np0+j];}
  29. private:
  30. sctl::Long Nt0, Np0;
  31. sctl::Vector<Real> X0_;
  32. };
  33. struct VTKData {
  34. typedef double VTKReal;
  35. static constexpr sctl::Integer COORD_DIM = 3;
  36. std::vector<VTKReal> point_coord;
  37. std::vector<VTKReal> point_value;
  38. std::vector<int32_t> line_connect;
  39. std::vector<int32_t> line_offset;
  40. std::vector<int32_t> poly_connect;
  41. std::vector<int32_t> poly_offset;
  42. void WriteVTK(const char* fname, const sctl::Comm& comm = sctl::Comm::Self()) {
  43. sctl::Integer np = comm.Size();
  44. sctl::Integer myrank = comm.Rank();
  45. std::vector<VTKReal>& coord=this->point_coord;
  46. std::vector<VTKReal>& value=this->point_value;
  47. std::vector<int32_t>& line_connect=this->line_connect;
  48. std::vector<int32_t>& line_offset=this->line_offset;
  49. std::vector<int32_t>& poly_connect=this->poly_connect;
  50. std::vector<int32_t>& poly_offset=this->poly_offset;
  51. sctl::Long pt_cnt=coord.size()/COORD_DIM;
  52. sctl::Long line_cnt=line_offset.size();
  53. sctl::Long poly_cnt=poly_offset.size();
  54. // Open file for writing.
  55. std::stringstream vtufname;
  56. vtufname<<fname<<"_"<<std::setfill('0')<<std::setw(6)<<myrank<<".vtp";
  57. std::ofstream vtufile;
  58. vtufile.open(vtufname.str().c_str());
  59. if(vtufile.fail()) return;
  60. bool isLittleEndian;
  61. { // Set isLittleEndian
  62. uint16_t number = 0x1;
  63. uint8_t *numPtr = (uint8_t*)&number;
  64. isLittleEndian=(numPtr[0] == 1);
  65. }
  66. // Proceed to write to file.
  67. sctl::Long data_size=0;
  68. vtufile<<"<?xml version=\"1.0\"?>\n";
  69. if(isLittleEndian) vtufile<<"<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
  70. else vtufile<<"<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"BigEndian\">\n";
  71. //===========================================================================
  72. vtufile<<" <PolyData>\n";
  73. vtufile<<" <Piece NumberOfPoints=\""<<pt_cnt<<"\" NumberOfVerts=\"0\" NumberOfLines=\""<<line_cnt<<"\" NumberOfStrips=\"0\" NumberOfPolys=\""<<poly_cnt<<"\">\n";
  74. //---------------------------------------------------------------------------
  75. vtufile<<" <Points>\n";
  76. vtufile<<" <DataArray type=\"Float"<<sizeof(VTKReal)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  77. data_size+=sizeof(uint32_t)+coord.size()*sizeof(VTKReal);
  78. vtufile<<" </Points>\n";
  79. //---------------------------------------------------------------------------
  80. if(value.size()){ // value
  81. vtufile<<" <PointData>\n";
  82. vtufile<<" <DataArray type=\"Float"<<sizeof(VTKReal)*8<<"\" NumberOfComponents=\""<<value.size()/pt_cnt<<"\" Name=\""<<"value"<<"\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  83. data_size+=sizeof(uint32_t)+value.size()*sizeof(VTKReal);
  84. vtufile<<" </PointData>\n";
  85. }
  86. //---------------------------------------------------------------------------
  87. vtufile<<" <Lines>\n";
  88. vtufile<<" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  89. data_size+=sizeof(uint32_t)+line_connect.size()*sizeof(int32_t);
  90. vtufile<<" <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  91. data_size+=sizeof(uint32_t)+line_offset.size() *sizeof(int32_t);
  92. vtufile<<" </Lines>\n";
  93. //---------------------------------------------------------------------------
  94. vtufile<<" <Polys>\n";
  95. vtufile<<" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  96. data_size+=sizeof(uint32_t)+poly_connect.size()*sizeof(int32_t);
  97. vtufile<<" <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\""<<data_size<<"\" />\n";
  98. data_size+=sizeof(uint32_t)+poly_offset.size() *sizeof(int32_t);
  99. vtufile<<" </Polys>\n";
  100. //---------------------------------------------------------------------------
  101. vtufile<<" </Piece>\n";
  102. vtufile<<" </PolyData>\n";
  103. //===========================================================================
  104. vtufile<<" <AppendedData encoding=\"raw\">\n";
  105. vtufile<<" _";
  106. int32_t block_size;
  107. block_size=coord.size()*sizeof(VTKReal); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&coord [0], coord.size()*sizeof(VTKReal));
  108. if(value.size()){ // value
  109. block_size=value.size()*sizeof(VTKReal); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&value [0], value.size()*sizeof(VTKReal));
  110. }
  111. block_size=line_connect.size()*sizeof(int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&line_connect[0], line_connect.size()*sizeof(int32_t));
  112. block_size=line_offset .size()*sizeof(int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&line_offset [0], line_offset .size()*sizeof(int32_t));
  113. block_size=poly_connect.size()*sizeof(int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&poly_connect[0], poly_connect.size()*sizeof(int32_t));
  114. block_size=poly_offset .size()*sizeof(int32_t); vtufile.write((char*)&block_size, sizeof(int32_t)); vtufile.write((char*)&poly_offset [0], poly_offset .size()*sizeof(int32_t));
  115. vtufile<<"\n";
  116. vtufile<<" </AppendedData>\n";
  117. //===========================================================================
  118. vtufile<<"</VTKFile>\n";
  119. vtufile.close();
  120. if(myrank) return;
  121. std::stringstream pvtufname;
  122. pvtufname<<fname<<".pvtp";
  123. std::ofstream pvtufile;
  124. pvtufile.open(pvtufname.str().c_str());
  125. if(pvtufile.fail()) return;
  126. pvtufile<<"<?xml version=\"1.0\"?>\n";
  127. pvtufile<<"<VTKFile type=\"PPolyData\">\n";
  128. pvtufile<<" <PPolyData GhostLevel=\"0\">\n";
  129. pvtufile<<" <PPoints>\n";
  130. pvtufile<<" <PDataArray type=\"Float"<<sizeof(VTKReal)*8<<"\" NumberOfComponents=\""<<COORD_DIM<<"\" Name=\"Position\"/>\n";
  131. pvtufile<<" </PPoints>\n";
  132. if(value.size()){ // value
  133. pvtufile<<" <PPointData>\n";
  134. pvtufile<<" <PDataArray type=\"Float"<<sizeof(VTKReal)*8<<"\" NumberOfComponents=\""<<value.size()/pt_cnt<<"\" Name=\""<<"value"<<"\"/>\n";
  135. pvtufile<<" </PPointData>\n";
  136. }
  137. {
  138. // Extract filename from path.
  139. std::stringstream vtupath;
  140. vtupath<<'/'<<fname;
  141. std::string pathname = vtupath.str();
  142. auto found = pathname.find_last_of("/\\");
  143. std::string fname_ = pathname.substr(found+1);
  144. for(sctl::Integer i=0;i<np;i++) pvtufile<<" <Piece Source=\""<<fname_<<"_"<<std::setfill('0')<<std::setw(6)<<i<<".vtp\"/>\n";
  145. }
  146. pvtufile<<" </PPolyData>\n";
  147. pvtufile<<"</VTKFile>\n";
  148. pvtufile.close();
  149. }
  150. };
  151. struct VTUData {
  152. typedef float VTKReal;
  153. // Point data
  154. sctl::Vector<VTKReal> coord; // always 3D
  155. sctl::Vector<VTKReal> value;
  156. // Cell data
  157. sctl::Vector<int32_t> connect;
  158. sctl::Vector<int32_t> offset;
  159. sctl::Vector<uint8_t> types;
  160. void WriteVTK(const std::string& fname, const sctl::Comm& comm = sctl::Comm::Self()) const {
  161. typedef typename VTUData::VTKReal VTKReal;
  162. sctl::Long value_dof = 0;
  163. { // Write vtu file.
  164. std::ofstream vtufile;
  165. { // Open file for writing.
  166. std::stringstream vtufname;
  167. vtufname << fname << std::setfill('0') << std::setw(6) << comm.Rank() << ".vtu";
  168. vtufile.open(vtufname.str().c_str());
  169. if (vtufile.fail()) return;
  170. }
  171. { // Write to file.
  172. sctl::Long pt_cnt = coord.Dim() / 3;
  173. sctl::Long cell_cnt = types.Dim();
  174. value_dof = (pt_cnt ? value.Dim() / pt_cnt : 0);
  175. sctl::Vector<int32_t> mpi_rank;
  176. { // Set mpi_rank
  177. sctl::Integer new_myrank = comm.Rank();
  178. mpi_rank.ReInit(pt_cnt);
  179. for (sctl::Long i = 0; i < mpi_rank.Dim(); i++) mpi_rank[i] = new_myrank;
  180. }
  181. bool isLittleEndian;
  182. { // Set isLittleEndian
  183. uint16_t number = 0x1;
  184. uint8_t *numPtr = (uint8_t *)&number;
  185. isLittleEndian = (numPtr[0] == 1);
  186. }
  187. sctl::Long data_size = 0;
  188. vtufile << "<?xml version=\"1.0\"?>\n";
  189. vtufile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"" << (isLittleEndian ? "LittleEndian" : "BigEndian") << "\">\n";
  190. // ===========================================================================
  191. vtufile << " <UnstructuredGrid>\n";
  192. vtufile << " <Piece NumberOfPoints=\"" << pt_cnt << "\" NumberOfCells=\"" << cell_cnt << "\">\n";
  193. //---------------------------------------------------------------------------
  194. vtufile << " <Points>\n";
  195. vtufile << " <DataArray type=\"Float" << sizeof(VTKReal) * 8 << "\" NumberOfComponents=\"3\" Name=\"Position\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  196. data_size += sizeof(uint32_t) + coord.Dim() * sizeof(VTKReal);
  197. vtufile << " </Points>\n";
  198. //---------------------------------------------------------------------------
  199. vtufile << " <PointData>\n";
  200. if (value_dof) { // value
  201. vtufile << " <DataArray type=\"Float" << sizeof(VTKReal) * 8 << "\" NumberOfComponents=\"" << value_dof << "\" Name=\"value\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  202. data_size += sizeof(uint32_t) + value.Dim() * sizeof(VTKReal);
  203. }
  204. { // mpi_rank
  205. vtufile << " <DataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  206. data_size += sizeof(uint32_t) + pt_cnt * sizeof(int32_t);
  207. }
  208. vtufile << " </PointData>\n";
  209. //---------------------------------------------------------------------------
  210. //---------------------------------------------------------------------------
  211. vtufile << " <Cells>\n";
  212. vtufile << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  213. data_size += sizeof(uint32_t) + connect.Dim() * sizeof(int32_t);
  214. vtufile << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  215. data_size += sizeof(uint32_t) + offset.Dim() * sizeof(int32_t);
  216. vtufile << " <DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"" << data_size << "\" />\n";
  217. data_size += sizeof(uint32_t) + types.Dim() * sizeof(uint8_t);
  218. vtufile << " </Cells>\n";
  219. //---------------------------------------------------------------------------
  220. vtufile << " </Piece>\n";
  221. vtufile << " </UnstructuredGrid>\n";
  222. // ===========================================================================
  223. vtufile << " <AppendedData encoding=\"raw\">\n";
  224. vtufile << " _";
  225. int32_t block_size;
  226. { // coord
  227. block_size = coord.Dim() * sizeof(VTKReal);
  228. vtufile.write((char *)&block_size, sizeof(int32_t));
  229. if (coord.Dim()) vtufile.write((char *)&coord[0], coord.Dim() * sizeof(VTKReal));
  230. }
  231. if (value_dof) { // value
  232. block_size = value.Dim() * sizeof(VTKReal);
  233. vtufile.write((char *)&block_size, sizeof(int32_t));
  234. if (value.Dim()) vtufile.write((char *)&value[0], value.Dim() * sizeof(VTKReal));
  235. }
  236. { // mpi_rank
  237. block_size = mpi_rank.Dim() * sizeof(int32_t);
  238. vtufile.write((char *)&block_size, sizeof(int32_t));
  239. if (mpi_rank.Dim()) vtufile.write((char *)&mpi_rank[0], mpi_rank.Dim() * sizeof(int32_t));
  240. }
  241. { // block_size
  242. block_size = connect.Dim() * sizeof(int32_t);
  243. vtufile.write((char *)&block_size, sizeof(int32_t));
  244. if (connect.Dim()) vtufile.write((char *)&connect[0], connect.Dim() * sizeof(int32_t));
  245. }
  246. { // offset
  247. block_size = offset.Dim() * sizeof(int32_t);
  248. vtufile.write((char *)&block_size, sizeof(int32_t));
  249. if (offset.Dim()) vtufile.write((char *)&offset[0], offset.Dim() * sizeof(int32_t));
  250. }
  251. { // types
  252. block_size = types.Dim() * sizeof(uint8_t);
  253. vtufile.write((char *)&block_size, sizeof(int32_t));
  254. if (types.Dim()) vtufile.write((char *)&types[0], types.Dim() * sizeof(uint8_t));
  255. }
  256. vtufile << "\n";
  257. vtufile << " </AppendedData>\n";
  258. // ===========================================================================
  259. vtufile << "</VTKFile>\n";
  260. }
  261. vtufile.close(); // close file
  262. }
  263. if (!comm.Rank()) { // Write pvtu file
  264. std::ofstream pvtufile;
  265. { // Open file for writing
  266. std::stringstream pvtufname;
  267. pvtufname << fname << ".pvtu";
  268. pvtufile.open(pvtufname.str().c_str());
  269. if (pvtufile.fail()) return;
  270. }
  271. { // Write to file.
  272. pvtufile << "<?xml version=\"1.0\"?>\n";
  273. pvtufile << "<VTKFile type=\"PUnstructuredGrid\">\n";
  274. pvtufile << " <PUnstructuredGrid GhostLevel=\"0\">\n";
  275. pvtufile << " <PPoints>\n";
  276. pvtufile << " <PDataArray type=\"Float" << sizeof(VTKReal) * 8 << "\" NumberOfComponents=\"3\" Name=\"Position\"/>\n";
  277. pvtufile << " </PPoints>\n";
  278. pvtufile << " <PPointData>\n";
  279. if (value_dof) { // value
  280. pvtufile << " <PDataArray type=\"Float" << sizeof(VTKReal) * 8 << "\" NumberOfComponents=\"" << value_dof << "\" Name=\"value\"/>\n";
  281. }
  282. { // mpi_rank
  283. pvtufile << " <PDataArray type=\"Int32\" NumberOfComponents=\"1\" Name=\"mpi_rank\"/>\n";
  284. }
  285. pvtufile << " </PPointData>\n";
  286. {
  287. // Extract filename from path.
  288. std::stringstream vtupath;
  289. vtupath << '/' << fname;
  290. std::string pathname = vtupath.str();
  291. std::string fname_ = pathname.substr(pathname.find_last_of("/\\") + 1);
  292. // char *fname_ = (char*)strrchr(vtupath.str().c_str(), '/') + 1;
  293. // std::string fname_ =
  294. // boost::filesystem::path(fname).filename().string().
  295. for (sctl::Integer i = 0; i < comm.Size(); i++) pvtufile << " <Piece Source=\"" << fname_ << std::setfill('0') << std::setw(6) << i << ".vtu\"/>\n";
  296. }
  297. pvtufile << " </PUnstructuredGrid>\n";
  298. pvtufile << "</VTKFile>\n";
  299. }
  300. pvtufile.close(); // close file
  301. }
  302. };
  303. };
  304. template <class Real> void WriteVTK(const char* fname, const sctl::Vector<Surface<Real>>& Svec, const sctl::Vector<Real> F = sctl::Vector<Real>(), const sctl::Comm& comm = sctl::Comm::Self()) {
  305. VTKData data;
  306. typedef VTKData::VTKReal VTKReal;
  307. auto& point_coord =data.point_coord ;
  308. auto& point_value =data.point_value ;
  309. auto& poly_connect=data.poly_connect;
  310. auto& poly_offset =data.poly_offset ;
  311. constexpr sctl::Integer COORD_DIM = VTKData::COORD_DIM;
  312. sctl::Long dof, offset = 0;
  313. { // Set dof
  314. sctl::Long Npt = 0;
  315. for (const auto& S : Svec) Npt += S.NTor() * S.NPol();
  316. dof = F.Dim() / Npt;
  317. SCTL_ASSERT(F.Dim() == dof * Npt);
  318. }
  319. for (auto& S: Svec) { // Set point_coord, point_value, poly_connect
  320. sctl::Long Nt = S.NTor();
  321. sctl::Long Np = S.NPol();
  322. sctl::Long N = Nt * Np;
  323. for (sctl::Long i = 0; i < Nt; i++) {
  324. for (sctl::Long j = 0; j < Np; j++) {
  325. sctl::Long i0 = (i + 0) % Nt;
  326. sctl::Long i1 = (i + 1) % Nt;
  327. sctl::Long j0 = (j + 0) % Np;
  328. sctl::Long j1 = (j + 1) % Np;
  329. poly_connect.push_back(point_coord.size() / COORD_DIM + Np*i0+j0);
  330. poly_connect.push_back(point_coord.size() / COORD_DIM + Np*i1+j0);
  331. poly_connect.push_back(point_coord.size() / COORD_DIM + Np*i1+j1);
  332. poly_connect.push_back(point_coord.size() / COORD_DIM + Np*i0+j1);
  333. poly_offset.push_back(poly_connect.size());
  334. }
  335. }
  336. const auto X = S.Coord();
  337. SCTL_ASSERT(X.Dim() == COORD_DIM * N);
  338. for (sctl::Long i = 0; i < N; i++){ // Set point_coord
  339. for (sctl::Integer k = 0; k < COORD_DIM; k++) point_coord.push_back((VTKReal)X[k * N + i]);
  340. }
  341. for (sctl::Long i = 0; i < N; i++){ // Set point_value
  342. for (sctl::Long k = 0; k < dof; k++) point_value.push_back((VTKReal)F[dof * offset + k * N + i]);
  343. }
  344. offset += N;
  345. }
  346. data.WriteVTK(fname, comm);
  347. }
  348. template <class Real> void WriteVTK(const char* fname, const Surface<Real>& S, const sctl::Vector<Real> F = sctl::Vector<Real>(), const sctl::Comm& comm = sctl::Comm::Self()) {
  349. WriteVTK(fname, sctl::Vector<Surface<Real>>(1,sctl::Ptr2Itr<Surface<Real>>((Surface<Real>*)&S,1),false), F, comm);
  350. }
  351. }
  352. #include <biest/surface.txx>
  353. #endif //_SURFACE_HPP_