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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 475 - (show annotations)
Mon Jan 30 05:36:15 2006 UTC (14 years, 9 months ago) by jgs
File size: 35313 byte(s)
rationalise #includes

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26