/[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 154 - (show annotations)
Mon Nov 7 05:51:17 2005 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/bruce/src/Bruce/Bruce.cpp
File size: 35271 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-11-07

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26