/[escript]/trunk/speckley/src/Rectangle.cpp
ViewVC logotype

Contents of /trunk/speckley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5128 - (show annotations)
Fri Aug 29 00:52:35 2014 UTC (4 years, 7 months ago) by sshaw
File size: 23859 byte(s)
fixing a bad index calculation
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <algorithm>
18 #include <limits>
19
20 #include <speckley/Rectangle.h>
21 #include <esysUtils/esysFileWriter.h>
22 #include <esysUtils/index.h>
23 #include <speckley/DefaultAssembler2D.h>
24 #include <boost/scoped_array.hpp>
25 #include <escript/FunctionSpaceFactory.h>
26 #include "esysUtils/EsysRandom.h"
27
28 #ifdef USE_NETCDF
29 #include <netcdfcpp.h>
30 #endif
31
32 #if USE_SILO
33 #include <silo.h>
34 #ifdef ESYS_MPI
35 #include <pmpio.h>
36 #endif
37 #endif
38
39 #include <iomanip>
40
41 using esysUtils::FileWriter;
42
43 namespace speckley {
44
45 Rectangle::Rectangle(int order, dim_t n0, dim_t n1, double x0, double y0, double x1,
46 double y1, int d0, int d1,
47 const std::vector<double>& points,
48 const std::vector<int>& tags,
49 const simap_t& tagnamestonums,
50 escript::SubWorld_ptr w
51 ) :
52 SpeckleyDomain(2, order, w)
53 {
54 if (static_cast<long>(n0 + 1) * static_cast<long>(n1 + 1)
55 > std::numeric_limits<dim_t>::max())
56 throw SpeckleyException("The number of elements has overflowed, this "
57 "limit may be raised in future releases.");
58
59 if (n0 <= 0 || n1 <= 0)
60 throw SpeckleyException("Number of elements in each spatial dimension "
61 "must be positive");
62
63 // ignore subdivision parameters for serial run
64 if (m_mpiInfo->size == 1) {
65 d0=1;
66 d1=1;
67 } else {
68 throw SpeckleyException("MPI not currently supported by Rectangle");
69 }
70
71 bool warn=false;
72 std::vector<int> factors;
73 int ranks = m_mpiInfo->size;
74 dim_t epr[2] = {n0,n1};
75 int d[2] = {d0,d1};
76 if (d0<=0 || d1<=0) {
77 for (int i = 0; i < 2; i++) {
78 if (d[i] < 1) {
79 d[i] = 1;
80 continue;
81 }
82 epr[i] = -1; // can no longer be max
83 //remove
84 if (ranks % d[i] != 0) {
85 throw SpeckleyException("Invalid number of spatial subdivisions");
86 }
87 ranks /= d[i];
88 }
89 factorise(factors, ranks);
90 if (factors.size() != 0) {
91 warn = true;
92 }
93 }
94
95 while (factors.size() > 0) {
96 int i = epr[0] > epr[1] ? 0 : 1;
97 int f = factors.back();
98 factors.pop_back();
99 d[i] *= f;
100 epr[i] /= f;
101 }
102 d0 = d[0]; d1 = d[1];
103
104 // ensure number of subdivisions is valid and nodes can be distributed
105 // among number of ranks
106 if (d0*d1 != m_mpiInfo->size)
107 throw SpeckleyException("Invalid number of spatial subdivisions");
108
109 if (warn) {
110 std::cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
111 << d1 << "). This may not be optimal!" << std::endl;
112 }
113
114 double l0 = x1-x0;
115 double l1 = y1-y0;
116 m_dx[0] = l0/n0;
117 m_dx[1] = l1/n1;
118
119 if (n0 % d0 > 0) {
120 n0 += d0 - (n0 % d0);
121 l0 = m_dx[0]*n0;
122 std::cout << "Warning: Adjusted number of elements and length. N0="
123 << n0 << ", l0=" << l0 << std::endl;
124 }
125 if (n1 % d1 > 0) {
126 n1 += d1 - (n1 % d1);
127 l1 = m_dx[1]*n1;
128 std::cout << "Warning: Adjusted number of elements and length. N1="
129 << n1 << ", l1=" << l1 << std::endl;
130 }
131
132 if (n0/d0 < 2 || n1/d1 < 2)
133 throw SpeckleyException("Too few elements for the number of ranks");
134
135 m_gNE[0] = n0;
136 m_gNE[1] = n1;
137 m_origin[0] = x0;
138 m_origin[1] = y0;
139 m_length[0] = l0;
140 m_length[1] = l1;
141 m_NX[0] = d0;
142 m_NX[1] = d1;
143
144 // local number of elements
145 m_NE[0] = m_gNE[0] / d0;
146 m_NE[1] = m_gNE[1] / d1;
147
148 // local number of nodes
149 m_NN[0] = m_NE[0]*m_order+1;
150 m_NN[1] = m_NE[1]*m_order+1;
151
152 // bottom-left node is at (offset0,offset1) in global mesh
153 m_offset[0] = n0/d0*(m_mpiInfo->rank%d0);
154 m_offset[1] = n1/d1*(m_mpiInfo->rank/d0);
155
156 populateSampleIds();
157
158 for (simap_t::const_iterator i = tagnamestonums.begin();
159 i != tagnamestonums.end(); i++) {
160 setTagMap(i->first, i->second);
161 }
162 addPoints(points, tags);
163 }
164
165 Rectangle::~Rectangle()
166 {
167 }
168
169 std::string Rectangle::getDescription() const
170 {
171 return "speckley::Rectangle";
172 }
173
174 bool Rectangle::operator==(const AbstractDomain& other) const
175 {
176 const Rectangle* o=dynamic_cast<const Rectangle*>(&other);
177 if (o) {
178 return (SpeckleyDomain::operator==(other) &&
179 m_order == o->m_order &&
180 m_gNE[0]==o->m_gNE[0] && m_gNE[1]==o->m_gNE[1]
181 && m_origin[0]==o->m_origin[0] && m_origin[1]==o->m_origin[1]
182 && m_length[0]==o->m_length[0] && m_length[1]==o->m_length[1]
183 && m_NX[0]==o->m_NX[0] && m_NX[1]==o->m_NX[1]);
184 }
185
186 return false;
187 }
188
189 const dim_t* Rectangle::borrowSampleReferenceIDs(int fsType) const
190 {
191 switch (fsType) {
192 case DegreesOfFreedom:
193 case ReducedDegreesOfFreedom: // FIXME: reduced
194 // return &m_dofId[0];
195 case Nodes:
196 case ReducedNodes: // FIXME: reduced
197 return &m_nodeId[0];
198 case Elements:
199 case ReducedElements:
200 return &m_elementId[0];
201 case FaceElements:
202 case ReducedFaceElements:
203 return &m_faceId[0];
204 case Points:
205 return &m_diracPointNodeIDs[0];
206 default:
207 break;
208 }
209
210 std::stringstream msg;
211 msg << "borrowSampleReferenceIDs: invalid function space type " << fsType;
212 throw SpeckleyException(msg.str());
213 }
214
215 bool Rectangle::ownSample(int fsType, index_t id) const
216 {
217 throw SpeckleyException("ownSample not implemented");
218 }
219
220 void Rectangle::setToNormal(escript::Data& out) const
221 {
222 throw SpeckleyException("setToNormal not implemented");
223 }
224
225 void Rectangle::setToSize(escript::Data& out) const
226 {
227 if (out.getFunctionSpace().getTypeCode() == Elements) {
228 out.requireWrite();
229 const dim_t numQuad = out.getNumDataPointsPerSample();
230 const dim_t numElements = getNumElements();
231 const double size=sqrt(m_dx[0]*m_dx[0]+m_dx[1]*m_dx[1]);
232 #pragma omp parallel for
233 for (index_t k = 0; k < numElements; ++k) {
234 double* o = out.getSampleDataRW(k);
235 std::fill(o, o+numQuad, size);
236 }
237 } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
238 out.requireWrite();
239 const dim_t numQuad=out.getNumDataPointsPerSample();
240 const dim_t NE0 = m_NE[0];
241 const dim_t NE1 = m_NE[1];
242 #pragma omp parallel
243 {
244 if (m_faceOffset[0] > -1) {
245 #pragma omp for nowait
246 for (index_t k1 = 0; k1 < NE1; ++k1) {
247 double* o = out.getSampleDataRW(m_faceOffset[0]+k1);
248 std::fill(o, o+numQuad, m_dx[1]);
249 }
250 }
251
252 if (m_faceOffset[1] > -1) {
253 #pragma omp for nowait
254 for (index_t k1 = 0; k1 < NE1; ++k1) {
255 double* o = out.getSampleDataRW(m_faceOffset[1]+k1);
256 std::fill(o, o+numQuad, m_dx[1]);
257 }
258 }
259
260 if (m_faceOffset[2] > -1) {
261 #pragma omp for nowait
262 for (index_t k0 = 0; k0 < NE0; ++k0) {
263 double* o = out.getSampleDataRW(m_faceOffset[2]+k0);
264 std::fill(o, o+numQuad, m_dx[0]);
265 }
266 }
267
268 if (m_faceOffset[3] > -1) {
269 #pragma omp for nowait
270 for (index_t k0 = 0; k0 < NE0; ++k0) {
271 double* o = out.getSampleDataRW(m_faceOffset[3]+k0);
272 std::fill(o, o+numQuad, m_dx[0]);
273 }
274 }
275 } // end of parallel section
276
277 } else {
278 std::stringstream msg;
279 msg << "setToSize: invalid function space type "
280 << out.getFunctionSpace().getTypeCode();
281 throw SpeckleyException(msg.str());
282 }
283 }
284
285 void Rectangle::Print_Mesh_Info(const bool full) const
286 {
287 SpeckleyDomain::Print_Mesh_Info(full);
288 if (full) {
289 std::cout << " Id Coordinates" << std::endl;
290 std::cout.precision(15);
291 std::cout.setf(std::ios::scientific, std::ios::floatfield);
292 for (index_t i=0; i < getNumNodes(); i++) {
293 std::cout << " " << std::setw(5) << m_nodeId[i]
294 << " " << getLocalCoordinate(i%m_NN[0], 0)
295 << " " << getLocalCoordinate(i/m_NN[0], 1) << std::endl;
296 }
297 }
298 }
299
300
301 //protected
302 void Rectangle::assembleCoordinates(escript::Data& arg) const
303 {
304 escriptDataC x = arg.getDataC();
305 int numDim = m_numDim;
306 if (!isDataPointShapeEqual(&x, 1, &numDim))
307 throw SpeckleyException("setToX: Invalid Data object shape");
308 if (!numSamplesEqual(&x, 1, getNumNodes()))
309 throw SpeckleyException("setToX: Illegal number of samples in Data object");
310
311 const dim_t NN0 = m_NN[0];
312 const dim_t NN1 = m_NN[1];
313 arg.requireWrite();
314 #pragma omp parallel for
315 for (dim_t y = 0; y < NN1; y++) {
316 for (dim_t x = 0; x < NN0; x++) {
317 double *point = arg.getSampleDataRW(y*NN0 + x);
318 point[0] = getLocalCoordinate(x, 0);
319 point[1] = getLocalCoordinate(y, 1);
320 }
321 }
322 }
323
324 //protected
325 void Rectangle::assembleGradient(escript::Data& out, const escript::Data& in) const
326 {
327 escript::Data converted;
328
329 if (in.getFunctionSpace().getTypeCode() != Elements) {
330 converted = escript::Data(in, escript::function(*this));
331 } else {
332 converted = in;
333 }
334
335 if (m_order == 2) {
336 gradient_order2(out,converted);
337 } else if (m_order == 3) {
338 gradient_order3(out,converted);
339 } else if (m_order == 4) {
340 gradient_order4(out,converted);
341 } else if (m_order == 5) {
342 gradient_order5(out,converted);
343 } else if (m_order == 6) {
344 gradient_order6(out,converted);
345 } else if (m_order == 7) {
346 gradient_order7(out,converted);
347 } else if (m_order == 8) {
348 gradient_order8(out,converted);
349 } else if (m_order == 9) {
350 gradient_order9(out,converted);
351 } else if (m_order == 10) {
352 gradient_order10(out,converted);
353 }
354 }
355
356 //protected
357 void Rectangle::assembleIntegrate(std::vector<double>& integrals,
358 const escript::Data& arg) const
359 {
360 const int fs = arg.getFunctionSpace().getTypeCode();
361 if (fs != Elements)
362 throw new SpeckleyException("Speckley doesn't currently support integrals of non-Element functionspaces");
363 if (!arg.actsExpanded())
364 throw new SpeckleyException("Speckley doesn't currently support unexpanded data");
365
366 if (m_order == 2) {
367 integral_order2(integrals, arg);
368 } else if (m_order == 3) {
369 integral_order3(integrals, arg);
370 } else if (m_order == 4) {
371 integral_order4(integrals, arg);
372 } else if (m_order == 5) {
373 integral_order5(integrals, arg);
374 } else if (m_order == 6) {
375 integral_order6(integrals, arg);
376 } else if (m_order == 7) {
377 integral_order7(integrals, arg);
378 } else if (m_order == 8) {
379 integral_order8(integrals, arg);
380 } else if (m_order == 9) {
381 integral_order9(integrals, arg);
382 } else if (m_order == 10) {
383 integral_order10(integrals, arg);
384 }
385 }
386
387 //protected
388 dim_t Rectangle::insertNeighbourNodes(IndexVector& index, index_t node) const
389 {
390 std::cerr << "insertNeighbourNodes not updated\n";
391 const dim_t nDOF0 = (m_gNE[0]+1)/m_NX[0];
392 const dim_t nDOF1 = (m_gNE[1]+1)/m_NX[1];
393 const int x=node%nDOF0;
394 const int y=node/nDOF0;
395 dim_t num=0;
396 // loop through potential neighbours and add to index if positions are
397 // within bounds
398 for (int i1=-1; i1<2; i1++) {
399 for (int i0=-1; i0<2; i0++) {
400 // skip node itself
401 if (i0==0 && i1==0)
402 continue;
403 // location of neighbour node
404 const int nx=x+i0;
405 const int ny=y+i1;
406 if (nx>=0 && ny>=0 && nx<nDOF0 && ny<nDOF1) {
407 index.push_back(ny*nDOF0+nx);
408 num++;
409 }
410 }
411 }
412
413 return num;
414 }
415
416 /* This is a wrapper for filtered (and non-filtered) randoms
417 * For detailed doco see randomFillWorker
418 */
419 escript::Data Rectangle::randomFill(const escript::DataTypes::ShapeType& shape,
420 const escript::FunctionSpace& fs,
421 long seed, const boost::python::tuple& filter) const {
422 int numvals=escript::DataTypes::noValues(shape);
423 int per_element = (m_order+1)*(m_order+1)*numvals;
424 if (len(filter)>0) {
425 throw SpeckleyException("Speckley does not support filters.");
426 }
427
428 double* src=new double[m_NE[0]*m_NE[1]*per_element*numvals];
429 esysUtils::randomFillArray(seed, src, m_NE[0]*m_NE[1]*per_element);
430 escript::Data res(0, shape, escript::function(*this), true);
431 int current = 0;
432 for (int ei = 0; ei < m_NE[1]; ++ei) {
433 for (int ej = 0; ej < m_NE[0]; ++ej) {
434 double *e = res.getSampleDataRW(INDEX2(ej,ei,m_NE[0]));
435 memcpy(e, &src[current], sizeof(double)*per_element);
436 current += per_element;
437 }
438 }
439 if (res.getFunctionSpace() != fs) {
440 return escript::Data(res, fs);
441 }
442 return res;
443 }
444
445 //private
446 void Rectangle::populateSampleIds()
447 {
448 // degrees of freedom are numbered from left to right, bottom to top in
449 // each rank, continuing on the next rank (ranks also go left-right,
450 // bottom-top).
451 // This means rank 0 has id 0...n0-1, rank 1 has id n0...n1-1 etc. which
452 // helps when writing out data rank after rank.
453
454 // build node distribution vector first.
455 // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes which is
456 // constant for all ranks in this implementation
457 m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
458 const index_t left = (m_offset[0]==0 ? 0 : 1);
459 const index_t bottom = (m_offset[1]==0 ? 0 : 1);
460 for (dim_t k=1; k<m_mpiInfo->size; k++) {
461 index_t rank_left = (k-1)%m_NX[0] == 0 ? 0 : 1;
462 index_t rank_bottom = (k-1)/m_NX[0] == 0 ? 0 : 1;
463 m_nodeDistribution[k] = m_nodeDistribution[k-1]
464 + (m_NN[0]-rank_left)*(m_NN[1]-rank_bottom);
465 }
466 m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
467 try {
468 m_nodeId.resize(getNumNodes());
469 m_elementId.resize(getNumElements());
470 } catch (const std::length_error& le) {
471 throw SpeckleyException("The system does not have sufficient memory for a domain of this size.");
472 }
473
474 // populate face element counts
475 //left
476 if (m_offset[0]==0)
477 m_faceCount[0]=m_NE[1];
478 else
479 m_faceCount[0]=0;
480 //right
481 if (m_mpiInfo->rank%m_NX[0]==m_NX[0]-1)
482 m_faceCount[1]=m_NE[1];
483 else
484 m_faceCount[1]=0;
485 //bottom
486 if (m_offset[1]==0)
487 m_faceCount[2]=m_NE[0];
488 else
489 m_faceCount[2]=0;
490 //top
491 if (m_mpiInfo->rank/m_NX[0]==m_NX[1]-1)
492 m_faceCount[3]=m_NE[0];
493 else
494 m_faceCount[3]=0;
495
496 const dim_t NFE = getNumFaceElements();
497 m_faceId.resize(NFE);
498
499
500 if (bottom && left) {
501 //get lower-left node
502 m_nodeId[0] = m_nodeDistribution[m_mpiInfo->rank - m_NX[0]] - 1;
503 }
504 if (bottom) {
505 //DOF size of the neighbouring rank
506 index_t rankDOF = m_nodeDistribution[m_mpiInfo->rank - m_NX[0] + 1]
507 - m_nodeDistribution[m_mpiInfo->rank - m_NX[0]];
508 //beginning of last row of neighbouring rank
509 index_t begin = m_nodeDistribution[m_mpiInfo->rank - m_NX[0]]
510 + rankDOF - m_NN[0];
511 for (index_t x = left; x < m_NN[0]; x++) {
512 m_nodeId[x] = begin + x;
513 }
514 }
515 if (left) {
516 //is the rank to the left itself right of another rank
517 index_t rank_left = (m_mpiInfo->rank - 1) % m_NX[0] == 0 ? 0 : 1;
518 //end of first owned row of neighbouring rank
519 index_t end = m_nodeDistribution[m_mpiInfo->rank - 1]
520 + m_NN[0] - rank_left - 1;
521 for (index_t y = bottom; y < m_NN[1]; y++) {
522 m_nodeId[y*m_NN[0]] = end + (y-bottom)*(m_NN[0]-rank_left);
523 }
524 }
525
526 #pragma omp parallel for
527 for (index_t y = bottom; y < m_NN[1]; y++) {
528 for (index_t x = left; x < m_NN[0]; x++) {
529 m_nodeId[y*m_NN[0]+x] = m_nodeDistribution[m_mpiInfo->rank] + (y-bottom)*(m_NN[0]-left) + (x-left);
530 }
531 }
532
533 m_nodeTags.assign(getNumNodes(), 0);
534 updateTagsInUse(Nodes);
535
536 m_elementTags.assign(getNumElements(), 0);
537 updateTagsInUse(Elements);
538
539 // generate face offset vector and set face tags
540 const index_t LEFT=1, RIGHT=2, BOTTOM=10, TOP=20;
541 const index_t faceTag[] = { LEFT, RIGHT, BOTTOM, TOP };
542 m_faceOffset.assign(4, -1);
543 m_faceTags.clear();
544 index_t offset=0;
545 for (size_t i=0; i<4; i++) {
546 if (m_faceCount[i]>0) {
547 m_faceOffset[i]=offset;
548 offset+=m_faceCount[i];
549 m_faceTags.insert(m_faceTags.end(), m_faceCount[i], faceTag[i]);
550 }
551 }
552 setTagMap("left", LEFT);
553 setTagMap("right", RIGHT);
554 setTagMap("bottom", BOTTOM);
555 setTagMap("top", TOP);
556 updateTagsInUse(FaceElements);
557 }
558
559 //private
560 void Rectangle::addToMatrixAndRHS(escript::AbstractSystemMatrix* S, escript::Data& F,
561 const std::vector<double>& EM_S, const std::vector<double>& EM_F, bool addS,
562 bool addF, index_t firstNode, dim_t nEq, dim_t nComp) const
563 {
564 throw SpeckleyException("Rectangle::addToMatrixAndRHS, adding to matrix not supported");
565 }
566
567 //protected
568 void Rectangle::interpolateNodesOnElements(escript::Data& out,
569 const escript::Data& in,
570 bool reduced) const
571 {
572 if (reduced) {
573 throw SpeckleyException("Speckley domains do not support reduced function spaces");
574 }
575 const dim_t numComp = in.getDataPointSize();
576 const dim_t NE0 = m_NE[0];
577 const dim_t NE1 = m_NE[1];
578 const int quads = m_order + 1;
579 const int max_x = m_NN[0];
580 out.requireWrite();
581 #pragma omp parallel for
582 for (dim_t ey = 0; ey < NE1; ey++) {
583 for (dim_t ex = 0; ex < NE0; ex++) {
584 double *e_out = out.getSampleDataRW(ex + ey*NE0);
585 int start = ex*m_order + ey*max_x*m_order;
586 int quad = 0;
587 for (int qy = 0; qy < quads; qy++) {
588 for (int qx = 0; qx < quads; qx++, quad++) {
589 const double *n_in = in.getSampleDataRO(start + max_x*qy + qx);
590 memcpy(e_out+quad*numComp, n_in, sizeof(double) * numComp);
591 }
592 }
593 }
594 }
595 }
596
597 //protected
598 void Rectangle::interpolateElementsOnNodes(escript::Data& out,
599 const escript::Data& in, bool reduced) const {
600 const dim_t numComp = in.getDataPointSize();
601 const dim_t NE0 = m_NE[0];
602 const dim_t NE1 = m_NE[1];
603 const int quads = m_order + 1;
604 const dim_t max_x = (m_order*NE0) + 1;
605 const dim_t max_y = (m_order*NE1) + 1;
606 out.requireWrite();
607 if (reduced) {
608 throw SpeckleyException("Speckley domains do not support reduced function spaces");
609 }
610 // the summation portion
611 for (dim_t colouring = 0; colouring < 2; colouring++) {
612 #pragma omp parallel for
613 for (dim_t ey = colouring; ey < NE1; ey += 2) {
614 for (dim_t ex = 0; ex < NE0; ex++) {
615 int start = ex*m_order + ey*max_x*m_order;
616 const double *e_in = in.getSampleDataRO(ex + ey*NE0);
617 for (int qy = 0; qy < quads; qy++) {
618 for (int qx = 0; qx < quads; qx++) {
619 double *n_out = out.getSampleDataRW(start + max_x*qy + qx);
620 for (int comp = 0; comp < numComp; comp++) {
621 n_out[comp] += e_in[INDEX3(comp, qx, qy, numComp, quads)];
622 }
623 }
624 }
625 }
626 }
627 }
628 // the averaging out
629 // for every non-border edge in x
630 #pragma omp parallel for
631 for (dim_t qy = 0; qy < max_y; qy++) {
632 for (dim_t qx = m_order; qx < max_x - m_order; qx += m_order) {
633 double *n_out = out.getSampleDataRW(qx + qy*max_x);
634 for (int comp = 0; comp < numComp; comp++) {
635 n_out[comp] /= 2;
636 }
637 }
638 }
639
640 // for every non-border edge in y
641 #pragma omp parallel for
642 for (dim_t qy = m_order; qy < max_y - m_order; qy += m_order) {
643 for (dim_t qx = 0; qx < max_x; qx ++) {
644 double *n_out = out.getSampleDataRW(qx + qy*max_x);
645 for (int comp = 0; comp < numComp; comp++) {
646 n_out[comp] /= 2;
647 }
648 }
649 }
650 }
651
652
653 //protected
654 void Rectangle::interpolateNodesOnFaces(escript::Data& out,
655 const escript::Data& in,
656 bool reduced) const
657 {
658 throw SpeckleyException("Rectangle::interpolateNodesOnFaces not implemented");
659 }
660
661
662
663 dim_t Rectangle::findNode(const double *coords) const
664 {
665 const dim_t NOT_MINE = -1;
666 //is the found element even owned by this rank
667 for (int dim = 0; dim < m_numDim; dim++) {
668 double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
669 - m_dx[dim]/2.; //allows for point outside mapping onto node
670 double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
671 + m_dx[dim]/2.;
672 if (min > coords[dim] || max < coords[dim]) {
673 return NOT_MINE;
674 }
675 }
676 // get distance from origin
677 double x = coords[0] - m_origin[0];
678 double y = coords[1] - m_origin[1];
679
680 //check if the point is even inside the domain
681 if (x < 0 || y < 0 || x > m_length[0] || y > m_length[1])
682 return NOT_MINE;
683
684 // distance in elements
685 dim_t ex = (dim_t) floor(x / m_dx[0] + 0.01*m_dx[0]);
686 dim_t ey = (dim_t) floor(y / m_dx[1] + 0.01*m_dx[1]);
687 dim_t start = ex*m_order + ey*m_order*m_NN[0];
688 // set the min distance high enough to be outside the element plus a bit
689 dim_t closest = NOT_MINE;
690 double minDist = 1;
691 for (int dim = 0; dim < m_numDim; dim++) {
692 minDist += m_dx[dim]*m_dx[dim];
693 }
694 //find the closest node
695 for (int dx = 0; dx < 2; dx++) {
696 double xdist = x - (ex + dx)*m_dx[0];
697 for (int dy = 0; dy < 2; dy++) {
698 double ydist = y - (ey + dy)*m_dx[1];
699 double total = xdist*xdist + ydist*ydist;
700 if (total < minDist) {
701 closest = start + dx*m_order + dy*m_NN[0]*m_order;
702 std::cerr << "Rectangle::findNode not updated for MPI\n";
703 minDist = total;
704 }
705 }
706 }
707 //if this happens, we've let a dirac point slip through, which is awful
708 if (closest == NOT_MINE) {
709 throw SpeckleyException("Unable to map appropriate dirac point to a node,"
710 " implementation problem in Rectangle::findNode()");
711 }
712 return closest;
713 }
714
715 Assembler_ptr Rectangle::createAssembler(std::string type,
716 const DataMap& options) const {
717 if (type.compare("DefaultAssembler") == 0) {
718 return Assembler_ptr(new DefaultAssembler2D(shared_from_this(), m_dx, m_NX, m_NE, m_NN));
719 }
720 throw SpeckleyException("Speckley::Rectangle does not support the"
721 " requested assembler");
722 }
723
724 } // end of namespace speckley
725

  ViewVC Help
Powered by ViewVC 1.1.26