IsoExtractorMPI.h 9.33 KB
Newer Older
1
2
3
4
5
6
7
8
9
// File       : IsoExtractorMPI.h
// Date       : Tue Nov 29 10:18:26 2016
// Author     : Fabian Wermelinger
// Description: MPI'ed isosurface extractor
// Copyright 2016 ETH Zurich. All Rights Reserved.
#ifndef ISOEXTRACTORMPI_H_XABTB0WN
#define ISOEXTRACTORMPI_H_XABTB0WN

#include <mpi.h>
10
#include <cstring>
11
12
13
#include "IsoExtractor.h"
#include "InterpolatorMPI.h"

fabianw's avatar
fabianw committed
14
15
16
17
18
19
struct TriangleConnect
{
    unsigned char n;
    int v[3];
} __attribute__((packed));

20
21
22
23
24
25
26
27
28
29
30
template <typename Tkernel, template<typename> class Tinterpolator=InterpolatorMPI>
class IsoExtractorMPI : public IsoExtractor<Tkernel,Tinterpolator>
{
public:
    IsoExtractorMPI(ArgumentParser& parser) : IsoExtractor<Tkernel,Tinterpolator>(parser) {}
    ~IsoExtractorMPI() {}

    void extract(const Real isoval, const std::string fout)
    {
        const Real cscale = this->m_parser("cube_scale").asDouble(1.0);
        assert(cscale > 0.0);
fabianw's avatar
fabianw committed
31
32
33
        int gNx = static_cast<Real>(this->m_interp.getNx()) / cscale;
        int gNy = static_cast<Real>(this->m_interp.getNy()) / cscale;
        int gNz = static_cast<Real>(this->m_interp.getNz()) / cscale;
34
35
        int pesize[3];
        this->m_interp.getPESize(pesize);
fabianw's avatar
fabianw committed
36
37
38
39
40
41
42
43
44
45
46
47
        gNx += (gNx % pesize[0]);
        gNy += (gNy % pesize[1]);
        gNz += (gNz % pesize[2]);
        const int Nx = gNx/pesize[0];
        const int Ny = gNy/pesize[1];
        const int Nz = gNz/pesize[2];
        assert(Nx > 0);
        assert(Ny > 0);
        assert(Nz > 0);
        const Real chx = this->m_interp.getExtentX() / gNx;
        const Real chy = this->m_interp.getExtentY() / gNy;
        const Real chz = this->m_interp.getExtentZ() / gNz;
48
49
50

        int peidx[3];
        this->m_interp.getPEIndex(peidx);
fabianw's avatar
fabianw committed
51
52
53
54
55
56
57
58
59
60
61
62
63
64
        const Real Ox = (this->m_interp.getExtentX()/pesize[0])*peidx[0];
        const Real Oy = (this->m_interp.getExtentY()/pesize[1])*peidx[1];
        const Real Oz = (this->m_interp.getExtentZ()/pesize[2])*peidx[2];
        const Real data_invh = 1.0 / this->m_interp.getH();
        Real origin[3];
        origin[0] = this->m_interp.getDataPos(ceil(Ox*data_invh)) - Ox;
        origin[1] = this->m_interp.getDataPos(ceil(Oy*data_invh)) - Oy;
        origin[2] = this->m_interp.getDataPos(ceil(Oz*data_invh)) - Oz;
        assert(origin[0] >= 0.0);
        assert(origin[1] >= 0.0);
        assert(origin[2] >= 0.0);
        this->m_interp.setOrigin(origin);

        this->_extractIso(isoval, Nx, Ny, Nz, chx, chy, chz, Ox, Oy, Oz, this->m_interp.isroot());
fabianw's avatar
fabianw committed
65
        MPI_Barrier(this->m_interp.getComm());
66
67
68
69
70
71
72
73
74

        _writePLY_MPI(fout, this->m_triList);
    }

private:

    void _writePLY_MPI(const std::string basename, const std::vector<Triangle>& tList)
    {
        const MPI_Comm comm = this->m_interp.getComm();
75
        unsigned long long mytriangles = tList.size();
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
        unsigned long long alltriangles;
        MPI_Reduce(&mytriangles, &alltriangles, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, comm);

        const bool isroot = this->m_interp.isroot();

        time_t rawtime;
        std::time(&rawtime);
        struct tm* timeinfo = std::localtime(&rawtime);
        char buf[256];
        std::strftime(buf, 256, "%A, %h %d %Y, %r %Z %z", timeinfo);

        const unsigned int one = 1;
        const unsigned char* const lsb = (const unsigned char*)(&one);
        const std::string endianess((*lsb) ? "binary_little_endian" : "binary_big_endian");

        // write header
        if (isroot)
        {
            std::ofstream ply(basename+".ply");
            ply << "ply" << std::endl;
            ply << "format ";
            if (this->m_parser("binary").asBool(true))
                ply << endianess;
            else
                ply << "ascii";
            ply << " 1.0" << std::endl;
            ply << "comment generated by IsoExtractor.h on " << buf << std::endl;
            ply << "element vertex " << alltriangles*3 << std::endl;
            ply << "property float x" << std::endl;
            ply << "property float y" << std::endl;
            ply << "property float z" << std::endl;
            ply << "element face " << alltriangles << std::endl;
            ply << "property list uchar int vertex_index" << std::endl;
            ply << "end_header" << std::endl;
            ply.close();
        }
        MPI_Barrier(comm);

        // TODO: (fabianw@mavt.ethz.ch; Tue Nov 29 14:04:45 2016) write proper
        // MPI I/O
        // write content
        // _writePLY_binary_MPI(comm, tList);

fabianw's avatar
fabianw committed
119
120
        // _writePLY_binary_MPI_lazy(tList, basename+".ply");
        _writePLY_binary_MPI(comm, tList, basename+".ply");
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    }

    void _writePLY_binary_MPI_lazy(const std::vector<Triangle>& tList, const std::string fname)
    {
        const MPI_Comm comm = this->m_interp.getComm();
        int size, rank;
        MPI_Comm_size(comm, &size);
        MPI_Comm_rank(comm, &rank);
        for (int r = 0; r < size; ++r)
        {
            if (r == rank)
            {
                std::ofstream ply(fname, std::ios::out | std::ios::app | std::ios::binary);
                for (auto t : tList)
                {
                    float tc[3][3];
                    for (int j = 0; j < 3; ++j)
                    {
                        tc[j][0] = t.vertices[j][0];
                        tc[j][1] = t.vertices[j][1];
                        tc[j][2] = t.vertices[j][2];
                    }
                    ply.write((char*)&tc[0][0], 9*sizeof(float));
                }
                ply.close();
            }
            MPI_Barrier(comm);
        }

        const char nVert = 3;
        std::vector<unsigned long long> allsize(size);
        unsigned long long mytri = tList.size();
        MPI_Allgather(&mytri, 1, MPI_UNSIGNED_LONG_LONG, allsize.data(), 1, MPI_UNSIGNED_LONG_LONG, comm);
        for (int r = 0; r < size; ++r)
        {
            if (r == rank)
            {
                unsigned long long start = 0;
                for (int i = 0; i < rank; ++i)
                    start += allsize[i];

                std::ofstream ply(fname, std::ios::out | std::ios::app | std::ios::binary);
                for (size_t i = 3*start; i < 3*start+tList.size()*3; i += 3)
                {
                    int tv[3] = {i, i+1, i+2};
                    ply.write(&nVert, sizeof(char));
                    ply.write((char*)&tv[0], 3*sizeof(int));
                }
                ply.close();
            }
            MPI_Barrier(comm);
        }
    }

fabianw's avatar
fabianw committed
175
176
177
178
179
180
181
182
    void _writePLY_binary_MPI(const MPI_Comm comm, const std::vector<Triangle>& tList, const std::string fname)
    {
        int rank, size;
        MPI_Status status;
        MPI_File file_id;
        MPI_Comm_rank(comm, &rank);
        MPI_Comm_size(comm, &size);

fabianw's avatar
fabianw committed
183
        int rc = MPI_File_open(MPI_COMM_SELF, (char*)fname.c_str(), MPI_MODE_APPEND | MPI_MODE_WRONLY, MPI_INFO_NULL, &file_id);
fabianw's avatar
fabianw committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
        if (rc)
        {
            if (0 == rank)
                std::cerr << "ERROR: IsoExtractorMPI: Can not open file." << std::endl;
            abort();
        }

        std::vector<unsigned long long> allsize(size);
        unsigned long long mytri = tList.size();
        MPI_Allgather(&mytri, 1, MPI_UNSIGNED_LONG_LONG, allsize.data(), 1, MPI_UNSIGNED_LONG_LONG, comm);

        long long offset = 0;
        MPI_File_get_position(file_id, &offset);

        assert(sizeof(Triangle) == 9*sizeof(float));
        for (int i = 0; i < rank; ++i)
            offset = offset + allsize[i]*sizeof(Triangle);

        Triangle* tarr = nullptr;
        if (allsize[rank] > 0)
        {
            tarr = new Triangle[allsize[rank]];
fabianw's avatar
fabianw committed
206
            char* const base = (char*)&tarr[0].vertices[0][0];
fabianw's avatar
fabianw committed
207
208
            for (size_t i=0; i<allsize[rank]; ++i)
            {
fabianw's avatar
fabianw committed
209
210
                char* const dst = base + i*sizeof(Triangle);
                const char* const src = (char*)&tList[i].vertices[0][0];
211
                memcpy(dst, src, sizeof(Triangle));
fabianw's avatar
fabianw committed
212
            }
213
214
215
216
217
218
219
220
221
            // for (size_t i=0; i<allsize[rank]; ++i)
            // {
            //     for (int j = 0; j < 3; ++j)
            //     {
            //         tarr[i].vertices[j][0] = tList[i].vertices[j][0];
            //         tarr[i].vertices[j][1] = tList[i].vertices[j][1];
            //         tarr[i].vertices[j][2] = tList[i].vertices[j][2];
            //     }
            // }
fabianw's avatar
fabianw committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258

            MPI_File_write_at(file_id, offset, (char *)&tarr[0].vertices[0][0], allsize[rank]*sizeof(Triangle), MPI_CHAR, &status);
        }
        offset = offset + allsize[rank]*sizeof(Triangle);
        MPI_Barrier(comm);
        MPI_Bcast(&offset, 1, MPI_LONG_LONG, size-1, comm);

        if (tarr != nullptr)
            delete [] tarr;

        for (int i = 0; i < rank; ++i)
            offset = offset + allsize[i]*sizeof(TriangleConnect);

        TriangleConnect* tcarr = nullptr;
        if (allsize[rank] > 0)
        {
            tcarr = new TriangleConnect[allsize[rank]];
            int base = 0;
            for (int i = 0; i < rank; ++i)
                base += 3*allsize[i];

            for (size_t i=0; i<allsize[rank]; ++i)
            {
                tcarr[i].n = 3;
                tcarr[i].v[0] = base + 3*i + 0;
                tcarr[i].v[1] = base + 3*i + 1;
                tcarr[i].v[2] = base + 3*i + 2;
            }

            MPI_File_write_at(file_id, offset, (char *)&tcarr[0].n, allsize[rank]*sizeof(TriangleConnect), MPI_CHAR, &status);
        }

        if (tcarr != nullptr)
            delete [] tcarr;

        MPI_File_close(&file_id);
    }
259
260
261
};

#endif /* ISOEXTRACTORMPI_H_XABTB0WN */