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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 209 - (show annotations)
Wed Nov 23 06:32:25 2005 UTC (15 years ago) by robwdcock
File size: 35551 byte(s)
PARTIAL WIN32 BUILD SYSTEM AND PORT
+ bruce, escript build system and library now ported
+ other no longer necessary directories removed from this branches subversion repository

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26