/[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 210 - (show annotations)
Wed Nov 23 09:54:02 2005 UTC (15 years ago) by robwdcock
File size: 36072 byte(s)
PARTIAL WIN32 BUILD SYSTEM AND PORT
+ All libraries now build under new build system. No unit test routines yet
+ Reversed some incorrect refactorings in Bruce.cpp

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 #include "BruceException.h"
18
19 #include "Paso/Common.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 /* win32 refactor */
902 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
903 for(int i=0;i<num_data;i++)
904 {
905 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
906 }
907
908 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
909 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
910 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
911
912 boost::python::list keys=dataDict.keys();
913 for (int i=0; i<num_data; i++) {
914 Data& d=boost::python::extract<Data&>(dataDict[keys[i]]);
915 if (dynamic_cast<const Bruce&>(d.getFunctionSpace().getDomain()) != *this) {
916 throw BruceException("Error: Bruce::saveVTK: Data must be defined on same Domain");
917 }
918 data[i]=d.getDataC();
919 ptr_data[i]=&(data[i]);
920 std::string n=boost::python::extract<std::string>(keys[i]);
921 c_names[i]=&(names[i][0]);
922 if (n.length()>MAX_namelength-1) {
923 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
924 c_names[i][MAX_namelength-1]='\0';
925 } else {
926 strcpy(c_names[i],n.c_str());
927 }
928 }
929
930 //
931 // open archive file
932 ofstream archiveFile;
933 archiveFile.open(filename.data(), ios::out);
934 if (!archiveFile.good()) {
935 throw BruceException("Error in Bruce::saveVTK: File could not be opened for writing");
936 }
937
938 //
939 // determine mesh type to be written
940 std::vector<bool> isCellCentered(num_data); // MSVC refactor of c99 array construct
941 bool write_celldata=false, write_pointdata=false;
942 for (int i_data=0; i_data<num_data; i_data++) {
943 if (!isEmpty(ptr_data[i_data])) {
944 switch(getFunctionSpaceType(ptr_data[i_data])) {
945 case ContinuousFunction:
946 isCellCentered[i_data]=false;
947 break;
948 case Function:
949 isCellCentered[i_data]=true;
950 break;
951 }
952 if (isCellCentered[i_data]) {
953 write_celldata=true;
954 } else {
955 write_pointdata=true;
956 }
957 }
958 }
959
960 //
961 // determine number of points and cells
962 int numPoints=getDataShape(ContinuousFunction).second;
963 int numCells=getDataShape(Function).second;
964
965 //
966 // determine VTK element type
967 int cellType;
968 char elemTypeStr[32];
969 int numVTKNodesPerElement;
970
971 int nDim = getDim();
972 switch (nDim) {
973
974 case 0:
975 cellType = VTK_VERTEX;
976 numVTKNodesPerElement = 1;
977 strcpy(elemTypeStr, "VTK_VERTEX");
978 break;
979
980 case 1:
981 cellType = VTK_LINE;
982 numVTKNodesPerElement = 2;
983 strcpy(elemTypeStr, "VTK_LINE");
984 break;
985
986 case 2:
987 cellType = VTK_QUAD;
988 numVTKNodesPerElement = 4;
989 strcpy(elemTypeStr, "VTK_QUAD");
990 break;
991
992 case 3:
993 cellType = VTK_HEXAHEDRON;
994 numVTKNodesPerElement = 8;
995 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
996 break;
997
998 }
999
1000 //
1001 // write xml header
1002 archiveFile << "<?xml version=\"1.0\"?>" << endl;
1003
1004 //
1005 // determine grid extent
1006
1007 // ??
1008
1009 //
1010 // write grid type and extent
1011 archiveFile << "\t<VTKFile type=\"StructuredGrid\" version=\"0.1\">" << endl;
1012 archiveFile << "\t\t<StructuredGrid WholeExtent=\"x1 x2 y1 y2 z1 z2\">" << endl;
1013 archiveFile << "\t\t\t<Piece Extent=\"x1 x2 y1 y2 z1 z2\">" << endl;
1014
1015 //
1016 // start to write out point definitions
1017 archiveFile << "\t\t\t\t<Points>" << endl;
1018
1019 //
1020 // determine grid cooordinates
1021
1022 // ??
1023
1024 // vtk/mayavi doesn't like 2D data, it likes 3D data with
1025 // a degenerate third dimension to handle 2D data.
1026 // So if nDim is less than 3, must pad all empty dimensions,
1027 // so that the total number of dims is 3.
1028
1029 archiveFile << "\t\t\t\t\t<DataArray NumberOfComponents=" << max(3,nDim) << " type=\"Float32\" format=\"ascii\">" << endl;
1030 /*
1031 for (int i=0; i<numPoints; i++) {
1032 for (int j=0; j<nDim; j++)
1033 archiveFile << "%e "; //mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)])
1034 for (int k=0; !(k>=3-nDim); k++)
1035 archiveFile << "0. ";
1036 archiveFile << endl;
1037 }
1038 */
1039 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1040 archiveFile << "\t\t\t\t</Points>" << endl;
1041
1042 //
1043 // cell data
1044 if (write_celldata) {
1045
1046 //
1047 // mark the active cell-data data arrays
1048 bool set_scalar=false, set_vector=false, set_tensor=false;
1049 archiveFile << "\t\t\t\t<CellData";
1050 for (int i_data=0; i_data<num_data; i_data++) {
1051 if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1052
1053 switch(getDataPointRank(ptr_data[i_data])) {
1054
1055 case 0:
1056 if (!set_scalar) {
1057 archiveFile << " Scalars=" << names[i_data];
1058 set_scalar=true;
1059 }
1060 break;
1061
1062 case 1:
1063 if (!set_vector) {
1064 archiveFile << " Vectors=" << names[i_data];
1065 set_vector=true;
1066 }
1067 break;
1068
1069 case 2:
1070 if (!set_tensor) {
1071 archiveFile << " Tensors=" << names[i_data];
1072 set_tensor=true;
1073 }
1074 break;
1075
1076 default:
1077 archiveFile.close();
1078 throw BruceException("saveVTK: VTK can't handle objects with rank greater than 2.");
1079
1080 }
1081 }
1082 }
1083 archiveFile << ">" << endl;
1084
1085 //
1086 // write the cell-data data arrays
1087 for (int i_data=0; i_data<num_data; i_data++) {
1088 if (!isEmpty(ptr_data[i_data]) && isCellCentered[i_data]) {
1089
1090 int numPointsPerSample=1;
1091 int rank=getDataPointRank(ptr_data[i_data]);
1092 int nComp=getDataPointSize(ptr_data[i_data]);
1093 int nCompReqd; // the number of components required by vtk
1094 int shape=0;
1095
1096 switch (rank) {
1097
1098 case 0:
1099 nCompReqd=1;
1100 break;
1101
1102 case 1:
1103 shape=getDataPointShape(ptr_data[i_data], 0);
1104 if (shape>3) {
1105 archiveFile.close();
1106 throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1107 }
1108 nCompReqd=3;
1109 break;
1110
1111 case 2:
1112 shape=getDataPointShape(ptr_data[i_data], 0);
1113 if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1114 archiveFile.close();
1115 throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1116 }
1117 nCompReqd=9;
1118 break;
1119
1120 }
1121
1122 archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1123
1124 /*
1125 //
1126 // write out the cell data
1127
1128 // if the number of required components is more than the number
1129 // of actual components, pad with zeros
1130 // probably only need to get shape of first element
1131 // write the data different ways for scalar, vector and tensor
1132
1133 for (int i=0; i<numCells; i++) {
1134
1135 double *values=getSampleData(ptr_data[i_data], i);
1136 double sampleAvg[nComp];
1137
1138 // average over the number of points in the sample (ie: 1)
1139 for (int k=0; k<nComp; k++) {
1140 sampleAvg[k] = values[INDEX2(k,0,nComp)];
1141 }
1142
1143 if (nCompReqd == 1) {
1144
1145 archiveFile << sampleAvg[0] << endl;
1146
1147 } else if (nCompReqd == 3) {
1148
1149 for (int m=0; m<shape; m++) {
1150 archiveFile << sampleAvg[m] << endl;
1151 }
1152 for (int m=0; m<nCompReqd-shape; m++) {
1153 archiveFile << "0." << endl;
1154 }
1155
1156 } else if (nCompReqd == 9) {
1157
1158 int count=0;
1159 for (int m=0; m<shape; m++) {
1160 for (int n=0; n<shape; n++) {
1161 archiveFile << sampleAvg[count] << endl;
1162 count++;
1163 }
1164 for (int n=0; n<3-shape; n++) {
1165 archiveFile << "0." << endl;
1166 }
1167 }
1168 for (int m=0; m<3-shape; m++) {
1169 for (int n=0; n<3; n++) {
1170 archiveFile << "0." << endl;
1171 }
1172 }
1173
1174 }
1175 }
1176 */
1177 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1178
1179 }
1180 }
1181
1182 archiveFile << "\t\t\t\t</CellData>" << endl;
1183 }
1184
1185 //
1186 // point data
1187 if (write_pointdata) {
1188
1189 //
1190 // mark the active point-data data arrays
1191 bool set_scalar=false, set_vector=false, set_tensor=false;
1192 archiveFile << "\t\t\t\t<PointData";
1193 for (int i_data=0; i_data<num_data; i_data++) {
1194 if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1195
1196 switch(getDataPointRank(ptr_data[i_data])) {
1197
1198 case 0:
1199 if (!set_scalar) {
1200 archiveFile << " Scalars=" << names[i_data];
1201 set_scalar=true;
1202 }
1203 break;
1204
1205 case 1:
1206 if (!set_vector) {
1207 archiveFile << " Vectors=" << names[i_data];
1208 set_vector=true;
1209 }
1210 break;
1211
1212 case 2:
1213 if (!set_tensor) {
1214 archiveFile << " Tensors=" << names[i_data];
1215 set_tensor=true;
1216 }
1217 break;
1218
1219 default:
1220 archiveFile.close();
1221 throw BruceException("saveVTK: Vtk can't handle objects with rank greater than 2");
1222 }
1223
1224 }
1225 }
1226 archiveFile << ">" << endl;
1227
1228 //
1229 // write the point-data data arrays
1230 for (int i_data=0; i_data<num_data; i_data++) {
1231 if (!isEmpty(ptr_data[i_data]) && !isCellCentered[i_data]) {
1232
1233 int numPointsPerSample=1;
1234 int rank=getDataPointRank(ptr_data[i_data]);
1235 int nComp=getDataPointSize(ptr_data[i_data]);
1236 int nCompReqd; // the number of components required by vtk
1237 int shape=0;
1238
1239 switch (rank) {
1240
1241 case 0:
1242 nCompReqd=1;
1243
1244 case 1:
1245 shape=getDataPointShape(ptr_data[i_data], 0);
1246 if (shape>3) {
1247 archiveFile.close();
1248 throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1249 }
1250 nCompReqd=3;
1251
1252 case 2:
1253 shape=getDataPointShape(ptr_data[i_data], 0);
1254 if (shape>3 || shape!=getDataPointShape(ptr_data[i_data], 1)) {
1255 archiveFile.close();
1256 throw BruceException("saveVTK: rank 2 object must have less than 4x4 components and must have a square shape");
1257 }
1258 nCompReqd=9;
1259
1260 }
1261
1262 archiveFile << "\t\t\t\t\t<DataArray Name=" << names[i_data] << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1263
1264 /*
1265 //
1266 // write out the point data
1267
1268 // if the number of required components is more than
1269 // the number of actual components, pad with zeros
1270 // write the data different ways for scalar, vector and tensor
1271
1272 for (int i=0; i<numNodes; i++) {
1273
1274 values=getSampleData(data[i_data], i);
1275
1276 if (nCompReqd==1) {
1277
1278 archiveFile << values[0];
1279
1280 } else if (nCompReqd==3) {
1281
1282 for (int m=0; m<shape; m++) {
1283 archiveFile << values[m];
1284 }
1285 for (int m=0; m<nCompReqd-shape; m++) {
1286 archiveFile << "0.";
1287 }
1288
1289 } else if (nCompReqd==9) {
1290
1291 int count=0;
1292 for (int m=0; m<shape; m++) {
1293 for (int n=0; n<shape; n++) {
1294 archiveFile << values[count];
1295 count++;
1296 }
1297 for (int n=0; n<3-shape; n++) {
1298 archiveFile << "0.";
1299 }
1300 }
1301 for (int m=0; m<3-shape; m++) {
1302 for (int n=0; n<3; n++) {
1303 archiveFile << "0.";
1304 }
1305 }
1306
1307 }
1308
1309 archiveFile << endl;
1310 }
1311 */
1312
1313 archiveFile << "\t\t\t\t\t</DataArray>" << endl;
1314 }
1315 }
1316
1317 archiveFile << "\t\t\t\t</PointData>" << endl;
1318 }
1319
1320 //
1321 // finish off the grid definition
1322 archiveFile << "\t\t\t</Piece>" << endl;
1323 archiveFile << "\t\t</StructuredGrid>" << endl;
1324
1325 //
1326 // write the xml footer
1327 archiveFile << "\t</VTKFile>" << endl;
1328
1329 //
1330 // Close archive file
1331 archiveFile.close();
1332
1333 /* win32 refactor */
1334 TMPMEMFREE(c_names);
1335 TMPMEMFREE(data);
1336 TMPMEMFREE(ptr_data);
1337 for(int i=0;i<num_data;i++)
1338 {
1339 TMPMEMFREE(names[i]);
1340 }
1341 TMPMEMFREE(names);
1342
1343 if (!archiveFile.good()) {
1344 throw BruceException("Error in Bruce::saveVTK: problem closing file");
1345 }
1346
1347 }
1348
1349 void
1350 Bruce::interpolateOnDomain(escript::Data& target,
1351 const escript::Data& source) const
1352 {
1353 }
1354
1355 bool
1356 Bruce::probeInterpolationOnDomain(int functionSpaceType_source,
1357 int functionSpaceType_target) const
1358 {
1359 return true;
1360 }
1361
1362 void
1363 Bruce::interpolateACross(escript::Data& target,
1364 const escript::Data& source) const
1365 {
1366 }
1367
1368 bool
1369 Bruce::probeInterpolationACross(int functionSpaceType_source,
1370 const AbstractDomain& targetDomain,
1371 int functionSpaceType_target) const
1372 {
1373 return true;
1374 }
1375
1376 bool
1377 Bruce::isZero(DimVec vec)
1378 {
1379 for (int i=0; i<vec.size(); i++) {
1380 if (vec[i] != 0) {
1381 return false;
1382 }
1383 }
1384 return true;
1385 }
1386
1387 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26