/[escript]/branches/trilinos_from_5897/ripley/src/MultiRectangle.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/ripley/src/MultiRectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6100 - (show annotations)
Tue Mar 29 04:56:51 2016 UTC (2 years, 10 months ago) by caltinay
File size: 38771 byte(s)
Fixed node IDs for 2D multiresolution ripley domains and added trilinos tests.

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

  ViewVC Help
Powered by ViewVC 1.1.26