/[escript]/branches/intelc_win32/bruce/src/Bruce.cpp
ViewVC logotype

Contents of /branches/intelc_win32/bruce/src/Bruce.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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