/[escript]/trunk/bruce/src/Bruce.cpp
ViewVC logotype

Contents of /trunk/bruce/src/Bruce.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (show annotations)
Mon Mar 27 02:43:09 2006 UTC (14 years, 6 months ago) by robwdcock
File size: 35321 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


1 // $Id$
2 /*
3 ******************************************************************************
4 * *
5 * COPYRIGHT ACcESS 2005 - All Rights Reserved *
6 * *
7 * This software is the property of ACcESS. No part of this code *
8 * may be copied in any form or by any means without the expressed written *
9 * consent of ACcESS. Copying, use or modification of this software *
10 * by any unauthorised person is illegal unless that person has a software *
11 * license agreement with ACcESS. *
12 * *
13 ******************************************************************************
14 */
15
16 #include "Bruce.h"
17
18 #include "BruceException.h"
19
20 #include "escript/FunctionSpaceFactory.h"
21
22 #include <boost/python/extract.hpp>
23
24 #include "vtkCellType.h"
25
26 using namespace std;
27 using namespace escript;
28
29 namespace bruce {
30
31 const int Bruce::ContinuousFunction=0;
32 const int Bruce::Function=1;
33
34 Bruce::FunctionSpaceNamesMapType Bruce::m_functionSpaceTypeNames;
35
36 Bruce::Bruce():
37 m_n0(0), m_n1(0), m_n2(0)
38 {
39 setFunctionSpaceTypeNames();
40 }
41
42 Bruce::Bruce(DimVec v0, DimVec v1, DimVec v2,
43 int n0, int n1, int n2,
44 DimVec origin):
45 m_v0(v0), m_v1(v1), m_v2(v2),
46 m_n0(n0), m_n1(n1), m_n2(n2),
47 m_origin(origin)
48 {
49 if (!checkParameters()) {
50 stringstream temp;
51 temp << "Error - Invalid parameters supplied to constructor.";
52 throw BruceException(temp.str());
53 }
54 setFunctionSpaceTypeNames();
55 }
56
57 Bruce::Bruce(const Bruce& other):
58 m_v0(other.m_v0), m_v1(other.m_v1), m_v2(other.m_v2),
59 m_n0(other.m_n0), m_n1(other.m_n1), m_n2(other.m_n2),
60 m_origin(other.m_origin)
61 {
62 setFunctionSpaceTypeNames();
63 }
64
65 Bruce::~Bruce()
66 {
67 m_n0=-1;
68 m_n1=-1;
69 m_n2=-1;
70 m_origin.clear();
71 m_v0.clear();
72 m_v1.clear();
73 m_v2.clear();
74 }
75
76 void
77 Bruce::setFunctionSpaceTypeNames()
78 {
79 m_functionSpaceTypeNames.insert
80 (FunctionSpaceNamesMapType::value_type(Function,"Bruce_Function"));
81 m_functionSpaceTypeNames.insert
82 (FunctionSpaceNamesMapType::value_type(ContinuousFunction,"Bruce_ContinuousFunction"));
83 }
84
85 bool
86 Bruce::checkParameters()
87 {
88 if (m_origin.size()>3) {
89 return false;
90 }
91
92 if (m_n0<1) {
93 m_n0=1;
94 }
95 if (m_n1<1) {
96 m_n1=1;
97 }
98 if (m_n2<1) {
99 m_n2=1;
100 }
101
102 // reorder vectors and point counts according to point counts
103
104 //
105 // domains in 3d space
106 if (m_origin.size()==3) {
107
108 if (m_v0.size()==0) {
109
110 // 0d domain in 3d space
111 if ( (m_v1.size()!=0) || (m_v2.size()!=0) ) {
112 return false;
113 }
114
115 m_v0.clear();
116 m_v0.push_back(0);
117 m_v0.push_back(0);
118 m_v0.push_back(0);
119 m_v1.clear();
120 m_v1.push_back(0);
121 m_v1.push_back(0);
122 m_v1.push_back(0);
123 m_v2.clear();
124 m_v2.push_back(0);
125 m_v2.push_back(0);
126 m_v2.push_back(0);
127
128 m_n0=1;
129 m_n1=1;
130 m_n2=1;
131
132 } else if (m_v1.size()==0) {
133
134 // 1d domain in 3d space
135 if ( (m_v0.size()!=3) || (m_v2.size()!=0) ) {
136 return false;
137 }
138
139 m_v1.clear();
140 m_v1.push_back(0);
141 m_v1.push_back(0);
142 m_v1.push_back(0);
143 m_v2.clear();
144 m_v2.push_back(0);
145 m_v2.push_back(0);
146 m_v2.push_back(0);
147
148 m_n1=1;
149 m_n2=1;
150
151 } else if (m_v2.size()==0) {
152
153 // 2d domain in 3d space
154 if ( (m_v0.size()!=3) || (m_v1.size()!=3) ) {
155 return false;
156 }
157
158 m_v2.clear();
159 m_v2.push_back(0);
160 m_v2.push_back(0);
161 m_v2.push_back(0);
162
163 m_n2=1;
164
165 } else {
166
167 // 3d domain in 3d space
168 if ( (m_v0.size()!=3) || (m_v1.size()!=3) || (m_v2.size()!=3) ) {
169 return false;
170 }
171
172 }
173
174 }
175
176 //
177 // domains in 2d space
178 if (m_origin.size()==2) {
179
180 if (m_v0.size()==0) {
181
182 // 0d domain in 2d space
183 if (m_v1.size()!=0) {
184 return false;
185 }
186
187 m_v0.clear();
188 m_v0.push_back(0);
189 m_v0.push_back(0);
190 m_v1.clear();
191 m_v1.push_back(0);
192 m_v1.push_back(0);
193 m_v2.clear();
194
195 m_n0=1;
196 m_n1=1;
197 m_n2=1;
198
199 } else if (m_v1.size()==0) {
200
201 // 1d domain in 2d space
202 if (m_v0.size()!=2) {
203 return false;
204 }
205
206 m_v1.clear();
207 m_v1.push_back(0);
208 m_v1.push_back(0);
209 m_v2.clear();
210
211 m_n1=1;
212 m_n2=1;
213
214 } else {
215
216 // 2d domain in 2d space
217 if ( (m_v0.size()!=2) || (m_v1.size()!=2) ) {
218 return false;
219 }
220
221 m_v2.clear();
222
223 m_n2=1;
224
225 }
226
227 }
228
229 //
230 // domains in 1d space
231 if (m_origin.size()==1) {
232
233 if (m_v0.size()==0) {
234
235 // 0d domain in 1d space
236 m_v0.clear();
237 m_v0.push_back(0);
238 m_v1.clear();
239 m_v2.clear();
240
241 m_n0=1;
242 m_n1=1;
243 m_n2=1;
244
245 } else {
246
247 // 1d domain in 1d space
248 if (m_v0.size()!=1) {
249 return false;
250 }
251
252 m_v1.clear();
253 m_v2.clear();
254
255 m_n1=1;
256 m_n2=1;
257
258 }
259
260 }
261
262 //
263 // domains in 0d space
264 if (m_origin.size()==0) {
265
266 // point (0d) domain in 0d space
267 m_v0.clear();
268 m_v1.clear();
269 m_v2.clear();
270
271 m_n0=1;
272 m_n1=1;
273 m_n2=1;
274
275 }
276
277 return true;
278 }
279
280 bool
281 Bruce::isValidFunctionSpaceType(int functionSpaceCode) const
282 {
283 FunctionSpaceNamesMapType::iterator loc;
284 loc=m_functionSpaceTypeNames.find(functionSpaceCode);
285 return (loc!=m_functionSpaceTypeNames.end());
286 }
287
288 std::string
289 Bruce::functionSpaceTypeAsString(int functionSpaceCode) const
290 {
291 FunctionSpaceNamesMapType::iterator loc;
292 loc=m_functionSpaceTypeNames.find(functionSpaceCode);
293 if (loc==m_functionSpaceTypeNames.end()) {
294 stringstream temp;
295 temp << "Error - Invalid function space type code.";
296 throw BruceException(temp.str());
297 } else {
298 return loc->second;
299 }
300 }
301
302 pair<int,int>
303 Bruce::getDataShape(int functionSpaceCode) const
304 {
305 int numDataPointsPerSample=1;
306 int numSamples;
307
308 switch (functionSpaceCode) {
309
310 //
311 // Continuous functions
312 case(ContinuousFunction):
313
314 numSamples = m_n0 * m_n1 * m_n2;
315 break;
316
317 //
318 // Functions
319 case(Function):
320
321 // 0d domains
322 if (getDim()==0) {
323
324 numSamples = 0;
325
326 // 1d domains
327 } else if (getDim()==1) {
328
329 if (isZero(m_v0)) {
330 numSamples = 0;
331 } else {
332 numSamples = m_n0-1;
333 }
334
335 // 2d domains
336 } else if (getDim()==2) {
337
338 if (isZero(m_v0)) {
339 numSamples = 0;
340 } else if (isZero(m_v1)) {
341 numSamples = m_n0-1;
342 } else {
343 numSamples = (m_n0-1) * (m_n1-1);
344 }
345
346 // 3d domains
347 } else {
348
349 if (isZero(m_v0)) {
350 numSamples = 0;
351 } else if (isZero(m_v1)) {
352 numSamples = m_n0-1;
353 } else if (isZero(m_v2)) {
354 numSamples = (m_n0-1) * (m_n1-1);
355 } else {
356 numSamples = (m_n0-1) * (m_n1-1) * (m_n2-1);
357 }
358
359 }
360
361 break;
362
363 default:
364 stringstream temp;
365 temp << "Error - Invalid function space type: "
366 << functionSpaceCode << " for domain: " << getDescription();
367 throw BruceException(temp.str());
368 break;
369
370 }
371
372 return pair<int,int>(numDataPointsPerSample,numSamples);
373 }
374
375 int
376 Bruce::getNumSamples(int functionSpaceCode) const
377 {
378 std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
379 return domainShape.second;
380 }
381
382 Data
383 Bruce::getX() const
384 {
385 FunctionSpace tempFunc = continuousFunction(asAbstractContinuousDomain());
386 return tempFunc.getX();
387 }
388
389 void
390 Bruce::setToX(escript::Data& out) const
391 {
392
393 //
394 // determine functionSpace type expected by supplied Data object
395 int functionSpaceCode = out.getFunctionSpace().getTypeCode();
396
397 //
398 // ensure supplied Data object has a matching number of data-points
399 // for this Bruce domain object
400 std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
401 if(domainShape.first!=out.getNumDataPointsPerSample() ||
402 domainShape.second!=out.getNumSamples()) {
403 stringstream temp;
404 temp << "Error - Incompatible number of data-points Data object supplied to Bruce::setToX";
405 throw BruceException(temp.str());
406 }
407
408 int dim = getDim();
409 int numSamples = domainShape.second;
410
411 //
412 // ensure shape of data-points in supplied Data object matches the
413 // shape needed to store the coordinates of each data-point in this
414 // Bruce domain
415 std::vector<int> dataShape = out.getDataPointShape();
416 if (dim>0 && (dataShape.size()!=1 || dataShape[0]!=dim)) {
417 stringstream temp;
418 temp << "Error - Incompatible shape Data object supplied to Bruce::setToX";
419 throw BruceException(temp.str());
420 } else if (dim==0 && dataShape.size()!=0) {
421 stringstream temp;
422 temp << "Error - Incompatible shape Data object supplied to Bruce::setToX";
423 throw BruceException(temp.str());
424 }
425
426 if (functionSpaceCode==ContinuousFunction) {
427
428 //
429 // calculate the spatial coordinates of data-points
430 // located on the nodes of this Bruce domain
431
432 if (dim==0) {
433
434 // Bruce domains in 0d space
435
436 int sampleNo=0;
437 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
438
439 } else if (dim==1) {
440
441 // Bruce domains in 1d space
442
443 for (int i=0; i<m_n0; i++) {
444 int sampleNo=i;
445 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
446 sampleData[0] = m_origin[0] + m_v0[0]*i;
447 }
448
449 } else if (dim==2) {
450
451 // Bruce domains in 2d space
452
453 for (int i=0; i<m_n0; i++) {
454 for (int j=0; j<m_n1; j++) {
455 int sampleNo=(m_n1*i)+j;
456 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
457 for (int d=0; d<dim; d++) {
458 sampleData[d] = m_origin[d] + m_v0[d]*i + m_v1[d]*j;
459 }
460 }
461 }
462
463 } else if (dim==3) {
464
465 // Bruce domains in 3d space
466
467 for (int i=0; i<m_n0; i++) {
468 for (int j=0; j<m_n1; j++) {
469 for (int k=0; k<m_n2; k++) {
470 int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
471 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
472 for (int d=0; d<dim; d++) {
473 sampleData[d] = m_origin[d] + m_v0[d]*i + m_v1[d]*j + m_v2[d]*k;
474 }
475 }
476 }
477 }
478
479 }
480
481 } else if (functionSpaceCode==Function) {
482
483 //
484 // calculate the spatial coordinates of data-points
485 // located on the cell centres of this Bruce domain
486
487 if (dim==0) {
488
489 // Bruce domains in 0d space
490
491 stringstream temp;
492 temp << "Error - Invalid function space type: "
493 << functionSpaceCode << " for Bruce::setToX";
494 throw BruceException(temp.str());
495
496 } else if (dim==1) {
497
498 // Bruce domains in 1d space
499
500 int n0max=m_n0-1;
501 if (isZero(m_v0)) {
502 stringstream temp;
503 temp << "Error - Invalid function space type: "
504 << functionSpaceCode << " for Bruce::setToX";
505 throw BruceException(temp.str());
506 } else {
507 for (int i=0; i<n0max; i++) {
508 int sampleNo=i;
509 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
510 sampleData[0] = m_origin[0] + m_v0[0]*(i + 0.5);
511 }
512 }
513
514 } else if (dim==2) {
515
516 // Bruce domains in 2d space
517
518 int n0max=m_n0-1;
519 int n1max=m_n1-1;
520 if (isZero(m_v0)) {
521 stringstream temp;
522 temp << "Error - Invalid function space type: "
523 << functionSpaceCode << " for Bruce::setToX";
524 throw BruceException(temp.str());
525 } else if (isZero(m_v1)) {
526 for (int i=0; i<n0max; i++) {
527 int sampleNo=i;
528 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
529 for (int d=0; d<dim; d++) {
530 sampleData[d] = m_origin[d] + m_v0[d]*(i + 0.5);
531 }
532 }
533 } else {
534 for (int i=0; i<n0max; i++) {
535 for (int j=0; j<n1max; j++) {
536 int sampleNo=(n1max*i)+j;
537 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
538 for (int d=0; d<dim; d++) {
539 sampleData[d] = m_origin[d] + m_v0[d]*(i + 0.5) + m_v1[d]*(j + 0.5);
540 }
541 }
542 }
543 }
544
545 } else if (dim==3) {
546
547 // Bruce domains in 3d space
548
549 int n0max=m_n0-1;
550 int n1max=m_n1-1;
551 int n2max=m_n2-1;
552 if (isZero(m_v0)) {
553 stringstream temp;
554 temp << "Error - Invalid function space type: "
555 << functionSpaceCode << " for Bruce::setToX";
556 throw BruceException(temp.str());
557 } else if (isZero(m_v1)) {
558 for (int i=0; i<n0max; i++) {
559 int sampleNo=i;
560 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
561 for (int d=0; d<dim; d++) {
562 sampleData[d] = m_origin[d] + m_v0[d]*(i + 0.5);
563 }
564 }
565 } else if (isZero(m_v2)) {
566 for (int i=0; i<n0max; i++) {
567 for (int j=0; j<n1max; j++) {
568 int sampleNo=(n1max*i)+j;
569 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
570 for (int d=0; d<dim; d++) {
571 sampleData[d] = m_origin[d] + m_v0[d]*(i + 0.5) + m_v1[d]*(j + 0.5);
572 }
573 }
574 }
575 } else {
576 for (int i=0; i<n0max; i++) {
577 for (int j=0; j<n1max; j++) {
578 for (int k=0; k<n2max; k++) {
579 int sampleNo=(n1max*n2max*i)+(n2max*j)+k;
580 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
581 for (int d=0; d<dim; d++) {
582 sampleData[d] = m_origin[d] + m_v0[d]*(i + 0.5) + m_v1[d]*(j + 0.5) + m_v2[d]*(k + 0.5);
583 }
584 }
585 }
586 }
587 }
588
589 }
590
591 } else {
592 stringstream temp;
593 temp << "Error - Invalid function space type: "
594 << functionSpaceCode << " for domain: " << getDescription();
595 throw BruceException(temp.str());
596 }
597
598 }
599
600 Data
601 Bruce::getSize() const
602 {
603 FunctionSpace tempFunc = function(asAbstractContinuousDomain());
604 return tempFunc.getSize();
605 }
606
607 void
608 Bruce::setToSize(escript::Data& out) const
609 {
610
611 //
612 // determine functionSpace type expected by supplied Data object
613 int functionSpaceCode = out.getFunctionSpace().getTypeCode();
614
615 //
616 // ensure supplied Data object has a matching number of data-points
617 // for this Bruce domain object
618 std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
619 if(domainShape.first!=out.getNumDataPointsPerSample() ||
620 domainShape.second!=out.getNumSamples()) {
621 stringstream temp;
622 temp << "Error - Incompatible number of data-points Data object supplied to Bruce::setToSize";
623 throw BruceException(temp.str());
624 }
625
626 int dim = getDim();
627 int numSamples = domainShape.second;
628
629 //
630 // ensure shape of data-points in supplied Data object matches the
631 // shape needed to store the size of each data-point in this Bruce domain
632 std::vector<int> dataShape = out.getDataPointShape();
633 // this check should be satisfied by Data objects passed to setToSize, but
634 // FunctionSpace::getSize() seems to create an object which is larger than
635 // this... either way, this method can deal with this
636 /*
637 if (dataShape.size()!=1 || dataShape[0]!=1) {
638 stringstream temp;
639 temp << "Error - Incompatible shape Data object supplied to Bruce::setToSize";
640 throw BruceException(temp.str());
641 }
642 */
643
644 double dp_size;
645
646 if (functionSpaceCode==ContinuousFunction) {
647
648 //
649 // calculate the size of data-points
650 // located on the nodes of this Bruce domain
651
652 if (dim==0) {
653
654 // Bruce domains in 0d space
655
656 dp_size = 0.0;
657
658 int sampleNo=0;
659 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
660 sampleData[0] = dp_size;
661
662 } else if (dim==1) {
663
664 // Bruce domains in 1d space
665
666 dp_size = m_v0[0];
667
668 for (int i=0; i<m_n0; i++) {
669 int sampleNo=i;
670 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
671 sampleData[0] = dp_size;
672 }
673
674 } else if (dim==2) {
675
676 // Bruce domains in 2d space
677
678 double x = m_v0[0] + m_v1[0];
679 double y = m_v0[1] + m_v1[1];
680 dp_size = sqrt(pow(x,2)+pow(y,2));
681
682 for (int i=0; i<m_n0; i++) {
683 for (int j=0; j<m_n1; j++) {
684 int sampleNo=(m_n1*i)+j;
685 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
686 sampleData[0] = dp_size;
687 }
688 }
689
690 } else if (dim==3) {
691
692 // Bruce domains in 3d space
693
694 double x = m_v0[0] + m_v1[0] + m_v2[0];
695 double y = m_v0[1] + m_v1[1] + m_v2[1];
696 double z = m_v0[2] + m_v1[2] + m_v2[2];
697 dp_size = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
698
699 for (int i=0; i<m_n0; i++) {
700 for (int j=0; j<m_n1; j++) {
701 for (int k=0; k<m_n2; k++) {
702 int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
703 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
704 sampleData[0] = dp_size;
705 }
706 }
707 }
708
709 }
710
711 } else if (functionSpaceCode==Function) {
712
713 //
714 // calculate the sizes of data-points
715 // located on the cell centres of this Bruce domain
716
717 if (dim==0) {
718
719 // Bruce domains in 0d space
720
721 stringstream temp;
722 temp << "Error - Invalid function space type: "
723 << functionSpaceCode << " for Bruce::setToSize";
724 throw BruceException(temp.str());
725
726 } else if (dim==1) {
727
728 // Bruce domains in 1d space
729
730 dp_size = m_v0[0];
731
732 int n0max=m_n0-1;
733 if (isZero(m_v0)) {
734 stringstream temp;
735 temp << "Error - Invalid function space type: "
736 << functionSpaceCode << " for Bruce::setToSize";
737 throw BruceException(temp.str());
738 } else {
739 for (int i=0; i<n0max; i++) {
740 int sampleNo=i;
741 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
742 sampleData[0] = dp_size;
743 }
744 }
745
746 } else if (dim==2) {
747
748 // Bruce domains in 2d space
749
750 double x = m_v0[0] + m_v1[0];
751 double y = m_v0[1] + m_v1[1];
752 dp_size = sqrt(pow(x,2)+pow(y,2));
753
754 int n0max=m_n0-1;
755 int n1max=m_n1-1;
756 if (isZero(m_v0)) {
757 stringstream temp;
758 temp << "Error - Invalid function space type: "
759 << functionSpaceCode << " for Bruce::setToSize";
760 throw BruceException(temp.str());
761 } else if (isZero(m_v1)) {
762 for (int i=0; i<n0max; i++) {
763 int sampleNo=i;
764 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
765 sampleData[0] = dp_size;
766 }
767 } else {
768 for (int i=0; i<n0max; i++) {
769 for (int j=0; j<n1max; j++) {
770 int sampleNo=(n1max*i)+j;
771 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
772 sampleData[0] = dp_size;
773 }
774 }
775 }
776
777 } else if (dim==3) {
778
779 // Bruce domains in 3d space
780
781 double x = m_v0[0] + m_v1[0] + m_v2[0];
782 double y = m_v0[1] + m_v1[1] + m_v2[1];
783 double z = m_v0[2] + m_v1[2] + m_v2[2];
784 dp_size = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
785
786 int n0max=m_n0-1;
787 int n1max=m_n1-1;
788 int n2max=m_n2-1;
789 if (isZero(m_v0)) {
790 stringstream temp;
791 temp << "Error - Invalid function space type: "
792 << functionSpaceCode << " for Bruce::setToSize";
793 throw BruceException(temp.str());
794 } else if (isZero(m_v1)) {
795 for (int i=0; i<n0max; i++) {
796 int sampleNo=i;
797 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
798 sampleData[0] = dp_size;
799 }
800 } else if (isZero(m_v2)) {
801 for (int i=0; i<n0max; i++) {
802 for (int j=0; j<n1max; j++) {
803 int sampleNo=(n1max*i)+j;
804 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
805 sampleData[0] = dp_size;
806 }
807 }
808 } else {
809 for (int i=0; i<n0max; i++) {
810 for (int j=0; j<n1max; j++) {
811 for (int k=0; k<n2max; k++) {
812 int sampleNo=(n2max*n1max*i)+(n1max*j)+k;
813 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
814 sampleData[0] = dp_size;
815 }
816 }
817 }
818 }
819
820 }
821
822 } else {
823 stringstream temp;
824 temp << "Error - Invalid function space type: "
825 << functionSpaceCode << " for domain: " << getDescription();
826 throw BruceException(temp.str());
827 }
828
829 }
830
831 void
832 Bruce::setToGradient(escript::Data& grad,
833 const escript::Data& arg) const
834 {
835 // pre-conditions: arg must be rank 0->3 only, rank 4 cannot be accomodated
836 // grad will be rank of arg, +1
837 // grad is calculated by, for each data-point in this Bruce domain, retrieving
838 // the associated value from arg, calculating the grad, and loading this value
839 // into the corresponding data-point in grad
840 }
841
842 bool
843 Bruce::operator==(const AbstractDomain& other) const
844 {
845 const Bruce* temp=dynamic_cast<const Bruce*>(&other);
846 if (temp!=0) {
847 if ((m_v0 != temp->m_v0) || (m_v1 != temp->m_v1) || (m_v2 != temp->m_v2)) {
848 return false;
849 }
850 if ((m_n0 != temp->m_n0) || (m_n1 != temp->m_n1) || (m_n2 != temp->m_n2)) {
851 return false;
852 }
853 if (m_origin != temp->m_origin) {
854 return false;
855 }
856 return true;
857 } else {
858 return false;
859 }
860 }
861
862 bool
863 Bruce::operator!=(const AbstractDomain& other) const
864 {
865 return !(operator==(other));
866 }
867
868 int
869 Bruce::getReferenceNoFromSampleNo(int functionSpaceCode,
870 int sampleNo) const
871 {
872 // ensure supplied sampleNo is valid
873 std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
874 int numSamples = domainShape.second;
875 if ( (sampleNo>=numSamples) || (sampleNo < 0)) {
876 stringstream temp;
877 temp << "Bruce::getReferenceNoFromSampleNo: Error - invalid sampleNo supplied";
878 throw BruceException(temp.str());
879 }
880 // we just set the reference number to be the sample number for all samples
881 return sampleNo;
882 }
883
884 void
885 Bruce::saveVTK(const std::string& filename,
886 const boost::python::dict& dataDict) const
887 {
888
889 //
890 // extract Data objects and associated names from dictionary object
891 const int num_data=boost::python::extract<int>(dataDict.attr("__len__")());
892
893 int MAX_namelength=256;
894 char names[num_data][MAX_namelength];
895 char* c_names[num_data];
896
897 escriptDataC data[num_data];
898 escriptDataC* ptr_data[num_data];
899
900 boost::python::list keys=dataDict.keys();
901 for (int i=0; i<num_data; i++) {
902 Data& d=boost::python::extract<Data&>(dataDict[keys[i]]);
903 if (dynamic_cast<const Bruce&>(d.getFunctionSpace().getDomain()) != *this) {
904 throw BruceException("Error: Bruce::saveVTK: Data must be defined on same Domain");
905 }
906 data[i]=d.getDataC();
907 ptr_data[i]=&(data[i]);
908 std::string n=boost::python::extract<std::string>(keys[i]);
909 c_names[i]=&(names[i][0]);
910 if (n.length()>MAX_namelength-1) {
911 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
912 c_names[i][MAX_namelength-1]='\0';
913 } else {
914 strcpy(c_names[i],n.c_str());
915 }
916 }
917
918 //
919 // open archive file
920 ofstream archiveFile;
921 archiveFile.open(filename.data(), ios::out);
922 if (!archiveFile.good()) {
923 throw BruceException("Error in Bruce::saveVTK: File could not be opened for writing");
924 }
925
926 //
927 // determine mesh type to be written
928 bool isCellCentered[num_data], write_celldata=false, write_pointdata=false;
929 for (int i_data=0; i_data<num_data; i_data++) {
930 if (!isEmpty(ptr_data[i_data])) {
931 switch(getFunctionSpaceType(ptr_data[i_data])) {
932 case ContinuousFunction:
933 isCellCentered[i_data]=false;
934 break;
935 case Function:
936 isCellCentered[i_data]=true;
937 break;
938 }
939 if (isCellCentered[i_data]) {
940 write_celldata=true;
941 } else {
942 write_pointdata=true;
943 }
944 }
945 }
946
947 //
948 // determine number of points and cells
949 int numPoints=getDataShape(ContinuousFunction).second;
950 int numCells=getDataShape(Function).second;
951
952 //
953 // determine VTK element type
954 int cellType;
955 char elemTypeStr[32];
956 int numVTKNodesPerElement;
957
958 int nDim = getDim();
959 switch (nDim) {
960
961 case 0:
962 cellType = VTK_VERTEX;
963 numVTKNodesPerElement = 1;
964 strcpy(elemTypeStr, "VTK_VERTEX");
965 break;
966
967 case 1:
968 cellType = VTK_LINE;
969 numVTKNodesPerElement = 2;
970 strcpy(elemTypeStr, "VTK_LINE");
971 break;
972
973 case 2:
974 cellType = VTK_QUAD;
975 numVTKNodesPerElement = 4;
976 strcpy(elemTypeStr, "VTK_QUAD");
977 break;
978
979 case 3:
980 cellType = VTK_HEXAHEDRON;
981 numVTKNodesPerElement = 8;
982 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
983 break;
984
985 }
986
987 //
988 // write xml header
989 archiveFile << "<?xml version=\"1.0\"?>" << endl;
990
991 //
992 // determine grid extent
993
994 // ??
995
996 //
997 // write grid type and extent
998 archiveFile << "\t<VTKFile type=\"StructuredGrid\" version=\"0.1\">" << endl;
999 archiveFile << "\t\t<StructuredGrid WholeExtent=\"x1 x2 y1 y2 z1 z2\">" << endl;
1000 archiveFile << "\t\t\t<Piece Extent=\"x1 x2 y1 y2 z1 z2\">" << endl;
1001
1002 //
1003 // start to write out point definitions
1004 archiveFile << "\t\t\t\t<Points>" << endl;
1005
1006 //
1007 // determine grid cooordinates
1008
1009 // ??
1010
1011 // vtk/mayavi doesn't like 2D data, it likes 3D data with
1012 // a degenerate third dimension to handle 2D data.
1013 // So if nDim is less than 3, must pad all empty dimensions,
1014 // so that the total number of dims is 3.
1015
1016 archiveFile << "\t\t\t\t\t<DataArray NumberOfComponents=" << max(3,nDim) << " type=\"Float32\" format=\"ascii\">" << endl;
1017 /*
1018 for (int i=0; i<numPoints; i++) {
1019 for (int j=0; j<nDim; j++)
1020 archiveFile << "%e "; //mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)])
1021 for (int k=0; !(k>=3-nDim); k++)
1022 archiveFile << "0. ";
1023 archiveFile << endl;
1024 }
1025 */
1026 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1027 archiveFile << "\t\t\t\t</Points>" << endl;
1028
1029 //
1030 // cell data
1031 if (write_celldata) {
1032
1033 //
1034 // mark the active cell-data data arrays
1035 bool set_scalar=false, set_vector=false, set_tensor=false;
1036 archiveFile << "\t\t\t\t<CellData";
1037 for (int i_data=0; i_data<num_data; i_data++) {
1038 if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1039
1040 switch(getDataPointRank(ptr_data[i_data])) {
1041
1042 case 0:
1043 if (!set_scalar) {
1044 archiveFile << " Scalars=" << names[i_data];
1045 set_scalar=true;
1046 }
1047 break;
1048
1049 case 1:
1050 if (!set_vector) {
1051 archiveFile << " Vectors=" << names[i_data];
1052 set_vector=true;
1053 }
1054 break;
1055
1056 case 2:
1057 if (!set_tensor) {
1058 archiveFile << " Tensors=" << names[i_data];
1059 set_tensor=true;
1060 }
1061 break;
1062
1063 default:
1064 archiveFile.close();
1065 throw BruceException("saveVTK: VTK can't handle objects with rank greater than 2.");
1066
1067 }
1068 }
1069 }
1070 archiveFile << ">" << endl;
1071
1072 //
1073 // write the cell-data data arrays
1074 for (int i_data=0; i_data<num_data; i_data++) {
1075 if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1076
1077 int numPointsPerSample=1;
1078 int rank=getDataPointRank(ptr_data[i_data]);
1079 int nComp=getDataPointSize(ptr_data[i_data]);
1080 int nCompReqd; // the number of components required by vtk
1081 int shape=0;
1082
1083 switch (rank) {
1084
1085 case 0:
1086 nCompReqd=1;
1087 break;
1088
1089 case 1:
1090 shape=getDataPointShape(ptr_data[i_data], 0);
1091 if (shape>3) {
1092 archiveFile.close();
1093 throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1094 }
1095 nCompReqd=3;
1096 break;
1097
1098 case 2:
1099 shape=getDataPointShape(ptr_data[i_data], 0);
1100 if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1101 archiveFile.close();
1102 throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1103 }
1104 nCompReqd=9;
1105 break;
1106
1107 }
1108
1109 archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1110
1111 /*
1112 //
1113 // write out the cell data
1114
1115 // if the number of required components is more than the number
1116 // of actual components, pad with zeros
1117 // probably only need to get shape of first element
1118 // write the data different ways for scalar, vector and tensor
1119
1120 for (int i=0; i<numCells; i++) {
1121
1122 double *values=getSampleData(ptr_data[i_data], i);
1123 double sampleAvg[nComp];
1124
1125 // average over the number of points in the sample (ie: 1)
1126 for (int k=0; k<nComp; k++) {
1127 sampleAvg[k] = values[INDEX2(k,0,nComp)];
1128 }
1129
1130 if (nCompReqd == 1) {
1131
1132 archiveFile << sampleAvg[0] << endl;
1133
1134 } else if (nCompReqd == 3) {
1135
1136 for (int m=0; m<shape; m++) {
1137 archiveFile << sampleAvg[m] << endl;
1138 }
1139 for (int m=0; m<nCompReqd-shape; m++) {
1140 archiveFile << "0." << endl;
1141 }
1142
1143 } else if (nCompReqd == 9) {
1144
1145 int count=0;
1146 for (int m=0; m<shape; m++) {
1147 for (int n=0; n<shape; n++) {
1148 archiveFile << sampleAvg[count] << endl;
1149 count++;
1150 }
1151 for (int n=0; n<3-shape; n++) {
1152 archiveFile << "0." << endl;
1153 }
1154 }
1155 for (int m=0; m<3-shape; m++) {
1156 for (int n=0; n<3; n++) {
1157 archiveFile << "0." << endl;
1158 }
1159 }
1160
1161 }
1162 }
1163 */
1164 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1165
1166 }
1167 }
1168
1169 archiveFile << "\t\t\t\t</CellData>" << endl;
1170 }
1171
1172 //
1173 // point data
1174 if (write_pointdata) {
1175
1176 //
1177 // mark the active point-data data arrays
1178 bool set_scalar=false, set_vector=false, set_tensor=false;
1179 archiveFile << "\t\t\t\t<PointData";
1180 for (int i_data=0; i_data<num_data; i_data++) {
1181 if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1182
1183 switch(getDataPointRank(ptr_data[i_data])) {
1184
1185 case 0:
1186 if (!set_scalar) {
1187 archiveFile << " Scalars=" << names[i_data];
1188 set_scalar=true;
1189 }
1190 break;
1191
1192 case 1:
1193 if (!set_vector) {
1194 archiveFile << " Vectors=" << names[i_data];
1195 set_vector=true;
1196 }
1197 break;
1198
1199 case 2:
1200 if (!set_tensor) {
1201 archiveFile << " Tensors=" << names[i_data];
1202 set_tensor=true;
1203 }
1204 break;
1205
1206 default:
1207 archiveFile.close();
1208 throw BruceException("saveVTK: Vtk can't handle objects with rank greater than 2");
1209 }
1210
1211 }
1212 }
1213 archiveFile << ">" << endl;
1214
1215 //
1216 // write the point-data data arrays
1217 for (int i_data=0; i_data<num_data; i_data++) {
1218 if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1219
1220 int numPointsPerSample=1;
1221 int rank=getDataPointRank(ptr_data[i_data]);
1222 int nComp=getDataPointSize(ptr_data[i_data]);
1223 int nCompReqd; // the number of components required by vtk
1224 int shape=0;
1225
1226 switch (rank) {
1227
1228 case 0:
1229 nCompReqd=1;
1230
1231 case 1:
1232 shape=getDataPointShape(ptr_data[i_data], 0);
1233 if (shape>3) {
1234 archiveFile.close();
1235 throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1236 }
1237 nCompReqd=3;
1238
1239 case 2:
1240 shape=getDataPointShape(ptr_data[i_data], 0);
1241 if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1242 archiveFile.close();
1243 throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1244 }
1245 nCompReqd=9;
1246
1247 }
1248
1249 archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1250
1251 /*
1252 //
1253 // write out the point data
1254
1255 // if the number of required components is more than
1256 // the number of actual components, pad with zeros
1257 // write the data different ways for scalar, vector and tensor
1258
1259 for (int i=0; i<numNodes; i++) {
1260
1261 values=getSampleData(data[i_data], i);
1262
1263 if (nCompReqd==1) {
1264
1265 archiveFile << values[0];
1266
1267 } else if (nCompReqd==3) {
1268
1269 for (int m=0; m<shape; m++) {
1270 archiveFile << values[m];
1271 }
1272 for (int m=0; m<nCompReqd-shape; m++) {
1273 archiveFile << "0.";
1274 }
1275
1276 } else if (nCompReqd==9) {
1277
1278 int count=0;
1279 for (int m=0; m<shape; m++) {
1280 for (int n=0; n<shape; n++) {
1281 archiveFile << values[count];
1282 count++;
1283 }
1284 for (int n=0; n<3-shape; n++) {
1285 archiveFile << "0.";
1286 }
1287 }
1288 for (int m=0; m<3-shape; m++) {
1289 for (int n=0; n<3; n++) {
1290 archiveFile << "0.";
1291 }
1292 }
1293
1294 }
1295
1296 archiveFile << endl;
1297 }
1298 */
1299
1300 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1301 }
1302 }
1303
1304 archiveFile << "\t\t\t\t</PointData>" << endl;
1305 }
1306
1307 //
1308 // finish off the grid definition
1309 archiveFile << "\t\t\t</Piece>" << endl;
1310 archiveFile << "\t\t</StructuredGrid>" << endl;
1311
1312 //
1313 // write the xml footer
1314 archiveFile << "\t</VTKFile>" << endl;
1315
1316 //
1317 // Close archive file
1318 archiveFile.close();
1319
1320 if (!archiveFile.good()) {
1321 throw BruceException("Error in Bruce::saveVTK: problem closing file");
1322 }
1323
1324 }
1325
1326 void
1327 Bruce::interpolateOnDomain(escript::Data& target,
1328 const escript::Data& source) const
1329 {
1330 }
1331
1332 bool
1333 Bruce::probeInterpolationOnDomain(int functionSpaceType_source,
1334 int functionSpaceType_target) const
1335 {
1336 return true;
1337 }
1338
1339 void
1340 Bruce::interpolateACross(escript::Data& target,
1341 const escript::Data& source) const
1342 {
1343 }
1344
1345 bool
1346 Bruce::probeInterpolationACross(int functionSpaceType_source,
1347 const AbstractDomain& targetDomain,
1348 int functionSpaceType_target) const
1349 {
1350 return true;
1351 }
1352
1353 bool
1354 Bruce::isZero(DimVec vec)
1355 {
1356 for (int i=0; i<vec.size(); i++) {
1357 if (vec[i] != 0) {
1358 return false;
1359 }
1360 }
1361 return true;
1362 }
1363
1364 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26