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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6141 - (show annotations)
Wed Apr 6 03:51:30 2016 UTC (2 years, 8 months ago) by caltinay
File size: 38751 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/MultiRectangle.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 <algorithm>
28 #include <iomanip>
29 #include <iostream>
30 #include <limits>
31
32 using std::vector;
33 using std::string;
34
35 namespace ripley {
36
37 MultiRectangle::MultiRectangle(dim_t n0, dim_t n1, double x0, double y0,
38 double x1, double y1, int d0, int d1,
39 const vector<double>& points,
40 const vector<int>& tags,
41 const TagMap& tagnamestonums,
42 escript::SubWorld_ptr w, unsigned int subdivisions)
43 : Rectangle(n0,n1, x0,y0, x1,y1, d0,d1, points, tags, tagnamestonums, w),
44 m_subdivisions(subdivisions)
45 {
46 if (subdivisions == 0 || (subdivisions & (subdivisions - 1)) != 0)
47 throw RipleyException("Element subdivisions must be a power of two");
48
49 dim_t oldNN[2] = {0};
50
51 if (d0 == 0 || d1 == 0)
52 throw RipleyException("Domain subdivisions must be positive");
53
54 for (int i = 0; i < 2; i++) {
55 m_NE[i] *= subdivisions;
56 oldNN[i] = m_NN[i];
57 m_NN[i] = m_NE[i] + 1;
58 m_gNE[i] *= subdivisions;
59 m_ownNE[i] *= subdivisions;
60 m_dx[i] /= subdivisions;
61 m_faceCount[i] *= subdivisions;
62 m_faceCount[2+i] *= subdivisions;
63 m_offset[i] *= subdivisions;
64 }
65 populateSampleIds();
66
67 const dim_t nDirac = m_diracPoints.size();
68 #pragma omp parallel for
69 for (int i = 0; i < nDirac; i++) {
70 const dim_t node = m_diracPoints[i].node;
71 const dim_t x = node % oldNN[0];
72 const dim_t y = node / oldNN[0];
73 m_diracPoints[i].node = INDEX2(x*subdivisions, y*subdivisions, m_NN[0]);
74 m_diracPointNodeIDs[i] = m_diracPoints[i].node;
75 }
76 }
77
78 MultiRectangle::~MultiRectangle()
79 {
80 }
81
82 void MultiRectangle::validateInterpolationAcross(int fsType_source,
83 const escript::AbstractDomain& domain, int fsType_target) const
84 {
85 const MultiRectangle *other = dynamic_cast<const MultiRectangle *>(&domain);
86 if (other == NULL)
87 throw RipleyException("Invalid interpolation: domains must both be instances of MultiRectangle");
88
89 const double *len = other->getLength();
90 const int *subdivs = other->getNumSubdivisionsPerDim();
91 const dim_t *elements = other->getNumElementsPerDim();
92 const unsigned int level = other->getNumSubdivisionsPerElement();
93 const unsigned int factor = m_subdivisions > level ? m_subdivisions/level : level/m_subdivisions;
94 if ((factor & (factor - 1)) != 0) //factor == 2**x
95 throw RipleyException("Invalid interpolation: elemental subdivisions of each domain must be powers of two");
96
97 if (other->getMPIComm() != m_mpiInfo->comm)
98 throw RipleyException("Invalid interpolation: Domains are on different communicators");
99 for (int i = 0; i < m_numDim; i++) {
100 if (m_length[i] != len[i]) {
101 throw RipleyException("Invalid interpolation: domain length mismatch");
102 }
103 if (m_NX[i] != subdivs[i]) {
104 throw RipleyException("Invalid interpolation: domain process subdivision mismatch");
105 }
106 if (m_subdivisions > level) {
107 if (m_NE[i]/elements[i] != factor) {
108 std::cerr << "m_ownNE[i]/elements[i] = "
109 << m_ownNE[i]/elements[i] << " != " << factor << std::endl;
110 throw RipleyException("Invalid interpolation: element factor mismatch");
111 }
112 } else {
113 if (elements[i]/m_NE[i] != factor) {
114 throw RipleyException("Invalid interpolation: element factor mismatch");
115 }
116 }
117 }
118 }
119
120 void MultiRectangle::interpolateNodesToNodesFiner(const escript::Data& source,
121 escript::Data& target, const MultiRectangle& other) const
122 {
123 const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
124 const dim_t NN0 = m_NN[0], NN1 = m_NN[1], otherNN0 = other.getNumNodesPerDim()[0];
125 const dim_t numComp = source.getDataPointSize();
126 target.requireWrite();
127 #pragma omp parallel for
128 for (dim_t ny = 0; ny < NN1 - 1; ny++) { //source nodes
129 for (dim_t nx = 0; nx < NN0 - 1; nx++) {
130 const double *x0y0 = source.getSampleDataRO(ny*NN0 + nx);
131 const double *x0y1 = source.getSampleDataRO((ny+1)*NN0 + nx);
132 const double *x1y0 = source.getSampleDataRO(ny*NN0 + nx + 1);
133 const double *x1y1 = source.getSampleDataRO((ny+1)*NN0 + nx + 1);
134 const double origin[2] = {getLocalCoordinate(nx, 0), getLocalCoordinate(ny, 1)};
135 for (int sy = 0; sy < scaling + 1; sy++) { //target nodes
136 const double y = (other.getLocalCoordinate(ny*scaling+sy, 1) - origin[1]) / m_dx[1];
137 for (int sx = 0; sx < scaling + 1; sx++) {
138 const double x = (other.getLocalCoordinate(nx*scaling+sx, 0) - origin[0]) / m_dx[0];
139 double *out = target.getSampleDataRW(nx*scaling+sx + (ny*scaling+sy)*otherNN0);
140 for (int comp = 0; comp < numComp; comp++) {
141 out[comp] = x0y0[comp]*(1-x)*(1-y) + x1y0[comp]*x*(1-y) + x0y1[comp]*(1-x)*y + x1y1[comp]*x*y;
142 }
143 }
144 }
145 }
146 }
147 }
148
149 void MultiRectangle::interpolateReducedToElementsFiner(const escript::Data& source,
150 escript::Data& target, const MultiRectangle& other) const
151 {
152 const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
153 const dim_t numComp = source.getDataPointSize();
154 target.requireWrite();
155 //for each of ours
156 #pragma omp parallel for
157 for (dim_t ey = 0; ey < m_NE[1]; ey++) {
158 for (dim_t ex = 0; ex < m_NE[0]; ex++) {
159 const double *in = source.getSampleDataRO(ex + ey*m_NE[0]);
160 //for each subelement
161 for (dim_t sy = 0; sy < scaling; sy++) {
162 const dim_t ty = ey*scaling + sy;
163 for (dim_t sx = 0; sx < scaling; sx++) {
164 const dim_t tx = ex*scaling + sx;
165 double *out = target.getSampleDataRW(tx + ty*m_NE[0]*scaling);
166 for (dim_t comp = 0; comp < numComp; comp++) {
167 const double quadvalue = in[comp];
168 out[INDEX3(comp, 0, 0, numComp, 2)] = quadvalue;
169 out[INDEX3(comp, 0, 1, numComp, 2)] = quadvalue;
170 out[INDEX3(comp, 1, 0, numComp, 2)] = quadvalue;
171 out[INDEX3(comp, 1, 1, numComp, 2)] = quadvalue;
172 }
173 }
174 }
175 }
176 }
177 }
178
179 void MultiRectangle::interpolateReducedToReducedFiner(const escript::Data& source,
180 escript::Data& target, const MultiRectangle& other) const
181 {
182 const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
183 const dim_t numComp = source.getDataPointSize();
184 target.requireWrite();
185 //for each of ours
186 #pragma omp parallel for
187 for (dim_t ey = 0; ey < m_NE[1]; ey++) {
188 for (dim_t ex = 0; ex < m_NE[0]; ex++) {
189 const double *in = source.getSampleDataRO(ex + ey*m_NE[0]);
190 //for each subelement
191 for (dim_t sy = 0; sy < scaling; sy++) {
192 const dim_t ty = ey*scaling + sy;
193 for (dim_t sx = 0; sx < scaling; sx++) {
194 const dim_t tx = ex*scaling + sx;
195 double *out = target.getSampleDataRW(tx + ty*m_NE[0]*scaling);
196 for (dim_t comp = 0; comp < numComp; comp++) {
197 out[comp] = in[comp];
198 }
199 }
200 }
201 }
202 }
203 }
204
205 void MultiRectangle::interpolateNodesToElementsFiner(const escript::Data& source,
206 escript::Data& target, const MultiRectangle& other) const
207 {
208 const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
209 const dim_t NE0 = m_NE[0], NE1 = m_NE[1];
210 const dim_t numComp = source.getDataPointSize();
211 target.requireWrite();
212 #pragma omp parallel for
213 for (dim_t ey = 0; ey < NE1; ey++) { //source nodes
214 for (dim_t ex = 0; ex < NE0; ex++) {
215 const double *x0y0 = source.getSampleDataRO(ey*(NE0+1) + ex);
216 const double *x0y1 = source.getSampleDataRO((ey+1)*(NE0+1) + ex);
217 const double *x1y0 = source.getSampleDataRO(ey*(NE0+1) + ex + 1);
218 const double *x1y1 = source.getSampleDataRO((ey+1)*(NE0+1) + ex + 1);
219 const double origin[2] = {getLocalCoordinate(ex, 0), getLocalCoordinate(ey, 1)};
220 for (int sy = 0; sy < scaling; sy++) { //target elements
221 for (int sx = 0; sx < scaling; sx++) {
222 const double x1 = (other.getLocalCoordinate(ex*scaling+sx, 0) - origin[0]) / m_dx[0] + FIRST_QUAD/scaling;
223 const double x2 = x1 + (SECOND_QUAD - FIRST_QUAD)/scaling;
224 const double y1 = (other.getLocalCoordinate(ey*scaling+sy, 1) - origin[1]) / m_dx[1] + FIRST_QUAD/scaling;
225 const double y2 = y1 + (SECOND_QUAD - FIRST_QUAD)/scaling;
226 double *out = target.getSampleDataRW(ex*scaling+sx + (ey*scaling+sy)*NE0*scaling);
227 for (int comp = 0; comp < numComp; comp++) {
228 out[INDEX3(comp, 0, 0, numComp, 2)] = x0y0[comp]*(1-x1)*(1-y1) + x1y0[comp]*x1*(1-y1) + x0y1[comp]*(1-x1)*y1 + x1y1[comp]*x1*y1;
229 out[INDEX3(comp, 0, 1, numComp, 2)] = x0y0[comp]*(1-x1)*(1-y2) + x1y0[comp]*x1*(1-y2) + x0y1[comp]*(1-x1)*y2 + x1y1[comp]*x1*y2;
230 out[INDEX3(comp, 1, 0, numComp, 2)] = x0y0[comp]*(1-x2)*(1-y1) + x1y0[comp]*x2*(1-y1) + x0y1[comp]*(1-x2)*y1 + x1y1[comp]*x2*y1;
231 out[INDEX3(comp, 1, 1, numComp, 2)] = x0y0[comp]*(1-x2)*(1-y2) + x1y0[comp]*x2*(1-y2) + x0y1[comp]*(1-x2)*y2 + x1y1[comp]*x2*y2;
232 }
233 }
234 }
235 }
236 }
237 }
238
239 void MultiRectangle::interpolateElementsToElementsCoarser(const escript::Data& source,
240 escript::Data& target, const MultiRectangle& other) const
241 {
242 const int scaling = m_subdivisions/other.getNumSubdivisionsPerElement();
243 const double scaling_volume = (1./scaling)*(1./scaling);
244 const dim_t *theirNE = other.getNumElementsPerDim();
245 const dim_t numComp = source.getDataPointSize();
246
247 vector<double> points(scaling*2, 0);
248 vector<double> first_lagrange(scaling*2, 1);
249 vector<double> second_lagrange(scaling*2, 1);
250
251 for (int i = 0; i < scaling*2; i+=2) {
252 points[i] = (i/2 + FIRST_QUAD)/scaling;
253 points[i+1] = (i/2 + SECOND_QUAD)/scaling;
254 }
255
256 for (int i = 0; i < scaling*2; i++) {
257 first_lagrange[i] = (points[i] - SECOND_QUAD) / (FIRST_QUAD - SECOND_QUAD);
258 second_lagrange[i] = (points[i] - FIRST_QUAD) / (SECOND_QUAD - FIRST_QUAD);
259 }
260 target.requireWrite();
261 //for each of theirs
262 #pragma omp parallel for
263 for (dim_t ty = 0; ty < theirNE[1]; ty++) {
264 for (dim_t tx = 0; tx < theirNE[0]; tx++) {
265 double *out = target.getSampleDataRW(tx + ty*theirNE[0]);
266 //for each subelement
267 for (dim_t sy = 0; sy < scaling; sy++) {
268 const dim_t ey = ty*scaling + sy;
269 for (dim_t sx = 0; sx < scaling; sx++) {
270 const dim_t ex = tx*scaling + sx;
271 const double *in = source.getSampleDataRO(ex + ey*m_NE[0]);
272 for (int quad = 0; quad < 4; quad++) {
273 int lx = sx*2 + quad%2;
274 int ly = sy*2 + quad/2;
275 for (dim_t comp = 0; comp < numComp; comp++) {
276 const double quadvalue = scaling_volume * in[comp + quad*numComp];
277 out[INDEX3(comp, 0, 0, numComp, 2)] += quadvalue * first_lagrange[lx] * first_lagrange[ly];
278 out[INDEX3(comp, 0, 1, numComp, 2)] += quadvalue * first_lagrange[lx] * second_lagrange[ly];
279 out[INDEX3(comp, 1, 0, numComp, 2)] += quadvalue * second_lagrange[lx] * first_lagrange[ly];
280 out[INDEX3(comp, 1, 1, numComp, 2)] += quadvalue * second_lagrange[lx] * second_lagrange[ly];
281 }
282 }
283 }
284 }
285 }
286 }
287 }
288
289
290 void MultiRectangle::interpolateElementsToElementsFiner(const escript::Data& source,
291 escript::Data& target, const MultiRectangle& other) const
292 {
293 const int scaling = other.getNumSubdivisionsPerElement()/m_subdivisions;
294 const dim_t numComp = source.getDataPointSize();
295
296 vector<double> points(scaling*2, 0);
297 vector<double> lagranges(scaling*4, 1);
298
299 for (int i = 0; i < scaling*2; i+=2) {
300 points[i] = (i/2 + FIRST_QUAD)/scaling;
301 points[i+1] = (i/2 + SECOND_QUAD)/scaling;
302 }
303 for (int i = 0; i < scaling*2; i++) {
304 lagranges[i] = (points[i] - SECOND_QUAD) / (FIRST_QUAD - SECOND_QUAD);
305 lagranges[i + 2*scaling] = (points[i] - FIRST_QUAD) / (SECOND_QUAD - FIRST_QUAD);
306 }
307 target.requireWrite();
308 //for each of ours
309 #pragma omp parallel for
310 for (dim_t ey = 0; ey < m_NE[1]; ey++) {
311 for (dim_t ex = 0; ex < m_NE[0]; ex++) {
312 const double *in = source.getSampleDataRO(ex + ey*m_NE[0]);
313 //for each subelement
314 for (dim_t sy = 0; sy < scaling; sy++) {
315 const dim_t ty = ey*scaling + sy;
316 for (dim_t sx = 0; sx < scaling; sx++) {
317 const dim_t tx = ex*scaling + sx;
318 double *out = target.getSampleDataRW(tx + ty*m_NE[0]*scaling);
319 for (int quad = 0; quad < 4; quad++) {
320 const int lx = scaling*2*(quad%2) + sx*2;
321 const int ly = scaling*2*(quad/2) + sy*2;
322 for (dim_t comp = 0; comp < numComp; comp++) {
323 const double quadvalue = in[comp + quad*numComp];
324 out[INDEX3(comp, 0, 0, numComp, 2)] += quadvalue * lagranges[lx] * lagranges[ly];
325 out[INDEX3(comp, 0, 1, numComp, 2)] += quadvalue * lagranges[lx] * lagranges[ly+1];
326 out[INDEX3(comp, 1, 0, numComp, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly];
327 out[INDEX3(comp, 1, 1, numComp, 2)] += quadvalue * lagranges[lx+1] * lagranges[ly+1];
328 }
329 }
330 }
331 }
332 }
333 }
334 }
335
336 void MultiRectangle::interpolateAcross(escript::Data& target,
337 const escript::Data& source) const
338 {
339 const MultiRectangle *other =
340 dynamic_cast<const MultiRectangle *>(target.getDomain().get());
341 if (other == NULL)
342 throw RipleyException("Invalid interpolation: Domains must both be instances of MultiRectangle");
343 //shouldn't ever happen, but I want to know if it does
344 if (other == this)
345 throw RipleyException("interpolateAcross: this domain is the target");
346
347 validateInterpolationAcross(source.getFunctionSpace().getTypeCode(),
348 *(target.getDomain().get()), target.getFunctionSpace().getTypeCode());
349 int fsSource = source.getFunctionSpace().getTypeCode();
350 int fsTarget = target.getFunctionSpace().getTypeCode();
351
352 std::stringstream msg;
353 msg << "Invalid interpolation: interpolation not implemented for function space "
354 << functionSpaceTypeAsString(fsSource)
355 << " -> "
356 << functionSpaceTypeAsString(fsTarget);
357 if (other->getNumSubdivisionsPerElement() > getNumSubdivisionsPerElement()) {
358 switch (fsSource) {
359 case Nodes:
360 switch (fsTarget) {
361 case Nodes:
362 case ReducedNodes:
363 case DegreesOfFreedom:
364 case ReducedDegreesOfFreedom:
365 interpolateNodesToNodesFiner(source, target, *other);
366 return;
367 case Elements:
368 interpolateNodesToElementsFiner(source, target, *other);
369 return;
370 }
371 break;
372 case Elements:
373 switch (fsTarget) {
374 case Elements:
375 interpolateElementsToElementsFiner(source, target, *other);
376 return;
377 }
378 break;
379 case ReducedElements:
380 switch (fsTarget) {
381 case Elements:
382 interpolateReducedToElementsFiner(source, target, *other);
383 return;
384 }
385 break;
386 }
387 msg << " when target is a finer mesh";
388 } else {
389 switch (fsSource) {
390 case Nodes:
391 switch (fsTarget) {
392 case Elements:
393 escript::Data elements=escript::Vector(0., escript::function(*this), true);
394 interpolateNodesOnElements(elements, source, false);
395 interpolateElementsToElementsCoarser(elements, target, *other);
396 return;
397 }
398 break;
399 case Elements:
400 switch (fsTarget) {
401 case Elements:
402 interpolateElementsToElementsCoarser(source, target, *other);
403 return;
404 }
405 break;
406 }
407 msg << " when target is a coarser mesh";
408 }
409 throw RipleyException(msg.str());
410 }
411
412 string MultiRectangle::getDescription() const
413 {
414 return "ripley::MultiRectangle";
415 }
416
417 bool MultiRectangle::operator==(const AbstractDomain& other) const
418 {
419 const MultiRectangle* o=dynamic_cast<const MultiRectangle*>(&other);
420 if (o) {
421 return (RipleyDomain::operator==(other) &&
422 m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
423 && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
424 && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
425 && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]
426 && m_subdivisions==o->m_subdivisions);
427 }
428
429 return false;
430 }
431
432 void MultiRectangle::readNcGrid(escript::Data& out, string filename, string varname,
433 const ReaderParameters& params) const
434 {
435 if (m_subdivisions != 1)
436 throw RipleyException("Non-parent MultiRectangles cannot read datafiles");
437 Rectangle::readNcGrid(out, filename, varname, params);
438 }
439
440 void MultiRectangle::readBinaryGrid(escript::Data& out, string filename,
441 const ReaderParameters& params) const
442 {
443 if (m_subdivisions != 1)
444 throw RipleyException("Non-parent MultiRectangles cannot read datafiles");
445 Rectangle::readBinaryGrid(out, filename, params);
446 }
447
448 void MultiRectangle::readBinaryGridFromZipped(escript::Data& out, string filename,
449 const ReaderParameters& params) const
450 {
451 if (m_subdivisions != 1)
452 throw RipleyException("Non-parent MultiRectangles cannot read datafiles");
453 Rectangle::readBinaryGridFromZipped(out, filename, params);
454 }
455
456 void MultiRectangle::writeBinaryGrid(const escript::Data& in, string filename,
457 int byteOrder, int dataType) const
458 {
459 if (m_subdivisions != 1)
460 throw RipleyException("Non-parent MultiRectangles cannot read datafiles");
461 Rectangle::writeBinaryGrid(in, filename, byteOrder, dataType);
462 }
463
464 void MultiRectangle::dump(const string& fileName) const
465 {
466 if (m_subdivisions != 1)
467 throw RipleyException("Non-parent MultiRectangles dump not implemented");
468 Rectangle::dump(fileName);
469 }
470
471 void MultiRectangle::populateDofMap()
472 {
473 // build node distribution vector first.
474 // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes.
475 // Unlike regular ripley domains this is NOT constant for all ranks so
476 // we do an Allgather (we could have also computed per rank but it's a bit
477 // involved)
478 m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
479 dim_t numDOF = getNumDOF();
480 if (m_mpiInfo->size > 1) {
481 #if ESYS_MPI
482 MPI_Allgather(&numDOF, 1, MPI_DIM_T, &m_nodeDistribution[0], 1,
483 MPI_DIM_T, m_mpiInfo->comm);
484
485 // accumulate
486 dim_t accu = 0;
487 for (int rank = 0; rank < m_mpiInfo->size; rank++) {
488 const dim_t n = m_nodeDistribution[rank];
489 m_nodeDistribution[rank] = accu;
490 accu += n;
491 }
492 ESYS_ASSERT(accu == getNumDataPointsGlobal(),
493 "something went wrong computing the DOF distribution!");
494
495 m_nodeDistribution[m_mpiInfo->size] = accu;
496 #endif
497 } else {
498 m_nodeDistribution[m_mpiInfo->size] = numDOF;
499 }
500
501 // degrees of freedom are numbered from left to right, bottom to top in
502 // each rank, continuing on the next rank (ranks also go left-right,
503 // bottom-top).
504 // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
505 // helps when writing out data rank after rank.
506
507 try {
508 m_nodeId.assign(getNumNodes(), -1);
509 m_dofMap.assign(getNumNodes(), -1);
510 m_dofId.assign(numDOF, -1);
511 } catch (const std::length_error& le) {
512 throw RipleyException("The system does not have sufficient memory for a domain of this size.");
513 }
514
515 const index_t left = getFirstInDim(0);
516 const index_t bottom = getFirstInDim(1);
517 const dim_t nDOF0 = getNumDOFInAxis(0);
518 const dim_t nDOF1 = getNumDOFInAxis(1);
519 // populate node->DOF mapping, DOF IDs and own node IDs.
520 // The rest of the node IDs are communicated further down.
521 #pragma omp parallel for
522 for (dim_t i=0; i<nDOF1; i++) {
523 for (dim_t j=0; j<nDOF0; j++) {
524 const index_t nodeIdx = j+left + (i+bottom)*m_NN[0];
525 const index_t dofIdx = j + i*nDOF0;
526 m_dofMap[nodeIdx] = dofIdx;
527 m_dofId[dofIdx] = m_nodeId[nodeIdx]
528 = m_nodeDistribution[m_mpiInfo->rank] + dofIdx;
529 }
530 }
531
532 // build list of shared components and neighbours
533 m_colIndices.assign(numDOF, IndexVector());
534 m_rowIndices.assign(getNumNodes() - numDOF, IndexVector());
535
536 RankVector neighbour;
537 IndexVector offsetInSharedSend(1,0);
538 IndexVector offsetInSharedRecv(1,0);
539 IndexVector sendShared, recvShared;
540 const int x = m_mpiInfo->rank%m_NX[0];
541 const int y = m_mpiInfo->rank/m_NX[0];
542 // numShared will contain the number of shared DOFs after the following
543 // blocks
544 dim_t numShared = 0;
545 // sharing bottom edge
546 if (y > 0) {
547 neighbour.push_back((y-1)*m_NX[0] + x);
548 //joining edge, send and recv
549 offsetInSharedSend.push_back(offsetInSharedSend.back()+nDOF0);
550 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+nDOF0*m_subdivisions);
551 // add to send only
552 for (dim_t i=0; i < nDOF0; i++) {
553 sendShared.push_back(i);
554 }
555
556 for (unsigned sy = 0; sy < m_subdivisions; sy++) {
557 for (index_t i = 0; i < nDOF0; i++, numShared++) {
558 const index_t nodeIdx = left + i + sy*m_NN[0];
559 const index_t dofIdx = i;
560 recvShared.push_back(nodeIdx);
561 m_dofMap[nodeIdx] = numDOF + numShared;
562 if (i > 0)
563 doublyLink(m_colIndices, m_rowIndices, dofIdx - 1, numShared);
564 doublyLink(m_colIndices, m_rowIndices, dofIdx, numShared);
565 if (i < nDOF0 - 1)
566 doublyLink(m_colIndices, m_rowIndices, dofIdx + 1, numShared);
567 }
568 }
569 }
570 // sharing top edge
571 if (y < m_NX[1] - 1) {
572 neighbour.push_back((y+1)*m_NX[0] + x);
573 offsetInSharedSend.push_back(offsetInSharedSend.back()+nDOF0*m_subdivisions);
574 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+nDOF0);
575 // add to send only
576 for (unsigned sy = 0; sy < m_subdivisions; sy++) {
577 for (index_t i = 0; i < nDOF0; i++) {
578 sendShared.push_back(numDOF-nDOF0*(m_subdivisions - sy) + i);
579 }
580 }
581 for (index_t i = 0; i < nDOF0; i++, numShared++) {
582 const index_t nodeIdx = left + i + m_NN[0]*(m_NN[1]-1);
583 const index_t dofIdx = numDOF - nDOF0 + i;
584 recvShared.push_back(nodeIdx);
585 m_dofMap[nodeIdx] = numDOF+numShared;
586 if (i > 0)
587 doublyLink(m_colIndices, m_rowIndices, dofIdx - 1, numShared);
588 doublyLink(m_colIndices, m_rowIndices, dofIdx, numShared);
589 if (i < nDOF0 - 1)
590 doublyLink(m_colIndices, m_rowIndices, dofIdx + 1, numShared);
591 }
592 }
593 // sharing left edge
594 if (x > 0) {
595 neighbour.push_back(y*m_NX[0] + x-1);
596 offsetInSharedSend.push_back(offsetInSharedSend.back()+nDOF1);
597 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+nDOF1*m_subdivisions);
598 for (index_t i = 0; i < nDOF1; i++) {
599 const index_t dofIdx = i*nDOF0;
600 sendShared.push_back(dofIdx);
601 for (unsigned sx = 0; sx < m_subdivisions; sx++, numShared++) {
602 const index_t nodeIdx = (bottom+i)*m_NN[0] + sx;
603 recvShared.push_back(nodeIdx);
604 m_dofMap[nodeIdx] = numDOF + numShared;
605 if (i > 0)
606 doublyLink(m_colIndices, m_rowIndices, dofIdx - nDOF0, numShared);
607 doublyLink(m_colIndices, m_rowIndices, dofIdx, numShared);
608 if (i < nDOF1 - 1)
609 doublyLink(m_colIndices, m_rowIndices, dofIdx + nDOF0, numShared);
610 }
611 }
612 }
613 // sharing right edge
614 if (x < m_NX[0] - 1) {
615 neighbour.push_back(y*m_NX[0] + x+1);
616 offsetInSharedSend.push_back(offsetInSharedSend.back()+nDOF1*m_subdivisions);
617 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+nDOF1);
618 for (index_t i = 0; i < nDOF1; i++, numShared++) {
619 for (unsigned sx = 0; sx < m_subdivisions - 1; sx++) {
620 sendShared.push_back((i+1)*nDOF0-(m_subdivisions - sx));
621 }
622 const index_t nodeIdx = (bottom+1+i)*m_NN[0] - 1;
623 const index_t dofIdx = (i+1)*nDOF0 - 1;
624 sendShared.push_back(dofIdx);
625 recvShared.push_back(nodeIdx);
626 m_dofMap[nodeIdx] = numDOF + numShared;
627 if (i > 0)
628 doublyLink(m_colIndices, m_rowIndices, dofIdx - nDOF0, numShared);
629 doublyLink(m_colIndices, m_rowIndices, dofIdx, numShared);
630 if (i < nDOF1 - 1)
631 doublyLink(m_colIndices, m_rowIndices, dofIdx + nDOF0, numShared);
632 }
633 }
634 // sharing bottom-left block
635 if (x > 0 && y > 0) {
636 neighbour.push_back((y-1)*m_NX[0] + x-1);
637 // sharing a node
638 offsetInSharedSend.push_back(offsetInSharedSend.back()+1);
639 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+m_subdivisions*m_subdivisions);
640 for (unsigned sy = 0; sy < m_subdivisions; sy++) {
641 for (unsigned sx = 0; sx < m_subdivisions; sx++, numShared++) {
642 const index_t nodeIdx = sx + sy*m_NN[0];
643 m_dofMap[nodeIdx] = numDOF + numShared;
644 recvShared.push_back(nodeIdx);
645 doublyLink(m_colIndices, m_rowIndices, 0, numShared);
646 }
647 }
648 sendShared.push_back(0);
649 }
650 // sharing top-left block
651 if (x > 0 && y < m_NX[1]-1) {
652 neighbour.push_back((y+1)*m_NX[0] + x-1);
653 offsetInSharedSend.push_back(offsetInSharedSend.back()+m_subdivisions);
654 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+m_subdivisions);
655 for (int s = 0; s < m_subdivisions; s++, numShared++) {
656 const index_t nodeIdx = m_NN[0]*(m_NN[1]-1) + s;
657 const index_t dofIdx = numDOF - (m_subdivisions - s)*nDOF0;
658 sendShared.push_back(dofIdx);
659 recvShared.push_back(nodeIdx);
660 m_dofMap[nodeIdx] = numDOF + numShared;
661 if (s > 0)
662 doublyLink(m_colIndices, m_rowIndices, dofIdx - nDOF0, numShared);
663 doublyLink(m_colIndices, m_rowIndices, dofIdx, numShared);
664 if (s < m_subdivisions - 1)
665 doublyLink(m_colIndices, m_rowIndices, dofIdx + nDOF0, numShared);
666 }
667 }
668 // sharing bottom-right block
669 if (x < m_NX[0]-1 && y > 0) {
670 neighbour.push_back((y-1)*m_NX[0] + x+1);
671 offsetInSharedSend.push_back(offsetInSharedSend.back()+m_subdivisions);
672 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+m_subdivisions);
673 for (int s = 0; s < m_subdivisions; s++, numShared++) {
674 const index_t nodeIdx = (s+1)*m_NN[0] - 1;
675 const dim_t ind = nDOF0 - (m_subdivisions - s);
676 recvShared.push_back(nodeIdx);
677 m_dofMap[nodeIdx] = numDOF + numShared;
678 sendShared.push_back(ind);
679 if (s > 0)
680 doublyLink(m_colIndices, m_rowIndices, ind - 1, numShared);
681 doublyLink(m_colIndices, m_rowIndices, ind, numShared);
682 if (s < m_subdivisions - 1)
683 doublyLink(m_colIndices, m_rowIndices, ind + 1, numShared);
684 }
685 }
686 // sharing top-right block
687 if (x < m_NX[0]-1 && y < m_NX[1]-1) {
688 neighbour.push_back((y+1)*m_NX[0] + x+1);
689 offsetInSharedSend.push_back(offsetInSharedSend.back()+m_subdivisions*m_subdivisions);
690 offsetInSharedRecv.push_back(offsetInSharedRecv.back()+1);
691 for (unsigned sy = 0; sy < m_subdivisions; sy++) {
692 for (unsigned sx = 0; sx < m_subdivisions; sx++) {
693 sendShared.push_back(numDOF-(m_subdivisions - sy - 1)*nDOF0 - (m_subdivisions - sx));
694 }
695 }
696 const dim_t nodeIdx = m_NN[0]*m_NN[1] - 1;
697 recvShared.push_back(nodeIdx);
698 m_dofMap[nodeIdx] = numDOF+numShared;
699 doublyLink(m_colIndices, m_rowIndices, numDOF-1, numShared);
700 ++numShared;
701 }
702
703 #ifdef ESYS_MPI
704 if (m_mpiInfo->size > 1) {
705 // now send off shared DOF IDs so nodeId will become a global node
706 // labelling.
707 const dim_t numSend = offsetInSharedSend.back();
708 const dim_t numRecv = offsetInSharedRecv.back();
709 IndexVector recvBuffer(numRecv);
710 IndexVector sendBuffer(numSend);
711 std::vector<MPI_Request> reqs(2*neighbour.size());
712 std::vector<MPI_Status> stats(2*neighbour.size());
713
714 // prepare the send buffer
715 #pragma omp parallel for
716 for (index_t i = 0; i < numSend; ++i) {
717 sendBuffer[i] = m_dofId[sendShared[i]];
718 }
719 for (index_t i = 0; i < neighbour.size(); i++) {
720 MPI_Irecv(&recvBuffer[offsetInSharedRecv[i]],
721 offsetInSharedRecv[i+1] - offsetInSharedRecv[i],
722 MPI_DIM_T, neighbour[i],
723 m_mpiInfo->counter()+neighbour[i],
724 m_mpiInfo->comm, &reqs[2*i]);
725 MPI_Issend(&sendBuffer[offsetInSharedSend[i]],
726 offsetInSharedSend[i+1] - offsetInSharedSend[i],
727 MPI_DIM_T, neighbour[i],
728 m_mpiInfo->counter()+m_mpiInfo->rank, m_mpiInfo->comm,
729 &reqs[2*i+1]);
730 }
731 m_mpiInfo->incCounter(m_mpiInfo->size);
732
733 // do something else here...
734
735 MPI_Waitall(2*neighbour.size(), &reqs[0], &stats[0]);
736
737 // now populate rest of node IDs
738 #pragma omp parallel for
739 for (index_t i=0; i < numRecv; i++) {
740 const index_t nodeIdx = recvShared[i];
741 m_nodeId[nodeIdx] = recvBuffer[m_dofMap[nodeIdx]-numDOF];
742 }
743 }
744 #endif // ESYS_MPI
745
746 // TODO: paso::SharedComponents should take vectors to avoid this
747 int* neighPtr = NULL;
748 index_t* sendPtr = NULL;
749 index_t* recvPtr = NULL;
750 if (neighbour.size() > 0) {
751 neighPtr = &neighbour[0];
752 sendPtr = &sendShared[0];
753 recvPtr = &recvShared[0];
754 }
755 // create Paso connector
756 paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
757 numDOF, neighbour.size(), neighPtr, sendPtr,
758 &offsetInSharedSend[0], 1, 0, m_mpiInfo));
759 paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
760 numDOF, neighbour.size(), neighPtr, recvPtr,
761 &offsetInSharedRecv[0], 1, 0, m_mpiInfo));
762 m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
763 }
764
765 void MultiRectangle::populateSampleIds()
766 {
767 // label nodes and DOF first
768 populateDofMap();
769
770 m_elementId.assign(getNumElements(), -1);
771
772 // populate face element counts
773 //left
774 if (m_offset[0]==0)
775 m_faceCount[0]=m_NE[1];
776 else
777 m_faceCount[0]=0;
778 //right
779 if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
780 m_faceCount[1]=m_NE[1];
781 else
782 m_faceCount[1]=0;
783 //bottom
784 if (m_offset[1]==0)
785 m_faceCount[2]=m_NE[0];
786 else
787 m_faceCount[2]=0;
788 //top
789 if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
790 m_faceCount[3]=m_NE[0];
791 else
792 m_faceCount[3]=0;
793
794 const dim_t NFE = getNumFaceElements();
795 const dim_t NE0 = m_NE[0];
796 const dim_t NE1 = m_NE[1];
797 m_faceId.resize(NFE);
798
799 #pragma omp parallel
800 {
801 // populate element IDs
802 #pragma omp for nowait
803 for (index_t i1 = 0; i1 < NE1; i1++) {
804 for (index_t i0 = 0; i0 < NE0; i0++) {
805 m_elementId[i0+i1*NE0]=(m_offset[1]+i1)*m_gNE[0]+m_offset[0]+i0;
806 }
807 }
808
809 // face elements
810 #pragma omp for
811 for (index_t k = 0; k < NFE; k++)
812 m_faceId[k] = k;
813 } // end parallel section
814
815 m_nodeTags.assign(getNumNodes(), 0);
816 updateTagsInUse(Nodes);
817
818 m_elementTags.assign(getNumElements(), 0);
819 updateTagsInUse(Elements);
820
821 // generate face offset vector and set face tags
822 const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
823 const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
824 m_faceOffset.assign(4, -1);
825 m_faceTags.clear();
826 index_t offset = 0;
827 for (size_t i=0; i<4; i++) {
828 if (m_faceCount[i] > 0) {
829 m_faceOffset[i]=offset;
830 offset += m_faceCount[i];
831 m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
832 }
833 }
834 setTagMap("left", LEFT);
835 setTagMap("right", RIGHT);
836 setTagMap("bottom", BOTTOM);
837 setTagMap("top", TOP);
838 updateTagsInUse(FaceElements);
839
840 }
841
842 paso::SystemMatrixPattern_ptr MultiRectangle::getPasoMatrixPattern(
843 bool reducedRowOrder,
844 bool reducedColOrder) const
845 {
846 if (!m_pattern) {
847 // first call - create pattern, then return
848 const dim_t numDOF = getNumDOF();
849 const dim_t numShared = getNumNodes() - numDOF;
850 #pragma omp parallel for
851 for (index_t i = 0; i < numShared; i++) {
852 sort(m_rowIndices[i].begin(), m_rowIndices[i].end());
853 }
854
855 // create main and couple blocks
856 paso::Pattern_ptr mainPattern = createPasoPattern(getConnections(), numDOF);
857 paso::Pattern_ptr colPattern = createPasoPattern(m_colIndices, numShared);
858 paso::Pattern_ptr rowPattern = createPasoPattern(m_rowIndices, numDOF);
859
860 // allocate Paso distribution
861 paso::Distribution_ptr distribution(new paso::Distribution(
862 m_mpiInfo, m_nodeDistribution, 1, 0));
863
864 // finally create the system matrix pattern
865 m_pattern.reset(new paso::SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
866 distribution, distribution, mainPattern, colPattern, rowPattern,
867 m_connector, m_connector));
868 }
869 return m_pattern;
870 }
871
872 RankVector MultiRectangle::getOwnerVector(int fsType) const
873 {
874 if (m_subdivisions != 1)
875 throw RipleyException("Multiresolution domains only support ownership for the coarsest level");
876 return Rectangle::getOwnerVector(fsType);
877 }
878
879 dim_t MultiRectangle::findNode(const double *coords) const
880 {
881 const dim_t NOT_MINE = -1;
882 //is the found element even owned by this rank
883 // (inside owned or shared elements but will map to an owned element)
884 for (int dim = 0; dim < m_numDim; dim++) {
885 double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
886 - m_dx[dim]/2.; //allows for point outside mapping onto node
887 double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
888 + m_dx[dim]/2.;
889 if (min > coords[dim] || max < coords[dim]) {
890 return NOT_MINE;
891 }
892 }
893 // get distance from origin
894 double x = coords[0] - m_origin[0];
895 double y = coords[1] - m_origin[1];
896
897 //check if the point is even inside the domain
898 if (x < 0 || y < 0 || x > m_length[0] || y > m_length[1])
899 return NOT_MINE;
900
901 // distance in elements
902 dim_t ex = (dim_t) floor((x + 0.01*m_dx[0]) / m_dx[0]);
903 dim_t ey = (dim_t) floor((y + 0.01*m_dx[1]) / m_dx[1]);
904 // set the min distance high enough to be outside the element plus a bit
905 dim_t closest = NOT_MINE;
906 double minDist = 1;
907 for (int dim = 0; dim < m_numDim; dim++) {
908 minDist += m_dx[dim]*m_dx[dim];
909 }
910 //find the closest node
911 for (int dx = 0; dx < 1; dx++) {
912 double xdist = (x - (ex + dx)*m_dx[0]);
913 for (int dy = 0; dy < 1; dy++) {
914 double ydist = (y - (ey + dy)*m_dx[1]);
915 double total = xdist*xdist + ydist*ydist;
916 if (total < minDist) {
917 closest = INDEX2(ex+dx-m_offset[0], ey+dy-m_offset[1], m_NN[0]);
918 minDist = total;
919 }
920 }
921 }
922 //if this happens, we've let a dirac point slip through, which is awful
923 if (closest == NOT_MINE) {
924 throw RipleyException("Unable to map appropriate dirac point to a node,"
925 " implementation problem in MultiRectangle::findNode()");
926 }
927 return closest;
928 }
929
930 } // end of namespace ripley
931

  ViewVC Help
Powered by ViewVC 1.1.26