IsoExtractorMPI.h (9556B)
1 // File : IsoExtractorMPI.h 2 // Date : Tue Nov 29 10:18:26 2016 3 // Author : Fabian Wermelinger 4 // Description: MPI'ed isosurface extractor 5 // Copyright 2016 ETH Zurich. All Rights Reserved. 6 #ifndef ISOEXTRACTORMPI_H_XABTB0WN 7 #define ISOEXTRACTORMPI_H_XABTB0WN 8 9 #include <mpi.h> 10 #include <cstring> 11 #include "IsoExtractor.h" 12 #include "InterpolatorMPI.h" 13 14 struct TriangleConnect 15 { 16 unsigned char n; 17 int v[3]; 18 } __attribute__((packed)); 19 20 template <typename Tkernel, template<typename> class Tinterpolator=InterpolatorMPI> 21 class IsoExtractorMPI : public IsoExtractor<Tkernel,Tinterpolator> 22 { 23 public: 24 IsoExtractorMPI(ArgumentParser& parser) : IsoExtractor<Tkernel,Tinterpolator>(parser) {} 25 ~IsoExtractorMPI() {} 26 27 void extract(const Real isoval, const std::string fout) 28 { 29 const Real cscale = this->m_parser("cube_scale").asDouble(1.0); 30 assert(cscale > 0.0); 31 int gNx = static_cast<Real>(this->m_interp.getNx()) / cscale; 32 int gNy = static_cast<Real>(this->m_interp.getNy()) / cscale; 33 int gNz = static_cast<Real>(this->m_interp.getNz()) / cscale; 34 int pesize[3]; 35 this->m_interp.getPESize(pesize); 36 gNx += (gNx % pesize[0]); 37 gNy += (gNy % pesize[1]); 38 gNz += (gNz % pesize[2]); 39 const int Nx = gNx/pesize[0]; 40 const int Ny = gNy/pesize[1]; 41 const int Nz = gNz/pesize[2]; 42 assert(Nx > 0); 43 assert(Ny > 0); 44 assert(Nz > 0); 45 const Real chx = this->m_interp.getExtentX() / gNx; 46 const Real chy = this->m_interp.getExtentY() / gNy; 47 const Real chz = this->m_interp.getExtentZ() / gNz; 48 49 int peidx[3]; 50 this->m_interp.getPEIndex(peidx); 51 const Real Ox = (this->m_interp.getExtentX()/pesize[0])*peidx[0]; 52 const Real Oy = (this->m_interp.getExtentY()/pesize[1])*peidx[1]; 53 const Real Oz = (this->m_interp.getExtentZ()/pesize[2])*peidx[2]; 54 const Real data_invh = 1.0 / this->m_interp.getH(); 55 Real origin[3]; 56 origin[0] = this->m_interp.getDataPos(ceil(Ox*data_invh)) - Ox; 57 origin[1] = this->m_interp.getDataPos(ceil(Oy*data_invh)) - Oy; 58 origin[2] = this->m_interp.getDataPos(ceil(Oz*data_invh)) - Oz; 59 assert(origin[0] >= 0.0); 60 assert(origin[1] >= 0.0); 61 assert(origin[2] >= 0.0); 62 this->m_interp.setOrigin(origin); 63 64 this->_extractIso(isoval, Nx, Ny, Nz, chx, chy, chz, Ox, Oy, Oz, this->m_interp.isroot()); 65 MPI_Barrier(this->m_interp.getComm()); 66 67 _writePLY_MPI(fout, this->m_triList); 68 } 69 70 private: 71 72 void _writePLY_MPI(const std::string basename, const std::vector<Triangle>& tList) 73 { 74 const MPI_Comm comm = this->m_interp.getComm(); 75 unsigned long long mytriangles = tList.size(); 76 unsigned long long alltriangles; 77 MPI_Reduce(&mytriangles, &alltriangles, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, comm); 78 79 const bool isroot = this->m_interp.isroot(); 80 81 time_t rawtime; 82 std::time(&rawtime); 83 struct tm* timeinfo = std::localtime(&rawtime); 84 char buf[256]; 85 std::strftime(buf, 256, "%A, %h %d %Y, %r %Z %z", timeinfo); 86 87 const unsigned int one = 1; 88 const unsigned char* const lsb = (const unsigned char*)(&one); 89 const std::string endianess((*lsb) ? "binary_little_endian" : "binary_big_endian"); 90 91 // write header 92 if (isroot) 93 { 94 std::ofstream ply(basename+".ply"); 95 ply << "ply" << std::endl; 96 ply << "format "; 97 if (this->m_parser("binary").asBool(true)) 98 ply << endianess; 99 else 100 ply << "ascii"; 101 ply << " 1.0" << std::endl; 102 ply << "comment generated by IsoExtractor.h on " << buf << std::endl; 103 ply << "element vertex " << alltriangles*3 << std::endl; 104 ply << "property float x" << std::endl; 105 ply << "property float y" << std::endl; 106 ply << "property float z" << std::endl; 107 ply << "element face " << alltriangles << std::endl; 108 ply << "property list uchar int vertex_index" << std::endl; 109 ply << "end_header" << std::endl; 110 ply.close(); 111 } 112 MPI_Barrier(comm); 113 114 // TODO: (fabianw@mavt.ethz.ch; Tue Nov 29 14:04:45 2016) write proper 115 // MPI I/O 116 // write content 117 // _writePLY_binary_MPI(comm, tList); 118 119 // _writePLY_binary_MPI_lazy(tList, basename+".ply"); 120 _writePLY_binary_MPI(comm, tList, basename+".ply"); 121 } 122 123 void _writePLY_binary_MPI_lazy(const std::vector<Triangle>& tList, const std::string fname) 124 { 125 const MPI_Comm comm = this->m_interp.getComm(); 126 int size, rank; 127 MPI_Comm_size(comm, &size); 128 MPI_Comm_rank(comm, &rank); 129 for (int r = 0; r < size; ++r) 130 { 131 if (r == rank) 132 { 133 std::ofstream ply(fname, std::ios::out | std::ios::app | std::ios::binary); 134 for (auto t : tList) 135 { 136 float tc[3][3]; 137 for (int j = 0; j < 3; ++j) 138 { 139 tc[j][0] = t.vertices[j][0]; 140 tc[j][1] = t.vertices[j][1]; 141 tc[j][2] = t.vertices[j][2]; 142 } 143 ply.write((char*)&tc[0][0], 9*sizeof(float)); 144 } 145 ply.close(); 146 } 147 MPI_Barrier(comm); 148 } 149 150 const char nVert = 3; 151 std::vector<unsigned long long> allsize(size); 152 unsigned long long mytri = tList.size(); 153 MPI_Allgather(&mytri, 1, MPI_UNSIGNED_LONG_LONG, allsize.data(), 1, MPI_UNSIGNED_LONG_LONG, comm); 154 for (int r = 0; r < size; ++r) 155 { 156 if (r == rank) 157 { 158 unsigned long long start = 0; 159 for (int i = 0; i < rank; ++i) 160 start += allsize[i]; 161 162 std::ofstream ply(fname, std::ios::out | std::ios::app | std::ios::binary); 163 for (size_t i = 3*start; i < 3*start+tList.size()*3; i += 3) 164 { 165 int tv[3] = {i, i+1, i+2}; 166 ply.write(&nVert, sizeof(char)); 167 ply.write((char*)&tv[0], 3*sizeof(int)); 168 } 169 ply.close(); 170 } 171 MPI_Barrier(comm); 172 } 173 } 174 175 void _writePLY_binary_MPI(const MPI_Comm comm, const std::vector<Triangle>& tList, const std::string fname) 176 { 177 int rank, size; 178 MPI_Status status; 179 MPI_File file_id; 180 MPI_Comm_rank(comm, &rank); 181 MPI_Comm_size(comm, &size); 182 183 int rc = MPI_File_open(MPI_COMM_SELF, (char*)fname.c_str(), MPI_MODE_APPEND | MPI_MODE_WRONLY, MPI_INFO_NULL, &file_id); 184 if (rc) 185 { 186 if (0 == rank) 187 std::cerr << "ERROR: IsoExtractorMPI: Can not open file." << std::endl; 188 abort(); 189 } 190 191 std::vector<unsigned long long> allsize(size); 192 unsigned long long mytri = tList.size(); 193 MPI_Allgather(&mytri, 1, MPI_UNSIGNED_LONG_LONG, allsize.data(), 1, MPI_UNSIGNED_LONG_LONG, comm); 194 195 long long offset = 0; 196 MPI_File_get_position(file_id, &offset); 197 198 assert(sizeof(Triangle) == 9*sizeof(float)); 199 for (int i = 0; i < rank; ++i) 200 offset = offset + allsize[i]*sizeof(Triangle); 201 202 Triangle* tarr = nullptr; 203 if (allsize[rank] > 0) 204 { 205 tarr = new Triangle[allsize[rank]]; 206 char* const base = (char*)&tarr[0].vertices[0][0]; 207 for (size_t i=0; i<allsize[rank]; ++i) 208 { 209 char* const dst = base + i*sizeof(Triangle); 210 const char* const src = (char*)&tList[i].vertices[0][0]; 211 memcpy(dst, src, sizeof(Triangle)); 212 } 213 // for (size_t i=0; i<allsize[rank]; ++i) 214 // { 215 // for (int j = 0; j < 3; ++j) 216 // { 217 // tarr[i].vertices[j][0] = tList[i].vertices[j][0]; 218 // tarr[i].vertices[j][1] = tList[i].vertices[j][1]; 219 // tarr[i].vertices[j][2] = tList[i].vertices[j][2]; 220 // } 221 // } 222 223 MPI_File_write_at(file_id, offset, (char *)&tarr[0].vertices[0][0], allsize[rank]*sizeof(Triangle), MPI_CHAR, &status); 224 } 225 offset = offset + allsize[rank]*sizeof(Triangle); 226 MPI_Barrier(comm); 227 MPI_Bcast(&offset, 1, MPI_LONG_LONG, size-1, comm); 228 229 if (tarr != nullptr) 230 delete [] tarr; 231 232 for (int i = 0; i < rank; ++i) 233 offset = offset + allsize[i]*sizeof(TriangleConnect); 234 235 TriangleConnect* tcarr = nullptr; 236 if (allsize[rank] > 0) 237 { 238 tcarr = new TriangleConnect[allsize[rank]]; 239 int base = 0; 240 for (int i = 0; i < rank; ++i) 241 base += 3*allsize[i]; 242 243 for (size_t i=0; i<allsize[rank]; ++i) 244 { 245 tcarr[i].n = 3; 246 tcarr[i].v[0] = base + 3*i + 0; 247 tcarr[i].v[1] = base + 3*i + 1; 248 tcarr[i].v[2] = base + 3*i + 2; 249 } 250 251 MPI_File_write_at(file_id, offset, (char *)&tcarr[0].n, allsize[rank]*sizeof(TriangleConnect), MPI_CHAR, &status); 252 } 253 254 if (tcarr != nullptr) 255 delete [] tcarr; 256 257 MPI_File_close(&file_id); 258 } 259 }; 260 261 #endif /* ISOEXTRACTORMPI_H_XABTB0WN */