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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26