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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 1 month ago) by ksteube
File size: 34837 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26