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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26