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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 465 - (show annotations)
Wed Jan 25 01:08:17 2006 UTC (14 years, 10 months ago) by jgs
File size: 35186 byte(s)
reorganise bruce source tree:
move all from src/Bruce -> src
remove inc
adjust all #includes appropriately


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26