/[escript]/trunk/ripley/src/MultiBrick.cpp
ViewVC logotype

Annotation of /trunk/ripley/src/MultiBrick.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5540 - (hide annotations)
Tue Mar 17 06:55:35 2015 UTC (4 years, 2 months ago) by sshaw
File size: 68529 byte(s)
multi-resolution domains now support multi-process if the multi-dimension doesn't exceed 2 (MultiRectangle now supports MPI)
1 sshaw 5528
2     /*****************************************************************************
3     *
4     * Copyright (c) 2003-2015 by University of Queensland
5     * http://www.uq.edu.au
6     *
7     * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10     *
11     * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14     *
15     *****************************************************************************/
16    
17     #define ESNEEDPYTHON
18     #include "esysUtils/first.h"
19    
20     #include <ripley/MultiBrick.h>
21     #include <ripley/blocktools.h>
22     #include <ripley/domainhelpers.h>
23     #include <paso/SystemMatrix.h>
24     #include <esysUtils/esysFileWriter.h>
25     #include <esysUtils/EsysRandom.h>
26     #include <escript/DataFactory.h>
27     #include <escript/FunctionSpaceFactory.h>
28     #include <boost/scoped_array.hpp>
29 sshaw 5529 #include <boost/math/special_functions/fpclassify.hpp> // for isnan
30 sshaw 5528
31     #ifdef USE_NETCDF
32     #include <netcdfcpp.h>
33     #endif
34    
35     #if USE_SILO
36     #include <silo.h>
37     #ifdef ESYS_MPI
38     #include <pmpio.h>
39     #endif
40     #endif
41    
42     #define FIRST_QUAD 0.21132486540518711775
43     #define SECOND_QUAD 0.78867513459481288225
44    
45     #include <iomanip>
46     #include <limits>
47    
48    
49    
50     namespace bp = boost::python;
51     using esysUtils::FileWriter;
52     using escript::AbstractSystemMatrix;
53    
54 sshaw 5529 using std::vector;
55     using std::string;
56     using std::min;
57     using std::max;
58     using std::copy;
59     using std::ios;
60     using std::fill;
61    
62 sshaw 5528 namespace ripley {
63    
64     inline int indexOfMax(dim_t a, dim_t b, dim_t c)
65     {
66     if (a > b) {
67     if (c > a) {
68     return 2;
69     }
70     return 0;
71     } else if (b > c) {
72     return 1;
73     }
74     return 2;
75     }
76    
77     MultiBrick::MultiBrick(dim_t n0, dim_t n1, dim_t n2, double x0, double y0, double z0,
78     double x1, double y1, double z1, int d0, int d1, int d2,
79     const vector<double>& points, const vector<int>& tags,
80     const TagMap& tagnamestonums,
81     escript::SubWorld_ptr w,
82     unsigned int subdivisions
83     ) :
84     Brick(n0, n1, n2, x0, y0, z0, x1, y1, z1, d0, d1, d2, points, tags, tagnamestonums, w),
85     m_subdivisions(subdivisions)
86     {
87     if (m_mpiInfo->size != 1)
88     throw RipleyException("Multiresolution domains don't currently support multiple processes");
89    
90     if (subdivisions == 0 || (subdivisions & (subdivisions - 1)) != 0)
91     throw RipleyException("Element subdivisions must be a power of two");
92    
93     if (d0 == 0 || d1 == 0)
94     throw RipleyException("Domain subdivisions must be positive");
95    
96     dim_t oldNN[3] = {0};
97    
98     for (int i = 0; i < 3; i++) {
99     m_NE[i] *= subdivisions;
100     oldNN[i] = m_NN[i];
101     m_NN[i] = m_NE[i] + 1;
102     m_gNE[i] *= subdivisions;
103     m_ownNE[i] *= subdivisions;
104     m_dx[i] /= subdivisions;
105     m_faceCount[i] *= subdivisions;
106     m_faceCount[2+i] *= subdivisions;
107     }
108    
109     // bottom-left node is at (offset0,offset1) in global mesh
110     m_offset[0] = m_gNE[0]*subdivisions/d0*(m_mpiInfo->rank%d0);
111     m_offset[1] = m_gNE[1]*subdivisions/d1*(m_mpiInfo->rank/d0);
112     m_offset[2] = m_gNE[2]*subdivisions/d2*(m_mpiInfo->rank/(d0*d1));
113     populateSampleIds();
114    
115     const dim_t nDirac = m_diracPoints.size();
116     #pragma omp parallel for
117     for (int i = 0; i < nDirac; i++) {
118     const dim_t node = m_diracPoints[i].node;
119     const dim_t x = node % oldNN[0];
120     const dim_t y = (node % (oldNN[0]*oldNN[1])) / oldNN[0];
121     const dim_t z = node / (oldNN[0]*oldNN[1]);
122     m_diracPoints[i].node = INDEX3(x*subdivisions, y*subdivisions, z*subdivisions,
123     m_NN[0], m_NN[1]);
124     m_diracPointNodeIDs[i] = m_diracPoints[i].node;
125     }
126     }
127    
128     MultiBrick::~MultiBrick()
129     {
130     }
131    
132    
133     void MultiBrick::validateInterpolationAcross(int fsType_source,
134     const escript::AbstractDomain& domain, int fsType_target) const
135     {
136     const MultiBrick *other = dynamic_cast<const MultiBrick *>(&domain);
137     if (other == NULL)
138     throw RipleyException("Invalid interpolation: domains must both be instances of MultiBrick");
139    
140     const double *len = other->getLength();
141     const int *subdivs = other->getNumSubdivisionsPerDim();
142     const int *elements = other->getNumElementsPerDim();
143     const unsigned int level = other->getNumSubdivisionsPerElement();
144     const unsigned int factor = m_subdivisions > level ? m_subdivisions/level : level/m_subdivisions;
145     if ((factor & (factor - 1)) != 0) //factor == 2**x
146     throw RipleyException("Invalid interpolation: elemental subdivisions of each domain must be powers of two");
147    
148     if (other->getMPIComm() != m_mpiInfo->comm)
149     throw RipleyException("Invalid interpolation: Domains are on different communicators");
150    
151     for (int i = 0; i < m_numDim; i++) {
152     if (m_length[i] != len[i])
153     throw RipleyException("Invalid interpolation: domain length mismatch");
154     if (m_NX[i] != subdivs[i])
155     throw RipleyException("Invalid interpolation: domain process subdivision mismatch");
156     if (m_subdivisions > level) {
157     if (m_ownNE[i]/elements[i] != factor)
158     throw RipleyException("Invalid interpolation: element factor mismatch");
159     } else {
160     if (elements[i]/m_ownNE[i] != factor)
161     throw RipleyException("Invalid interpolation: element factor mismatch");
162     }
163     }
164     }
165    
166     void MultiBrick::interpolateNodesToNodesFiner(const escript::Data& source,
167     escript::Data& target, const MultiBrick& other) const
168     {
169     const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
170     const dim_t NN0 = m_NN[0], NN1 = m_NN[1], NN2 = m_NN[2], *otherNN = other.getNumNodesPerDim();
171     const dim_t numComp = source.getDataPointSize();
172     #pragma omp parallel for
173     for (dim_t nz = 0; nz < NN2 - 1; nz++) { //source nodes
174     for (dim_t ny = 0; ny < NN1 - 1; ny++) {
175     for (dim_t nx = 0; nx < NN0 - 1; nx++) {
176     const double *x0y0z0 = source.getSampleDataRO(INDEX3(nx, ny, nz, NN0, NN1));
177     const double *x0y1z0 = source.getSampleDataRO(INDEX3(nx, ny+1, nz, NN0, NN1));
178     const double *x1y0z0 = source.getSampleDataRO(INDEX3(nx+1, ny, nz, NN0, NN1));
179     const double *x1y1z0 = source.getSampleDataRO(INDEX3(nx+1, ny+1, nz, NN0, NN1));
180     const double *x0y0z1 = source.getSampleDataRO(INDEX3(nx, ny, nz+1, NN0, NN1));
181     const double *x0y1z1 = source.getSampleDataRO(INDEX3(nx, ny+1, nz+1, NN0, NN1));
182     const double *x1y0z1 = source.getSampleDataRO(INDEX3(nx+1, ny, nz+1, NN0, NN1));
183     const double *x1y1z1 = source.getSampleDataRO(INDEX3(nx+1, ny+1, nz+1, NN0, NN1));
184     const double origin[3] = {getLocalCoordinate(nx, 0),
185     getLocalCoordinate(ny, 1),
186     getLocalCoordinate(nz, 2)};
187     for (int sz = 0; sz < scaling + 1; sz++) { //target nodes
188     const double z = (other.getLocalCoordinate(nz*scaling+sz, 2) - origin[2]) / m_dx[2];
189     for (int sy = 0; sy < scaling + 1; sy++) {
190     const double y = (other.getLocalCoordinate(ny*scaling+sy, 1) - origin[1]) / m_dx[1];
191     for (int sx = 0; sx < scaling + 1; sx++) {
192     const double x = (other.getLocalCoordinate(nx*scaling+sx, 0) - origin[0]) / m_dx[0];
193     double *out = target.getSampleDataRW(INDEX3(nx*scaling+sx, ny*scaling+sy, nz*scaling+sz, otherNN[0], otherNN[1]));
194     for (int comp = 0; comp < numComp; comp++) {
195     out[comp] = x0y0z0[comp]*(1-x)*(1-y)*(1-z) + x1y0z0[comp]*x*(1-y)*(1-z) + x0y1z0[comp]*(1-x)*y*(1-z) + x1y1z0[comp]*x*y*(1-z)
196     + x0y0z1[comp]*(1-x)*(1-y)*z + x1y0z1[comp]*x*(1-y)*z + x0y1z1[comp]*(1-x)*y*z + x1y1z1[comp]*x*y*z;
197     }
198     }
199     }
200     }
201     }
202     }
203     }
204     }
205    
206     void MultiBrick::interpolateReducedToElementsFiner(const escript::Data& source,
207     escript::Data& target, const MultiBrick& other) const
208     {
209     const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
210     const dim_t numComp = source.getDataPointSize();
211     //for each of ours
212     #pragma omp parallel for
213     for (dim_t ez = 0; ez < m_NE[2]; ez++) {
214     for (dim_t ey = 0; ey < m_NE[1]; ey++) {
215     for (dim_t ex = 0; ex < m_NE[0]; ex++) {
216     const double *in = source.getSampleDataRO(INDEX3(ex, ey, ez, m_NE[0], m_NE[1]));
217     //for each subelement
218     for (dim_t sz = 0; sz < scaling; sz++) {
219     const dim_t tz = ez*scaling + sz;
220     for (dim_t sy = 0; sy < scaling; sy++) {
221     const dim_t ty = ey*scaling + sy;
222     for (dim_t sx = 0; sx < scaling; sx++) {
223     const dim_t tx = ex*scaling + sx;
224     double *out = target.getSampleDataRW(INDEX3(tx, ty, tz, m_NE[0]*scaling, m_NE[1]*scaling));
225     for (dim_t comp = 0; comp < numComp; comp++) {
226     const double quadvalue = in[comp];
227     for (int i = 0; i < 8; i++) {
228     out[comp + i*numComp] = quadvalue;
229     }
230     }
231     }
232     }
233     }
234     }
235     }
236     }
237     }
238    
239     void MultiBrick::interpolateReducedToReducedFiner(const escript::Data& source,
240     escript::Data& target, const MultiBrick& other) const
241     {
242     const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
243     const dim_t numComp = source.getDataPointSize();
244     //for each of ours
245     #pragma omp parallel for
246     for (dim_t ey = 0; ey < m_NE[1]; ey++) {
247     for (dim_t ex = 0; ex < m_NE[0]; ex++) {
248     const double *in = source.getSampleDataRO(ex + ey*m_NE[0]);
249     //for each subelement
250     for (dim_t sy = 0; sy < scaling; sy++) {
251     const dim_t ty = ey*scaling + sy;
252     for (dim_t sx = 0; sx < scaling; sx++) {
253     const dim_t tx = ex*scaling + sx;
254     double *out = target.getSampleDataRW(tx + ty*m_NE[0]*scaling);
255     for (dim_t comp = 0; comp < numComp; comp++) {
256     out[comp] = in[comp];
257     }
258     }
259     }
260     }
261     }
262     }
263    
264     void MultiBrick::interpolateNodesToElementsFiner(const escript::Data& source,
265     escript::Data& target, const MultiBrick& other) const
266     {
267     const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
268     const dim_t NE0 = m_NE[0], NE1 = m_NE[1], NE2 = m_NE[2], *theirNE = other.getNumElementsPerDim();
269     const dim_t numComp = source.getDataPointSize();
270     #pragma omp parallel for
271     for (dim_t ez = 0; ez < NE2; ez++) { //source nodes
272     for (dim_t ey = 0; ey < NE1; ey++) {
273     for (dim_t ex = 0; ex < NE0; ex++) {
274     const double *points[8] = {
275     source.getSampleDataRO(INDEX3(ex, ey, ez, NE0+1, NE1+1)),
276     source.getSampleDataRO(INDEX3(ex+1, ey, ez, NE0+1, NE1+1)),
277     source.getSampleDataRO(INDEX3(ex, ey+1, ez, NE0+1, NE1+1)),
278     source.getSampleDataRO(INDEX3(ex+1, ey+1, ez, NE0+1, NE1+1)),
279     source.getSampleDataRO(INDEX3(ex, ey, ez+1, NE0+1, NE1+1)),
280     source.getSampleDataRO(INDEX3(ex+1, ey, ez+1, NE0+1, NE1+1)),
281     source.getSampleDataRO(INDEX3(ex, ey+1, ez+1, NE0+1, NE1+1)),
282     source.getSampleDataRO(INDEX3(ex+1, ey+1, ez+1, NE0+1, NE1+1)),
283     };
284     const double origin[3] = {getLocalCoordinate(ex, 0),
285     getLocalCoordinate(ey, 1),
286     getLocalCoordinate(ez, 2)
287     };
288     for (int sz = 0; sz < scaling; sz++) { //target elements
289     for (int sy = 0; sy < scaling; sy++) {
290     for (int sx = 0; sx < scaling; sx++) {
291     const double x1 = (other.getLocalCoordinate(ex*scaling+sx, 0) - origin[0]) / m_dx[0] + FIRST_QUAD/scaling;
292     const double x2 = x1 + (SECOND_QUAD - FIRST_QUAD)/scaling;
293     const double y1 = (other.getLocalCoordinate(ey*scaling+sy, 1) - origin[1]) / m_dx[1] + FIRST_QUAD/scaling;
294     const double y2 = y1 + (SECOND_QUAD - FIRST_QUAD)/scaling;
295     const double z1 = (other.getLocalCoordinate(ez*scaling+sz, 2) - origin[2]) / m_dx[2] + FIRST_QUAD/scaling;
296     const double z2 = z1 + (SECOND_QUAD - FIRST_QUAD)/scaling;
297     double *out = target.getSampleDataRW(INDEX3(ex*scaling+sx, ey*scaling+sy, ez*scaling+sz, theirNE[0], theirNE[1]));
298     for (int comp = 0; comp < numComp; comp++) {
299     out[INDEX4(comp, 0, 0, 0, numComp, 2, 2)] =
300     points[0][comp]*(1-x1)*(1-y1)*(1-z1)
301     + points[1][comp]*x1*(1-y1)*(1-z1)
302     + points[2][comp]*(1-x1)*y1*(1-z1)
303     + points[3][comp]*x1*y1*(1-z1)
304     + points[4][comp]*(1-x1)*(1-y1)*z1
305     + points[5][comp]*x1*(1-y1)*z1
306     + points[6][comp]*(1-x1)*y1*z1
307     + points[7][comp]*x1*y1*z1;
308     out[INDEX4(comp, 1, 0, 0, numComp, 2, 2)] =
309     points[0][comp]*(1-x2)*(1-y1)*(1-z1)
310     + points[1][comp]*x2*(1-y1)*(1-z1)
311     + points[2][comp]*(1-x2)*y1*(1-z1)
312     + points[3][comp]*x2*y1*(1-z1)
313     + points[4][comp]*(1-x2)*(1-y1)*z1
314     + points[5][comp]*x2*(1-y1)*z1
315     + points[6][comp]*(1-x2)*y1*z1
316     + points[7][comp]*x2*y1*z1;
317     out[INDEX4(comp, 0, 1, 0, numComp, 2, 2)] =
318     points[0][comp]*(1-x1)*(1-y2)*(1-z1)
319     + points[1][comp]*x1*(1-y2)*(1-z1)
320     + points[2][comp]*(1-x1)*y2*(1-z1)
321     + points[3][comp]*x1*y2*(1-z1)
322     + points[4][comp]*(1-x1)*(1-y2)*z1
323     + points[5][comp]*x1*(1-y2)*z1
324     + points[6][comp]*(1-x1)*y2*z1
325     + points[7][comp]*x1*y2*z1;
326     out[INDEX4(comp, 1, 1, 0, numComp, 2, 2)] =
327     points[0][comp]*(1-x2)*(1-y2)*(1-z1)
328     + points[1][comp]*x2*(1-y2)*(1-z1)
329     + points[2][comp]*(1-x2)*y2*(1-z1)
330     + points[3][comp]*x2*y2*(1-z1)
331     + points[4][comp]*(1-x2)*(1-y2)*z1
332     + points[5][comp]*x2*(1-y2)*z1
333     + points[6][comp]*(1-x2)*y2*z1
334     + points[7][comp]*x2*y2*z1;
335     out[INDEX4(comp, 0, 0, 1, numComp, 2, 2)] =
336     points[0][comp]*(1-x1)*(1-y1)*(1-z2)
337     + points[1][comp]*x1*(1-y1)*(1-z2)
338     + points[2][comp]*(1-x1)*y1*(1-z2)
339     + points[3][comp]*x1*y1*(1-z2)
340     + points[4][comp]*(1-x1)*(1-y1)*z2
341     + points[5][comp]*x1*(1-y1)*z2
342     + points[6][comp]*(1-x1)*y1*z2
343     + points[7][comp]*x1*y1*z2;
344     out[INDEX4(comp, 1, 0, 1, numComp, 2, 2)] =
345     points[0][comp]*(1-x2)*(1-y1)*(1-z2)
346     + points[1][comp]*x2*(1-y1)*(1-z2)
347     + points[2][comp]*(1-x2)*y1*(1-z2)
348     + points[3][comp]*x2*y1*(1-z2)
349     + points[4][comp]*(1-x2)*(1-y1)*z2
350     + points[5][comp]*x2*(1-y1)*z2
351     + points[6][comp]*(1-x2)*y1*z2
352     + points[7][comp]*x2*y1*z2;
353     out[INDEX4(comp, 0, 1, 1, numComp, 2, 2)] =
354     points[0][comp]*(1-x1)*(1-y2)*(1-z2)
355     + points[1][comp]*x1*(1-y2)*(1-z2)
356     + points[2][comp]*(1-x1)*y2*(1-z2)
357     + points[3][comp]*x1*y2*(1-z2)
358     + points[4][comp]*(1-x1)*(1-y2)*z2
359     + points[5][comp]*x1*(1-y2)*z2
360     + points[6][comp]*(1-x1)*y2*z2
361     + points[7][comp]*x1*y2*z2;
362     out[INDEX4(comp, 1, 1, 1, numComp, 2, 2)] =
363     points[0][comp]*(1-x2)*(1-y2)*(1-z2)
364     + points[1][comp]*x2*(1-y2)*(1-z2)
365     + points[2][comp]*(1-x2)*y2*(1-z2)
366     + points[3][comp]*x2*y2*(1-z2)
367     + points[4][comp]*(1-x2)*(1-y2)*z2
368     + points[5][comp]*x2*(1-y2)*z2
369     + points[6][comp]*(1-x2)*y2*z2
370     + points[7][comp]*x2*y2*z2;
371     }
372     }
373     }
374     }
375     }
376     }
377     }
378     }
379    
380     void MultiBrick::interpolateElementsToElementsCoarser(const escript::Data& source,
381     escript::Data& target, const MultiBrick& other) const
382     {
383     const int scaling = m_subdivisions/other.getNumSubdivisionsPerElement();
384     const double scaling_volume = (1./scaling)*(1./scaling)*(1./scaling);
385     const dim_t *theirNE = other.getNumElementsPerDim();
386     const dim_t numComp = source.getDataPointSize();
387    
388     std::vector<double> points(scaling*2, 0);
389     std::vector<double> first_lagrange(scaling*2, 1);
390     std::vector<double> second_lagrange(scaling*2, 1);
391    
392     for (int i = 0; i < scaling*2; i+=2) {
393     points[i] = (i/2 + FIRST_QUAD)/scaling;
394     points[i+1] = (i/2 + SECOND_QUAD)/scaling;
395     }
396    
397     for (int i = 0; i < scaling*2; i++) {
398     first_lagrange[i] = (points[i] - SECOND_QUAD) / (FIRST_QUAD - SECOND_QUAD);
399     second_lagrange[i] = (points[i] - FIRST_QUAD) / (SECOND_QUAD - FIRST_QUAD);
400     }
401    
402     //for each of theirs
403     #pragma omp parallel for
404     for (dim_t tz = 0; tz < theirNE[2]; tz++) {
405     for (dim_t ty = 0; ty < theirNE[1]; ty++) {
406     for (dim_t tx = 0; tx < theirNE[0]; tx++) {
407     double *out = target.getSampleDataRW(INDEX3(tx, ty, tz, theirNE[0], theirNE[1]));
408     //for each subelement
409     for (dim_t sz = 0; sz < scaling; sz++) {
410     const dim_t ez = tz*scaling + sz;
411     for (dim_t sy = 0; sy < scaling; sy++) {
412     const dim_t ey = ty*scaling + sy;
413     for (dim_t sx = 0; sx < scaling; sx++) {
414     const dim_t ex = tx*scaling + sx;
415     const double *in = source.getSampleDataRO(INDEX3(ex, ey, ez, m_NE[0], m_NE[1]));
416     for (int quad = 0; quad < 8; quad++) {
417     int lx = sx*2 + quad%2;
418     int ly = sy*2 + (quad%4)/2;
419     int lz = sz*2 + quad/4;
420     for (dim_t comp = 0; comp < numComp; comp++) {
421     const double quadvalue = scaling_volume * in[comp + quad*numComp];
422     out[INDEX4(comp, 0, 0, 0, numComp, 2, 2)] += quadvalue * first_lagrange[lx] * first_lagrange[ly] * first_lagrange[lz];
423     out[INDEX4(comp, 1, 0, 0, numComp, 2, 2)] += quadvalue * second_lagrange[lx] * first_lagrange[ly] * first_lagrange[lz];
424     out[INDEX4(comp, 0, 1, 0, numComp, 2, 2)] += quadvalue * first_lagrange[lx] * second_lagrange[ly] * first_lagrange[lz];
425     out[INDEX4(comp, 1, 1, 0, numComp, 2, 2)] += quadvalue * second_lagrange[lx] * second_lagrange[ly] * first_lagrange[lz];
426     out[INDEX4(comp, 0, 0, 1, numComp, 2, 2)] += quadvalue * first_lagrange[lx] * first_lagrange[ly] * second_lagrange[lz];
427     out[INDEX4(comp, 1, 0, 1, numComp, 2, 2)] += quadvalue * second_lagrange[lx] * first_lagrange[ly] * second_lagrange[lz];
428     out[INDEX4(comp, 0, 1, 1, numComp, 2, 2)] += quadvalue * first_lagrange[lx] * second_lagrange[ly] * second_lagrange[lz];
429     out[INDEX4(comp, 1, 1, 1, numComp, 2, 2)] += quadvalue * second_lagrange[lx] * second_lagrange[ly] * second_lagrange[lz];
430     }
431     }
432     }
433     }
434     }
435     }
436     }
437     }
438     }
439    
440    
441     void MultiBrick::interpolateElementsToElementsFiner(const escript::Data& source,
442     escript::Data& target, const MultiBrick& other) const
443     {
444     const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
445     const dim_t numComp = source.getDataPointSize();
446    
447     std::vector<double> points(scaling*2, 0);
448     std::vector<double> lagranges(scaling*4, 1);
449    
450     for (int i = 0; i < scaling*2; i+=2) {
451     points[i] = (i/2 + FIRST_QUAD)/scaling;
452     points[i+1] = (i/2 + SECOND_QUAD)/scaling;
453     }
454     for (int i = 0; i < scaling*2; i++) {
455     lagranges[i] = (points[i] - SECOND_QUAD) / (FIRST_QUAD - SECOND_QUAD);
456     lagranges[i + 2*scaling] = (points[i] - FIRST_QUAD) / (SECOND_QUAD - FIRST_QUAD);
457     }
458     //for each of ours
459     #pragma omp parallel for
460     for (dim_t ez = 0; ez < m_NE[2]; ez++) {
461     for (dim_t ey = 0; ey < m_NE[1]; ey++) {
462     for (dim_t ex = 0; ex < m_NE[0]; ex++) {
463     const double *in = source.getSampleDataRO(INDEX3(ex, ey, ez, m_NE[0], m_NE[1]));
464     //for each subelement
465     for (dim_t sz = 0; sz < scaling; sz++) {
466     const dim_t tz = ez*scaling + sz;
467     for (dim_t sy = 0; sy < scaling; sy++) {
468     const dim_t ty = ey*scaling + sy;
469     for (dim_t sx = 0; sx < scaling; sx++) {
470     const dim_t tx = ex*scaling + sx;
471     double *out = target.getSampleDataRW(INDEX3(tx, ty, tz, m_NE[0]*scaling, m_NE[1]*scaling));
472     for (int quad = 0; quad < 8; quad++) {
473     const int lx = scaling*2*(quad%2) + sx*2;
474     const int ly = scaling*2*((quad%4)/2) + sy*2;
475     const int lz = scaling*2*(quad/4) + sz*2;
476     for (dim_t comp = 0; comp < numComp; comp++) {
477     const double quadvalue = in[comp + quad*numComp];
478     out[INDEX4(comp, 0, 0, 0, numComp, 2, 2)] += quadvalue * lagranges[lx] * lagranges[ly] * lagranges[lz];
479     out[INDEX4(comp, 0, 1, 0, numComp, 2, 2)] += quadvalue * lagranges[lx] * lagranges[ly+1] * lagranges[lz];
480     out[INDEX4(comp, 1, 0, 0, numComp, 2, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly] * lagranges[lz];
481     out[INDEX4(comp, 1, 1, 0, numComp, 2, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly+1] * lagranges[lz];
482     out[INDEX4(comp, 0, 0, 1, numComp, 2, 2)] += quadvalue * lagranges[lx] * lagranges[ly] * lagranges[lz+1];
483     out[INDEX4(comp, 0, 1, 1, numComp, 2, 2)] += quadvalue * lagranges[lx] * lagranges[ly+1] * lagranges[lz+1];
484     out[INDEX4(comp, 1, 0, 1, numComp, 2, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly] * lagranges[lz+1];
485     out[INDEX4(comp, 1, 1, 1, numComp, 2, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly+1] * lagranges[lz+1];
486     }
487     }
488     }
489     }
490     }
491     }
492     }
493     }
494     }
495    
496     void MultiBrick::interpolateAcross(escript::Data& target,
497     const escript::Data& source) const
498     {
499     const MultiBrick *other =
500     dynamic_cast<const MultiBrick *>(target.getDomain().get());
501     if (other == NULL)
502     throw RipleyException("Invalid interpolation: Domains must both be instances of MultiBrick");
503     //shouldn't ever happen, but I want to know if it does
504     if (other == this)
505     throw RipleyException("interpolateAcross: this domain is the target");
506    
507     validateInterpolationAcross(source.getFunctionSpace().getTypeCode(),
508     *(target.getDomain().get()), target.getFunctionSpace().getTypeCode());
509     int fsSource = source.getFunctionSpace().getTypeCode();
510     int fsTarget = target.getFunctionSpace().getTypeCode();
511    
512     std::stringstream msg;
513     msg << "Invalid interpolation: interpolation not implemented for function space "
514     << functionSpaceTypeAsString(fsSource)
515     << " -> "
516     << functionSpaceTypeAsString(fsTarget);
517     if (other->getNumSubdivisionsPerElement() > getNumSubdivisionsPerElement()) {
518     switch (fsSource) {
519     case Nodes:
520     switch (fsTarget) {
521     case Nodes:
522     case ReducedNodes:
523     case DegreesOfFreedom:
524     case ReducedDegreesOfFreedom:
525     interpolateNodesToNodesFiner(source, target, *other);
526     return;
527     case Elements:
528     interpolateNodesToElementsFiner(source, target, *other);
529     return;
530     }
531     break;
532     case Elements:
533     switch (fsTarget) {
534     case Elements:
535     interpolateElementsToElementsFiner(source, target, *other);
536     return;
537     }
538     break;
539     case ReducedElements:
540     switch (fsTarget) {
541     case Elements:
542     interpolateReducedToElementsFiner(source, target, *other);
543     return;
544     }
545     break;
546     }
547     msg << " when target is a finer mesh";
548     } else {
549     switch (fsSource) {
550     case Nodes:
551     switch (fsTarget) {
552     case Elements:
553     escript::Data elements=escript::Vector(0., escript::function(*this), true);
554     interpolateNodesOnElements(elements, source, false);
555     interpolateElementsToElementsCoarser(elements, target, *other);
556     return;
557     }
558     break;
559     case Elements:
560     switch (fsTarget) {
561     case Elements:
562     interpolateElementsToElementsCoarser(source, target, *other);
563     return;
564     }
565     break;
566     }
567     msg << " when target is a coarser mesh";
568     }
569     throw RipleyException(msg.str());
570     }
571    
572     string MultiBrick::getDescription() const
573     {
574     return "ripley::MultiBrick";
575     }
576    
577     bool MultiBrick::operator==(const AbstractDomain& other) const
578     {
579     const MultiBrick* o=dynamic_cast<const MultiBrick*>(&other);
580     if (o) {
581     return (RipleyDomain::operator==(other) &&
582     m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1] && m_gNE[2]==o->m_gNE[2]
583     && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1] && m_origin[2]==o->m_origin[2]
584     && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1] && m_length[2]==o->m_length[2]
585     && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1] && m_NX[2]==o->m_NX[2]
586     && m_subdivisions == o->m_subdivisions);
587     }
588    
589     return false;
590     }
591    
592     void MultiBrick::readNcGrid(escript::Data& out, string filename, string varname,
593     const ReaderParameters& params) const
594     {
595     if (m_subdivisions != 1)
596     throw RipleyException("Non-parent MultiBricks cannot read datafiles");
597     Brick::readNcGrid(out, filename, varname, params);
598     }
599    
600     void MultiBrick::readBinaryGrid(escript::Data& out, string filename,
601     const ReaderParameters& params) const
602     {
603     if (m_subdivisions != 1)
604     throw RipleyException("Non-parent MultiBricks cannot read datafiles");
605     Brick::readBinaryGrid(out, filename, params);
606     }
607    
608     #ifdef USE_BOOSTIO
609     void MultiBrick::readBinaryGridFromZipped(escript::Data& out, string filename,
610     const ReaderParameters& params) const
611     {
612     if (m_subdivisions != 1)
613     throw RipleyException("Non-parent MultiBricks cannot read datafiles");
614     Brick::readBinaryGridFromZipped(out, filename, params);
615     }
616     #endif
617    
618     void MultiBrick::writeBinaryGrid(const escript::Data& in, string filename,
619     int byteOrder, int dataType) const
620     {
621     if (m_subdivisions != 1)
622     throw RipleyException("Non-parent MultiBricks cannot read datafiles");
623     Brick::writeBinaryGrid(in, filename, byteOrder, dataType);
624     }
625    
626     void MultiBrick::dump(const string& fileName) const
627     {
628     if (m_subdivisions != 1)
629     throw RipleyException("Non-parent MultiBricks dump not implemented");
630     Brick::dump(fileName);
631     }
632    
633     const dim_t* MultiBrick::borrowSampleReferenceIDs(int fsType) const
634     {
635     switch (fsType) {
636     case Nodes:
637     case ReducedNodes: //FIXME: reduced
638     return &m_nodeId[0];
639     case DegreesOfFreedom:
640     case ReducedDegreesOfFreedom: //FIXME: reduced
641     return &m_dofId[0];
642     case Elements:
643     case ReducedElements:
644     return &m_elementId[0];
645     case FaceElements:
646     case ReducedFaceElements:
647     return &m_faceId[0];
648     case Points:
649     return &m_diracPointNodeIDs[0];
650     default:
651     break;
652     }
653    
654 sshaw 5529 std::stringstream msg;
655 sshaw 5528 msg << "borrowSampleReferenceIDs: invalid function space type "<<fsType;
656     throw RipleyException(msg.str());
657     }
658    
659     bool MultiBrick::ownSample(int fsType, index_t id) const
660     {
661     if (getMPISize()==1)
662     return true;
663    
664     switch (fsType) {
665     case Nodes:
666     case ReducedNodes: //FIXME: reduced
667     return (m_dofMap[id] < getNumDOF());
668     case DegreesOfFreedom:
669     case ReducedDegreesOfFreedom:
670     return true;
671     case Elements:
672     case ReducedElements:
673     {
674     // check ownership of element's _last_ node
675     const index_t x=id%m_NE[0] + 1;
676     const index_t y=id%(m_NE[0]*m_NE[1])/m_NE[0] + 1;
677     const index_t z=id/(m_NE[0]*m_NE[1]) + 1;
678     return (m_dofMap[x + m_NN[0]*y + m_NN[0]*m_NN[1]*z] < getNumDOF());
679     }
680     case FaceElements:
681     case ReducedFaceElements:
682     {
683     // check ownership of face element's last node
684     dim_t n=0;
685     for (size_t i=0; i<6; i++) {
686     n+=m_faceCount[i];
687     if (id<n) {
688     const index_t j=id-n+m_faceCount[i];
689     if (i>=4) { // front or back
690     const index_t first=(i==4 ? 0 : m_NN[0]*m_NN[1]*(m_NN[2]-1));
691     return (m_dofMap[first+j%m_NE[0]+1+(j/m_NE[0]+1)*m_NN[0]] < getNumDOF());
692     } else if (i>=2) { // bottom or top
693     const index_t first=(i==2 ? 0 : m_NN[0]*(m_NN[1]-1));
694     return (m_dofMap[first+j%m_NE[0]+1+(j/m_NE[0]+1)*m_NN[0]*m_NN[1]] < getNumDOF());
695     } else { // left or right
696     const index_t first=(i==0 ? 0 : m_NN[0]-1);
697     return (m_dofMap[first+(j%m_NE[1]+1)*m_NN[0]+(j/m_NE[1]+1)*m_NN[0]*m_NN[1]] < getNumDOF());
698     }
699     }
700     }
701     return false;
702     }
703     default:
704     break;
705     }
706    
707 sshaw 5529 std::stringstream msg;
708 sshaw 5528 msg << "ownSample: invalid function space type " << fsType;
709     throw RipleyException(msg.str());
710     }
711    
712     void MultiBrick::setToNormal(escript::Data& out) const
713     {
714     const dim_t NE0 = m_NE[0];
715     const dim_t NE1 = m_NE[1];
716     const dim_t NE2 = m_NE[2];
717    
718     if (out.getFunctionSpace().getTypeCode() == FaceElements) {
719     out.requireWrite();
720     #pragma omp parallel
721     {
722     if (m_faceOffset[0] > -1) {
723     #pragma omp for nowait
724     for (index_t k2 = 0; k2 < NE2; ++k2) {
725     for (index_t k1 = 0; k1 < NE1; ++k1) {
726     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
727     // set vector at four quadrature points
728     *o++ = -1.; *o++ = 0.; *o++ = 0.;
729     *o++ = -1.; *o++ = 0.; *o++ = 0.;
730     *o++ = -1.; *o++ = 0.; *o++ = 0.;
731     *o++ = -1.; *o++ = 0.; *o = 0.;
732     }
733     }
734     }
735    
736     if (m_faceOffset[1] > -1) {
737     #pragma omp for nowait
738     for (index_t k2 = 0; k2 < NE2; ++k2) {
739     for (index_t k1 = 0; k1 < NE1; ++k1) {
740     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
741     // set vector at four quadrature points
742     *o++ = 1.; *o++ = 0.; *o++ = 0.;
743     *o++ = 1.; *o++ = 0.; *o++ = 0.;
744     *o++ = 1.; *o++ = 0.; *o++ = 0.;
745     *o++ = 1.; *o++ = 0.; *o = 0.;
746     }
747     }
748     }
749    
750     if (m_faceOffset[2] > -1) {
751     #pragma omp for nowait
752     for (index_t k2 = 0; k2 < NE2; ++k2) {
753     for (index_t k0 = 0; k0 < NE0; ++k0) {
754     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
755     // set vector at four quadrature points
756     *o++ = 0.; *o++ = -1.; *o++ = 0.;
757     *o++ = 0.; *o++ = -1.; *o++ = 0.;
758     *o++ = 0.; *o++ = -1.; *o++ = 0.;
759     *o++ = 0.; *o++ = -1.; *o = 0.;
760     }
761     }
762     }
763    
764     if (m_faceOffset[3] > -1) {
765     #pragma omp for nowait
766     for (index_t k2 = 0; k2 < NE2; ++k2) {
767     for (index_t k0 = 0; k0 < NE0; ++k0) {
768     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
769     // set vector at four quadrature points
770     *o++ = 0.; *o++ = 1.; *o++ = 0.;
771     *o++ = 0.; *o++ = 1.; *o++ = 0.;
772     *o++ = 0.; *o++ = 1.; *o++ = 0.;
773     *o++ = 0.; *o++ = 1.; *o = 0.;
774     }
775     }
776     }
777    
778     if (m_faceOffset[4] > -1) {
779     #pragma omp for nowait
780     for (index_t k1 = 0; k1 < NE1; ++k1) {
781     for (index_t k0 = 0; k0 < NE0; ++k0) {
782     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
783     // set vector at four quadrature points
784     *o++ = 0.; *o++ = 0.; *o++ = -1.;
785     *o++ = 0.; *o++ = 0.; *o++ = -1.;
786     *o++ = 0.; *o++ = 0.; *o++ = -1.;
787     *o++ = 0.; *o++ = 0.; *o = -1.;
788     }
789     }
790     }
791    
792     if (m_faceOffset[5] > -1) {
793     #pragma omp for nowait
794     for (index_t k1 = 0; k1 < NE1; ++k1) {
795     for (index_t k0 = 0; k0 < NE0; ++k0) {
796     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
797     // set vector at four quadrature points
798     *o++ = 0.; *o++ = 0.; *o++ = 1.;
799     *o++ = 0.; *o++ = 0.; *o++ = 1.;
800     *o++ = 0.; *o++ = 0.; *o++ = 1.;
801     *o++ = 0.; *o++ = 0.; *o = 1.;
802     }
803     }
804     }
805     } // end of parallel section
806     } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
807     out.requireWrite();
808     #pragma omp parallel
809     {
810     if (m_faceOffset[0] > -1) {
811     #pragma omp for nowait
812     for (index_t k2 = 0; k2 < NE2; ++k2) {
813     for (index_t k1 = 0; k1 < NE1; ++k1) {
814     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
815     *o++ = -1.;
816     *o++ = 0.;
817     *o = 0.;
818     }
819     }
820     }
821    
822     if (m_faceOffset[1] > -1) {
823     #pragma omp for nowait
824     for (index_t k2 = 0; k2 < NE2; ++k2) {
825     for (index_t k1 = 0; k1 < NE1; ++k1) {
826     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
827     *o++ = 1.;
828     *o++ = 0.;
829     *o = 0.;
830     }
831     }
832     }
833    
834     if (m_faceOffset[2] > -1) {
835     #pragma omp for nowait
836     for (index_t k2 = 0; k2 < NE2; ++k2) {
837     for (index_t k0 = 0; k0 < NE0; ++k0) {
838     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
839     *o++ = 0.;
840     *o++ = -1.;
841     *o = 0.;
842     }
843     }
844     }
845    
846     if (m_faceOffset[3] > -1) {
847     #pragma omp for nowait
848     for (index_t k2 = 0; k2 < NE2; ++k2) {
849     for (index_t k0 = 0; k0 < NE0; ++k0) {
850     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
851     *o++ = 0.;
852     *o++ = 1.;
853     *o = 0.;
854     }
855     }
856     }
857    
858     if (m_faceOffset[4] > -1) {
859     #pragma omp for nowait
860     for (index_t k1 = 0; k1 < NE1; ++k1) {
861     for (index_t k0 = 0; k0 < NE0; ++k0) {
862     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
863     *o++ = 0.;
864     *o++ = 0.;
865     *o = -1.;
866     }
867     }
868     }
869    
870     if (m_faceOffset[5] > -1) {
871     #pragma omp for nowait
872     for (index_t k1 = 0; k1 < NE1; ++k1) {
873     for (index_t k0 = 0; k0 < NE0; ++k0) {
874     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
875     *o++ = 0.;
876     *o++ = 0.;
877     *o = 1.;
878     }
879     }
880     }
881     } // end of parallel section
882    
883     } else {
884 sshaw 5529 std::stringstream msg;
885 sshaw 5528 msg << "setToNormal: invalid function space type "
886     << out.getFunctionSpace().getTypeCode();
887     throw RipleyException(msg.str());
888     }
889     }
890    
891     void MultiBrick::setToSize(escript::Data& out) const
892     {
893     if (out.getFunctionSpace().getTypeCode() == Elements
894     || out.getFunctionSpace().getTypeCode() == ReducedElements) {
895     out.requireWrite();
896     const dim_t numQuad=out.getNumDataPointsPerSample();
897     const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]+m_dx[2]*m_dx[2]);
898     const dim_t NE = getNumElements();
899     #pragma omp parallel for
900     for (index_t k = 0; k < NE; ++k) {
901     double* o = out.getSampleDataRW(k);
902     fill(o, o+numQuad, size);
903     }
904     } else if (out.getFunctionSpace().getTypeCode() == FaceElements
905     || out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
906     out.requireWrite();
907     const dim_t numQuad=out.getNumDataPointsPerSample();
908     const dim_t NE0 = m_NE[0];
909     const dim_t NE1 = m_NE[1];
910     const dim_t NE2 = m_NE[2];
911     #pragma omp parallel
912     {
913     if (m_faceOffset[0] > -1) {
914     const double size=min(m_dx[1],m_dx[2]);
915     #pragma omp for nowait
916     for (index_t k2 = 0; k2 < NE2; ++k2) {
917     for (index_t k1 = 0; k1 < NE1; ++k1) {
918     double* o = out.getSampleDataRW(m_faceOffset[0]+INDEX2(k1,k2,m_NE[1]));
919     fill(o, o+numQuad, size);
920     }
921     }
922     }
923    
924     if (m_faceOffset[1] > -1) {
925     const double size=min(m_dx[1],m_dx[2]);
926     #pragma omp for nowait
927     for (index_t k2 = 0; k2 < NE2; ++k2) {
928     for (index_t k1 = 0; k1 < NE1; ++k1) {
929     double* o = out.getSampleDataRW(m_faceOffset[1]+INDEX2(k1,k2,m_NE[1]));
930     fill(o, o+numQuad, size);
931     }
932     }
933     }
934    
935     if (m_faceOffset[2] > -1) {
936     const double size=min(m_dx[0],m_dx[2]);
937     #pragma omp for nowait
938     for (index_t k2 = 0; k2 < NE2; ++k2) {
939     for (index_t k0 = 0; k0 < NE0; ++k0) {
940     double* o = out.getSampleDataRW(m_faceOffset[2]+INDEX2(k0,k2,m_NE[0]));
941     fill(o, o+numQuad, size);
942     }
943     }
944     }
945    
946     if (m_faceOffset[3] > -1) {
947     const double size=min(m_dx[0],m_dx[2]);
948     #pragma omp for nowait
949     for (index_t k2 = 0; k2 < NE2; ++k2) {
950     for (index_t k0 = 0; k0 < NE0; ++k0) {
951     double* o = out.getSampleDataRW(m_faceOffset[3]+INDEX2(k0,k2,m_NE[0]));
952     fill(o, o+numQuad, size);
953     }
954     }
955     }
956    
957     if (m_faceOffset[4] > -1) {
958     const double size=min(m_dx[0],m_dx[1]);
959     #pragma omp for nowait
960     for (index_t k1 = 0; k1 < NE1; ++k1) {
961     for (index_t k0 = 0; k0 < NE0; ++k0) {
962     double* o = out.getSampleDataRW(m_faceOffset[4]+INDEX2(k0,k1,m_NE[0]));
963     fill(o, o+numQuad, size);
964     }
965     }
966     }
967    
968     if (m_faceOffset[5] > -1) {
969     const double size=min(m_dx[0],m_dx[1]);
970     #pragma omp for nowait
971     for (index_t k1 = 0; k1 < NE1; ++k1) {
972     for (index_t k0 = 0; k0 < NE0; ++k0) {
973     double* o = out.getSampleDataRW(m_faceOffset[5]+INDEX2(k0,k1,m_NE[0]));
974     fill(o, o+numQuad, size);
975     }
976     }
977     }
978     } // end of parallel section
979    
980     } else {
981 sshaw 5529 std::stringstream msg;
982 sshaw 5528 msg << "setToSize: invalid function space type "
983     << out.getFunctionSpace().getTypeCode();
984     throw RipleyException(msg.str());
985     }
986     }
987    
988     void MultiBrick::Print_Mesh_Info(const bool full) const
989     {
990     RipleyDomain::Print_Mesh_Info(full);
991     if (full) {
992 sshaw 5529 std::cout << " Id Coordinates" << std::endl;
993     std::cout.precision(15);
994     std::cout.setf(ios::scientific, std::ios::floatfield);
995 sshaw 5528 for (index_t i=0; i < getNumNodes(); i++) {
996 sshaw 5529 std::cout << " " << std::setw(5) << m_nodeId[i]
997 sshaw 5528 << " " << getLocalCoordinate(i%m_NN[0], 0)
998     << " " << getLocalCoordinate(i%(m_NN[0]*m_NN[1])/m_NN[0], 1)
999 sshaw 5529 << " " << getLocalCoordinate(i/(m_NN[0]*m_NN[1]), 2) << std::endl;
1000 sshaw 5528 }
1001     }
1002     }
1003    
1004     //protected
1005     IndexVector MultiBrick::getDiagonalIndices(bool upperOnly) const
1006     {
1007     IndexVector ret;
1008     // only store non-negative indices if requested
1009     if (upperOnly)
1010     ret.resize(14);
1011     else
1012     ret.resize(27);
1013    
1014     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1015     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1016     size_t idx = 0;
1017     for (int i2=-1; i2<2; i2++) {
1018     for (int i1=-1; i1<2; i1++) {
1019     for (int i0=-1; i0<2; i0++) {
1020     const int index = i2*nDOF0*nDOF1 + i1*nDOF0 + i0;
1021     if (!upperOnly || index >= 0)
1022     ret[idx++] = index;
1023     }
1024     }
1025     }
1026    
1027     return ret;
1028     }
1029    
1030     //private
1031     void MultiBrick::populateSampleIds()
1032     {
1033     // degrees of freedom are numbered from left to right, bottom to top, front
1034     // to back in each rank, continuing on the next rank (ranks also go
1035     // left-right, bottom-top, front-back).
1036     // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
1037     // helps when writing out data rank after rank.
1038    
1039     // build node distribution vector first.
1040     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
1041     // constant for all ranks in this implementation
1042     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
1043     const dim_t numDOF=getNumDOF();
1044     for (dim_t k=1; k<m_mpiInfo->size; k++) {
1045     m_nodeDistribution[k]=k*numDOF;
1046     }
1047     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
1048    
1049     try {
1050     m_nodeId.resize(getNumNodes());
1051     m_dofId.resize(numDOF);
1052     m_elementId.resize(getNumElements());
1053 sshaw 5529 } catch (const std::length_error& le) {
1054 sshaw 5528 throw RipleyException("The system does not have sufficient memory for a domain of this size.");
1055     }
1056    
1057     // populate face element counts
1058     //left
1059     if (m_offset[0]==0)
1060     m_faceCount[0]=m_NE[1]*m_NE[2];
1061     else
1062     m_faceCount[0]=0;
1063     //right
1064     if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
1065     m_faceCount[1]=m_NE[1]*m_NE[2];
1066     else
1067     m_faceCount[1]=0;
1068     //bottom
1069     if (m_offset[1]==0)
1070     m_faceCount[2]=m_NE[0]*m_NE[2];
1071     else
1072     m_faceCount[2]=0;
1073     //top
1074     if (m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0]==m_NX[1]-1)
1075     m_faceCount[3]=m_NE[0]*m_NE[2];
1076     else
1077     m_faceCount[3]=0;
1078     //front
1079     if (m_offset[2]==0)
1080     m_faceCount[4]=m_NE[0]*m_NE[1];
1081     else
1082     m_faceCount[4]=0;
1083     //back
1084     if (m_mpiInfo->rank/(m_NX[0]*m_NX[1])==m_NX[2]-1)
1085     m_faceCount[5]=m_NE[0]*m_NE[1];
1086     else
1087     m_faceCount[5]=0;
1088    
1089     const dim_t NFE = getNumFaceElements();
1090     m_faceId.resize(NFE);
1091    
1092     const index_t left = (m_offset[0]==0 ? 0 : 1);
1093     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1094     const index_t front = (m_offset[2]==0 ? 0 : 1);
1095     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1096     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1097     const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
1098     const dim_t NN0 = m_NN[0];
1099     const dim_t NN1 = m_NN[1];
1100     const dim_t NN2 = m_NN[2];
1101     const dim_t NE0 = m_NE[0];
1102     const dim_t NE1 = m_NE[1];
1103     const dim_t NE2 = m_NE[2];
1104    
1105     // the following is a compromise between efficiency and code length to
1106     // set the node id's according to the order mentioned above.
1107     // First we set all the edge and corner id's in a rather slow way since
1108     // they might or might not be owned by this rank. Next come the own
1109     // node id's which are identical to the DOF id's (simple loop), and finally
1110     // the 6 faces are set but only if required...
1111    
1112     #define globalNodeId(x,y,z) \
1113     ((m_offset[0]+x)/nDOF0)*nDOF0*nDOF1*nDOF2+(m_offset[0]+x)%nDOF0\
1114     + ((m_offset[1]+y)/nDOF1)*nDOF0*nDOF1*nDOF2*m_NX[0]+((m_offset[1]+y)%nDOF1)*nDOF0\
1115     + ((m_offset[2]+z)/nDOF2)*nDOF0*nDOF1*nDOF2*m_NX[0]*m_NX[1]+((m_offset[2]+z)%nDOF2)*nDOF0*nDOF1
1116    
1117     #pragma omp parallel
1118     {
1119     // set edge id's
1120     // edges in x-direction, including corners
1121     #pragma omp for nowait
1122     for (dim_t i=0; i<NN0; i++) {
1123     m_nodeId[i] = globalNodeId(i, 0, 0); // LF
1124     m_nodeId[NN0*(NN1-1)+i] = globalNodeId(i, NN1-1, 0); // UF
1125     m_nodeId[NN0*NN1*(NN2-1)+i] = globalNodeId(i, 0, NN2-1); // LB
1126     m_nodeId[NN0*NN1*NN2-NN0+i] = globalNodeId(i, NN1-1, NN2-1); // UB
1127     }
1128     // edges in y-direction, without corners
1129     #pragma omp for nowait
1130     for (dim_t i=1; i<NN1-1; i++) {
1131     m_nodeId[NN0*i] = globalNodeId(0, i, 0); // FL
1132     m_nodeId[NN0*(i+1)-1] = globalNodeId(NN0-1, i, 0); // FR
1133     m_nodeId[NN0*NN1*(NN2-1)+NN0*i] = globalNodeId(0, i, NN2-1); // BL
1134     m_nodeId[NN0*NN1*(NN2-1)+NN0*(i+1)-1] = globalNodeId(NN0-1, i, NN2-1); // BR
1135     }
1136     // edges in z-direction, without corners
1137     #pragma omp for
1138     for (dim_t i=1; i<NN2-1; i++) {
1139     m_nodeId[NN0*NN1*i] = globalNodeId(0, 0, i); // LL
1140     m_nodeId[NN0*NN1*i+NN0-1] = globalNodeId(NN0-1, 0, i); // LR
1141     m_nodeId[NN0*NN1*(i+1)-NN0] = globalNodeId(0, NN1-1, i); // UL
1142     m_nodeId[NN0*NN1*(i+1)-1] = globalNodeId(NN0-1, NN1-1, i); // UR
1143     }
1144     // implicit barrier here because some node IDs will be overwritten
1145     // below
1146    
1147     // populate degrees of freedom and own nodes (identical id)
1148     #pragma omp for nowait
1149     for (dim_t i=0; i<nDOF2; i++) {
1150     for (dim_t j=0; j<nDOF1; j++) {
1151     for (dim_t k=0; k<nDOF0; k++) {
1152     const index_t nodeIdx=k+left+(j+bottom)*NN0+(i+front)*NN0*NN1;
1153     const index_t dofIdx=k+j*nDOF0+i*nDOF0*nDOF1;
1154     m_dofId[dofIdx] = m_nodeId[nodeIdx]
1155     = m_nodeDistribution[m_mpiInfo->rank]+dofIdx;
1156     }
1157     }
1158     }
1159    
1160     // populate the rest of the nodes (shared with other ranks)
1161     if (m_faceCount[0]==0) { // left plane
1162     #pragma omp for nowait
1163     for (dim_t i=0; i<nDOF2; i++) {
1164     for (dim_t j=0; j<nDOF1; j++) {
1165     const index_t nodeIdx=(j+bottom)*NN0+(i+front)*NN0*NN1;
1166     const index_t dofId=(j+1)*nDOF0-1+i*nDOF0*nDOF1;
1167     m_nodeId[nodeIdx]
1168     = m_nodeDistribution[m_mpiInfo->rank-1]+dofId;
1169     }
1170     }
1171     }
1172     if (m_faceCount[1]==0) { // right plane
1173     #pragma omp for nowait
1174     for (dim_t i=0; i<nDOF2; i++) {
1175     for (dim_t j=0; j<nDOF1; j++) {
1176     const index_t nodeIdx=(j+bottom+1)*NN0-1+(i+front)*NN0*NN1;
1177     const index_t dofId=j*nDOF0+i*nDOF0*nDOF1;
1178     m_nodeId[nodeIdx]
1179     = m_nodeDistribution[m_mpiInfo->rank+1]+dofId;
1180     }
1181     }
1182     }
1183     if (m_faceCount[2]==0) { // bottom plane
1184     #pragma omp for nowait
1185     for (dim_t i=0; i<nDOF2; i++) {
1186     for (dim_t k=0; k<nDOF0; k++) {
1187     const index_t nodeIdx=k+left+(i+front)*NN0*NN1;
1188     const index_t dofId=nDOF0*(nDOF1-1)+k+i*nDOF0*nDOF1;
1189     m_nodeId[nodeIdx]
1190     = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]]+dofId;
1191     }
1192     }
1193     }
1194     if (m_faceCount[3]==0) { // top plane
1195     #pragma omp for nowait
1196     for (dim_t i=0; i<nDOF2; i++) {
1197     for (dim_t k=0; k<nDOF0; k++) {
1198     const index_t nodeIdx=k+left+(i+front)*NN0*NN1+NN0*(NN1-1);
1199     const index_t dofId=k+i*nDOF0*nDOF1;
1200     m_nodeId[nodeIdx]
1201     = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]]+dofId;
1202     }
1203     }
1204     }
1205     if (m_faceCount[4]==0) { // front plane
1206     #pragma omp for nowait
1207     for (dim_t j=0; j<nDOF1; j++) {
1208     for (dim_t k=0; k<nDOF0; k++) {
1209     const index_t nodeIdx=k+left+(j+bottom)*NN0;
1210     const index_t dofId=k+j*nDOF0+nDOF0*nDOF1*(nDOF2-1);
1211     m_nodeId[nodeIdx]
1212     = m_nodeDistribution[m_mpiInfo->rank-m_NX[0]*m_NX[1]]+dofId;
1213     }
1214     }
1215     }
1216     if (m_faceCount[5]==0) { // back plane
1217     #pragma omp for nowait
1218     for (dim_t j=0; j<nDOF1; j++) {
1219     for (dim_t k=0; k<nDOF0; k++) {
1220     const index_t nodeIdx=k+left+(j+bottom)*NN0+NN0*NN1*(NN2-1);
1221     const index_t dofId=k+j*nDOF0;
1222     m_nodeId[nodeIdx]
1223     = m_nodeDistribution[m_mpiInfo->rank+m_NX[0]*m_NX[1]]+dofId;
1224     }
1225     }
1226     }
1227    
1228     // populate element id's
1229     #pragma omp for nowait
1230     for (dim_t i2=0; i2<NE2; i2++) {
1231     for (dim_t i1=0; i1<NE1; i1++) {
1232     for (dim_t i0=0; i0<NE0; i0++) {
1233     m_elementId[i0+i1*NE0+i2*NE0*NE1] =
1234     (m_offset[2]+i2)*m_gNE[0]*m_gNE[1]
1235     +(m_offset[1]+i1)*m_gNE[0]
1236     +m_offset[0]+i0;
1237     }
1238     }
1239     }
1240    
1241     // face elements
1242     #pragma omp for
1243     for (dim_t k=0; k<NFE; k++)
1244     m_faceId[k]=k;
1245     } // end parallel section
1246    
1247     #undef globalNodeId
1248    
1249     m_nodeTags.assign(getNumNodes(), 0);
1250     updateTagsInUse(Nodes);
1251    
1252     m_elementTags.assign(getNumElements(), 0);
1253     updateTagsInUse(Elements);
1254    
1255     // generate face offset vector and set face tags
1256     const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20, FRONT=100, BACK=200;
1257     const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP, FRONT, BACK };
1258     m_faceOffset.assign(6, -1);
1259     m_faceTags.clear();
1260     index_t offset=0;
1261     for (size_t i=0; i<6; i++) {
1262     if (m_faceCount[i]>0) {
1263     m_faceOffset[i]=offset;
1264     offset+=m_faceCount[i];
1265     m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
1266     }
1267     }
1268     setTagMap("left", LEFT);
1269     setTagMap("right", RIGHT);
1270     setTagMap("bottom", BOTTOM);
1271     setTagMap("top", TOP);
1272     setTagMap("front", FRONT);
1273     setTagMap("back", BACK);
1274     updateTagsInUse(FaceElements);
1275    
1276     populateDofMap();
1277     }
1278    
1279     //private
1280     vector<IndexVector> MultiBrick::getConnections() const
1281     {
1282     // returns a vector v of size numDOF where v[i] is a vector with indices
1283     // of DOFs connected to i (up to 27 in 3D)
1284     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1285     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1286     const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
1287     const dim_t M = nDOF0*nDOF1*nDOF2;
1288     vector<IndexVector> indices(M);
1289    
1290     #pragma omp parallel for
1291     for (index_t i=0; i < M; i++) {
1292     const index_t x = i % nDOF0;
1293     const index_t y = i % (nDOF0*nDOF1)/nDOF0;
1294     const index_t z = i / (nDOF0*nDOF1);
1295     // loop through potential neighbours and add to index if positions are
1296     // within bounds
1297     for (int i2=z-1; i2<z+2; i2++) {
1298     for (int i1=y-1; i1<y+2; i1++) {
1299     for (int i0=x-1; i0<x+2; i0++) {
1300     if (i0>=0 && i1>=0 && i2>=0
1301     && i0<nDOF0 && i1<nDOF1 && i2<nDOF2) {
1302     indices[i].push_back(i2*nDOF0*nDOF1 + i1*nDOF0 + i0);
1303     }
1304     }
1305     }
1306     }
1307     }
1308     return indices;
1309     }
1310    
1311     //private
1312     void MultiBrick::populateDofMap()
1313     {
1314     const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
1315     const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
1316     const dim_t nDOF2 = (m_gNE[2]+1)/m_NX[2];
1317     const index_t left = (m_offset[0]==0 ? 0 : 1);
1318     const index_t bottom = (m_offset[1]==0 ? 0 : 1);
1319     const index_t front = (m_offset[2]==0 ? 0 : 1);
1320    
1321     // populate node->DOF mapping with own degrees of freedom.
1322     // The rest is assigned in the loop further down
1323     m_dofMap.assign(getNumNodes(), 0);
1324     #pragma omp parallel for
1325     for (index_t i=front; i<front+nDOF2; i++) {
1326     for (index_t j=bottom; j<bottom+nDOF1; j++) {
1327     for (index_t k=left; k<left+nDOF0; k++) {
1328     m_dofMap[i*m_NN[0]*m_NN[1]+j*m_NN[0]+k]=(i-front)*nDOF0*nDOF1+(j-bottom)*nDOF0+k-left;
1329     }
1330     }
1331     }
1332    
1333     const dim_t numDOF=nDOF0*nDOF1*nDOF2;
1334     RankVector neighbour;
1335     IndexVector offsetInShared(1,0);
1336     IndexVector sendShared, recvShared;
1337     dim_t numShared=0;
1338     const int x=m_mpiInfo->rank%m_NX[0];
1339     const int y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
1340     const int z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
1341    
1342     // build list of shared components and neighbours by looping through
1343     // all potential neighbouring ranks and checking if positions are
1344     // within bounds
1345     for (int i2=-1; i2<2; i2++) {
1346     for (int i1=-1; i1<2; i1++) {
1347     for (int i0=-1; i0<2; i0++) {
1348     // skip this rank
1349     if (i0==0 && i1==0 && i2==0)
1350     continue;
1351     // location of neighbour rank
1352     const int nx=x+i0;
1353     const int ny=y+i1;
1354     const int nz=z+i2;
1355     if (nx>=0 && ny>=0 && nz>=0 && nx<m_NX[0] && ny<m_NX[1] && nz<m_NX[2]) {
1356     neighbour.push_back(nz*m_NX[0]*m_NX[1]+ny*m_NX[0]+nx);
1357     if (i0==0 && i1==0) {
1358     // sharing front or back plane
1359     offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF1);
1360     for (dim_t i=0; i<nDOF1; i++) {
1361     const dim_t firstDOF=(i2==-1 ? i*nDOF0
1362     : i*nDOF0 + nDOF0*nDOF1*(nDOF2-1));
1363     const dim_t firstNode=(i2==-1 ? left+(i+bottom)*m_NN[0]
1364     : left+(i+bottom)*m_NN[0]+m_NN[0]*m_NN[1]*(m_NN[2]-1));
1365     for (dim_t j=0; j<nDOF0; j++, numShared++) {
1366     sendShared.push_back(firstDOF+j);
1367     recvShared.push_back(numDOF+numShared);
1368     m_dofMap[firstNode+j]=numDOF+numShared;
1369     }
1370     }
1371     } else if (i0==0 && i2==0) {
1372     // sharing top or bottom plane
1373     offsetInShared.push_back(offsetInShared.back()+nDOF0*nDOF2);
1374     for (dim_t i=0; i<nDOF2; i++) {
1375     const dim_t firstDOF=(i1==-1 ? i*nDOF0*nDOF1
1376     : nDOF0*((i+1)*nDOF1-1));
1377     const dim_t firstNode=(i1==-1 ?
1378     left+(i+front)*m_NN[0]*m_NN[1]
1379     : left+m_NN[0]*((i+1+front)*m_NN[1]-1));
1380     for (dim_t j=0; j<nDOF0; j++, numShared++) {
1381     sendShared.push_back(firstDOF+j);
1382     recvShared.push_back(numDOF+numShared);
1383     m_dofMap[firstNode+j]=numDOF+numShared;
1384     }
1385     }
1386     } else if (i1==0 && i2==0) {
1387     // sharing left or right plane
1388     offsetInShared.push_back(offsetInShared.back()+nDOF1*nDOF2);
1389     for (dim_t i=0; i<nDOF2; i++) {
1390     const dim_t firstDOF=(i0==-1 ? i*nDOF0*nDOF1
1391     : nDOF0*(1+i*nDOF1)-1);
1392     const dim_t firstNode=(i0==-1 ?
1393     (bottom+(i+front)*m_NN[1])*m_NN[0]
1394     : (bottom+1+(i+front)*m_NN[1])*m_NN[0]-1);
1395     for (dim_t j=0; j<nDOF1; j++, numShared++) {
1396     sendShared.push_back(firstDOF+j*nDOF0);
1397     recvShared.push_back(numDOF+numShared);
1398     m_dofMap[firstNode+j*m_NN[0]]=numDOF+numShared;
1399     }
1400     }
1401     } else if (i0==0) {
1402     // sharing an edge in x direction
1403     offsetInShared.push_back(offsetInShared.back()+nDOF0);
1404     const dim_t firstDOF=(i1+1)/2*nDOF0*(nDOF1-1)
1405     +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
1406     const dim_t firstNode=left+(i1+1)/2*m_NN[0]*(m_NN[1]-1)
1407     +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
1408     for (dim_t i=0; i<nDOF0; i++, numShared++) {
1409     sendShared.push_back(firstDOF+i);
1410     recvShared.push_back(numDOF+numShared);
1411     m_dofMap[firstNode+i]=numDOF+numShared;
1412     }
1413     } else if (i1==0) {
1414     // sharing an edge in y direction
1415     offsetInShared.push_back(offsetInShared.back()+nDOF1);
1416     const dim_t firstDOF=(i0+1)/2*(nDOF0-1)
1417     +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
1418     const dim_t firstNode=bottom*m_NN[0]
1419     +(i0+1)/2*(m_NN[0]-1)
1420     +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
1421     for (dim_t i=0; i<nDOF1; i++, numShared++) {
1422     sendShared.push_back(firstDOF+i*nDOF0);
1423     recvShared.push_back(numDOF+numShared);
1424     m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1425     }
1426     } else if (i2==0) {
1427     // sharing an edge in z direction
1428     offsetInShared.push_back(offsetInShared.back()+nDOF2);
1429     const dim_t firstDOF=(i0+1)/2*(nDOF0-1)
1430     +(i1+1)/2*nDOF0*(nDOF1-1);
1431     const dim_t firstNode=front*m_NN[0]*m_NN[1]
1432     +(i0+1)/2*(m_NN[0]-1)
1433     +(i1+1)/2*m_NN[0]*(m_NN[1]-1);
1434     for (dim_t i=0; i<nDOF2; i++, numShared++) {
1435     sendShared.push_back(firstDOF+i*nDOF0*nDOF1);
1436     recvShared.push_back(numDOF+numShared);
1437     m_dofMap[firstNode+i*m_NN[0]*m_NN[1]]=numDOF+numShared;
1438     }
1439     } else {
1440     // sharing a node
1441     const dim_t dof = (i0+1)/2*(nDOF0-1)
1442     +(i1+1)/2*nDOF0*(nDOF1-1)
1443     +(i2+1)/2*nDOF0*nDOF1*(nDOF2-1);
1444     const dim_t node = (i0+1)/2*(m_NN[0]-1)
1445     +(i1+1)/2*m_NN[0]*(m_NN[1]-1)
1446     +(i2+1)/2*m_NN[0]*m_NN[1]*(m_NN[2]-1);
1447     offsetInShared.push_back(offsetInShared.back()+1);
1448     sendShared.push_back(dof);
1449     recvShared.push_back(numDOF+numShared);
1450     m_dofMap[node] = numDOF+numShared;
1451     ++numShared;
1452     }
1453     }
1454     }
1455     }
1456     }
1457    
1458     // TODO: paso::SharedComponents should take vectors to avoid this
1459     Esys_MPI_rank* neighPtr = NULL;
1460     index_t* sendPtr = NULL;
1461     index_t* recvPtr = NULL;
1462     if (neighbour.size() > 0) {
1463     neighPtr = &neighbour[0];
1464     sendPtr = &sendShared[0];
1465     recvPtr = &recvShared[0];
1466     }
1467     // create connector
1468     paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
1469     numDOF, neighbour.size(), neighPtr, sendPtr,
1470     &offsetInShared[0], 1, 0, m_mpiInfo));
1471     paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
1472     numDOF, neighbour.size(), neighPtr, recvPtr,
1473     &offsetInShared[0], 1, 0, m_mpiInfo));
1474     m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
1475    
1476     // useful debug output
1477     /*
1478 sshaw 5529 std::cout << "--- rcv_shcomp ---" << std::endl;
1479     std::cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << std::endl;
1480 sshaw 5528 for (size_t i=0; i<neighbour.size(); i++) {
1481 sshaw 5529 std::cout << "neighbor[" << i << "]=" << neighbour[i]
1482     << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << std::endl;
1483 sshaw 5528 }
1484     for (size_t i=0; i<recvShared.size(); i++) {
1485 sshaw 5529 std::cout << "shared[" << i << "]=" << recvShared[i] << std::endl;
1486 sshaw 5528 }
1487 sshaw 5529 std::cout << "--- snd_shcomp ---" << std::endl;
1488 sshaw 5528 for (size_t i=0; i<sendShared.size(); i++) {
1489 sshaw 5529 std::cout << "shared[" << i << "]=" << sendShared[i] << std::endl;
1490 sshaw 5528 }
1491 sshaw 5529 std::cout << "--- dofMap ---" << std::endl;
1492 sshaw 5528 for (size_t i=0; i<m_dofMap.size(); i++) {
1493 sshaw 5529 std::cout << "m_dofMap[" << i << "]=" << m_dofMap[i] << std::endl;
1494 sshaw 5528 }
1495     */
1496     }
1497    
1498 sshaw 5540 RankVector MultiBrick::getOwnerVector(int fsType) const
1499     {
1500     if (m_subdivisions != 1)
1501     throw RipleyException("Multiresolution domains only support ownership for the coarsest level");
1502     return Brick::getOwnerVector(fsType);
1503     }
1504    
1505 sshaw 5528 dim_t MultiBrick::findNode(const double *coords) const
1506     {
1507     const dim_t NOT_MINE = -1;
1508     //is the found element even owned by this rank
1509     // (inside owned or shared elements but will map to an owned element)
1510     for (int dim = 0; dim < m_numDim; dim++) {
1511     double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
1512     - m_dx[dim]/2.; //allows for point outside mapping onto node
1513     double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
1514     + m_dx[dim]/2.;
1515     if (min > coords[dim] || max < coords[dim]) {
1516     return NOT_MINE;
1517     }
1518     }
1519     // get distance from origin
1520     double x = coords[0] - m_origin[0];
1521     double y = coords[1] - m_origin[1];
1522     double z = coords[2] - m_origin[2];
1523    
1524     //check if the point is even inside the domain
1525     if (x < 0 || y < 0 || z < 0
1526     || x > m_length[0] || y > m_length[1] || z > m_length[2])
1527     return NOT_MINE;
1528    
1529     // distance in elements
1530     dim_t ex = (dim_t) floor(x / m_dx[0]);
1531     dim_t ey = (dim_t) floor(y / m_dx[1]);
1532     dim_t ez = (dim_t) floor(z / m_dx[2]);
1533     // set the min distance high enough to be outside the element plus a bit
1534     dim_t closest = NOT_MINE;
1535     double minDist = 1;
1536     for (int dim = 0; dim < m_numDim; dim++) {
1537     minDist += m_dx[dim]*m_dx[dim];
1538     }
1539     //find the closest node
1540     for (int dx = 0; dx < 1; dx++) {
1541     double xdist = x - (ex + dx)*m_dx[0];
1542     for (int dy = 0; dy < 1; dy++) {
1543     double ydist = y - (ey + dy)*m_dx[1];
1544     for (int dz = 0; dz < 1; dz++) {
1545     double zdist = z - (ez + dz)*m_dx[2];
1546     double total = xdist*xdist + ydist*ydist + zdist*zdist;
1547     if (total < minDist) {
1548     closest = INDEX3(ex+dy-m_offset[0], ey+dy-m_offset[1],
1549     ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
1550     minDist = total;
1551     }
1552     }
1553     }
1554     }
1555     if (closest == NOT_MINE) {
1556     throw RipleyException("Unable to map appropriate dirac point to a "
1557     "node, implementation problem in MultiBrick::findNode()");
1558     }
1559     return closest;
1560     }
1561    
1562     } // end of namespace ripley
1563    

  ViewVC Help
Powered by ViewVC 1.1.26