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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6141 - (show annotations)
Wed Apr 6 03:51:30 2016 UTC (2 years, 6 months ago) by caltinay
File size: 68739 byte(s)
more namespacing of defines.

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

  ViewVC Help
Powered by ViewVC 1.1.26