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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 11 months 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 jgs 149
2 ksteube 1312 /* $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 jgs 465 #include "Bruce.h"
17 jgs 149
18 jgs 475 #include "BruceException.h"
19    
20 robwdcock 682 #include "escript/FunctionSpaceFactory.h"
21 jgs 475
22     #include <boost/python/extract.hpp>
23 woo409 757 #include <vector>
24     #include <string>
25 jgs 475 #include "vtkCellType.h"
26    
27 jgs 149 using namespace std;
28     using namespace escript;
29    
30     namespace bruce {
31    
32 woo409 757 BRUCE_DLL_API
33 jgs 150 const int Bruce::ContinuousFunction=0;
34 woo409 757 BRUCE_DLL_API
35 jgs 150 const int Bruce::Function=1;
36 jgs 149
37 jgs 150 Bruce::FunctionSpaceNamesMapType Bruce::m_functionSpaceTypeNames;
38    
39 jgs 372 Bruce::Bruce():
40     m_n0(0), m_n1(0), m_n2(0)
41 jgs 149 {
42 jgs 150 setFunctionSpaceTypeNames();
43 jgs 149 }
44    
45 jgs 150 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 jgs 149 {
52 jgs 150 if (!checkParameters()) {
53     stringstream temp;
54     temp << "Error - Invalid parameters supplied to constructor.";
55     throw BruceException(temp.str());
56     }
57     setFunctionSpaceTypeNames();
58 jgs 149 }
59    
60 jgs 150 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 jgs 149 Bruce::~Bruce()
69     {
70 jgs 150 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 jgs 149 }
78    
79 jgs 150 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 jgs 149 bool
89 jgs 150 Bruce::checkParameters()
90 jgs 149 {
91 jgs 150 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 jgs 153 // reorder vectors and point counts according to point counts
106    
107 jgs 150 //
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 jgs 149 }
282    
283 jgs 150 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 jgs 149 pair<int,int>
306     Bruce::getDataShape(int functionSpaceCode) const
307     {
308 jgs 150 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 jgs 149 return pair<int,int>(numDataPointsPerSample,numSamples);
376     }
377    
378 jgs 154 int
379     Bruce::getNumSamples(int functionSpaceCode) const
380     {
381     std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
382     return domainShape.second;
383     }
384    
385 jgs 150 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 jgs 151 if (dim>0 && (dataShape.size()!=1 || dataShape[0]!=dim)) {
420 jgs 150 stringstream temp;
421     temp << "Error - Incompatible shape Data object supplied to Bruce::setToX";
422     throw BruceException(temp.str());
423 jgs 151 } 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 jgs 150 }
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 jgs 153 int sampleNo=i;
448 jgs 150 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 jgs 153 int sampleNo=(m_n1*i)+j;
459 jgs 150 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 jgs 153 int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
474 jgs 150 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 jgs 151 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 jgs 153 int sampleNo=i;
512 jgs 151 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
513     sampleData[0] = m_origin[0] + m_v0[0]*(i + 0.5);
514     }
515 jgs 150 }
516 jgs 151
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 jgs 153 int sampleNo=i;
531 jgs 151 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 jgs 153 int sampleNo=(n1max*i)+j;
540 jgs 151 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 jgs 153 int sampleNo=i;
563 jgs 151 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 jgs 153 int sampleNo=(n1max*i)+j;
572 jgs 151 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 jgs 153 int sampleNo=(n1max*n2max*i)+(n2max*j)+k;
583 jgs 151 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 jgs 150 }
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 jgs 151 //
615     // determine functionSpace type expected by supplied Data object
616 jgs 150 int functionSpaceCode = out.getFunctionSpace().getTypeCode();
617    
618 jgs 151 //
619     // ensure supplied Data object has a matching number of data-points
620     // for this Bruce domain object
621 jgs 150 std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
622     if(domainShape.first!=out.getNumDataPointsPerSample() ||
623     domainShape.second!=out.getNumSamples()) {
624     stringstream temp;
625 jgs 151 temp << "Error - Incompatible number of data-points Data object supplied to Bruce::setToSize";
626 jgs 150 throw BruceException(temp.str());
627     }
628    
629 jgs 151 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 jgs 150 std::vector<int> dataShape = out.getDataPointShape();
636 jgs 153 // 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 jgs 151 if (dataShape.size()!=1 || dataShape[0]!=1) {
641 jgs 150 stringstream temp;
642 jgs 151 temp << "Error - Incompatible shape Data object supplied to Bruce::setToSize";
643 jgs 150 throw BruceException(temp.str());
644     }
645 jgs 153 */
646 jgs 150
647 jgs 151 double dp_size;
648 jgs 150
649     if (functionSpaceCode==ContinuousFunction) {
650    
651 jgs 151 //
652     // calculate the size of data-points
653     // located on the nodes of this Bruce domain
654 jgs 150
655 jgs 151 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 jgs 153 int sampleNo=i;
673 jgs 151 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 jgs 153 int sampleNo=(m_n1*i)+j;
688 jgs 151 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 jgs 153 int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
706 jgs 151 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
707     sampleData[0] = dp_size;
708     }
709     }
710     }
711    
712 jgs 150 }
713    
714     } else if (functionSpaceCode==Function) {
715    
716 jgs 151 //
717     // calculate the sizes of data-points
718     // located on the cell centres of this Bruce domain
719 jgs 150
720 jgs 151 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 jgs 153 int sampleNo=i;
744 jgs 151 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 jgs 153 int sampleNo=i;
767 jgs 151 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 jgs 153 int sampleNo=(n1max*i)+j;
774 jgs 151 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 jgs 153 int sampleNo=i;
800 jgs 151 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 jgs 153 int sampleNo=(n1max*i)+j;
807 jgs 151 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 jgs 153 int sampleNo=(n2max*n1max*i)+(n1max*j)+k;
816 jgs 151 DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
817     sampleData[0] = dp_size;
818     }
819     }
820     }
821     }
822    
823 jgs 150 }
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 jgs 154 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 jgs 149 bool
846     Bruce::operator==(const AbstractDomain& other) const
847     {
848     const Bruce* temp=dynamic_cast<const Bruce*>(&other);
849     if (temp!=0) {
850 jgs 150 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 jgs 149 } else {
861     return false;
862     }
863     }
864    
865     bool
866     Bruce::operator!=(const AbstractDomain& other) const
867     {
868     return !(operator==(other));
869     }
870    
871 jgs 153 int
872 jgs 154 Bruce::getReferenceNoFromSampleNo(int functionSpaceCode,
873     int sampleNo) const
874 jgs 153 {
875 jgs 154 // 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 jgs 153 }
886    
887 jgs 154 void
888     Bruce::saveVTK(const std::string& filename,
889     const boost::python::dict& dataDict) const
890 jgs 153 {
891 jgs 154
892     const int num_data=boost::python::extract<int>(dataDict.attr("__len__")());
893     boost::python::list keys=dataDict.keys();
894    
895     //
896 woo409 757 // extract Data objects and associated names from dictionary object
897 phornby 870 const std::string::size_type MAX_namelength=256;
898 woo409 757
899     //
900 jgs 154 // 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 woo409 757 bool write_celldata=false, write_pointdata=false;
909 jgs 154 for (int i_data=0; i_data<num_data; i_data++) {
910 woo409 757 Data& d=boost::python::extract<Data&>(dataDict[keys[i_data]]);
911     if (!d.isEmpty()) {
912     switch(d.getFunctionSpace().getTypeCode()) {
913 jgs 154 case ContinuousFunction:
914 woo409 757 write_pointdata=true;
915 jgs 154 break;
916     case Function:
917 woo409 757 write_celldata=true;
918 jgs 154 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 woo409 757 std::string elemTypeStr;
932 jgs 154 int numVTKNodesPerElement;
933    
934     int nDim = getDim();
935     switch (nDim) {
936    
937     case 0:
938     cellType = VTK_VERTEX;
939     numVTKNodesPerElement = 1;
940 woo409 757 elemTypeStr = "VTK_VERTEX";
941 jgs 154 break;
942    
943     case 1:
944     cellType = VTK_LINE;
945     numVTKNodesPerElement = 2;
946 woo409 757 elemTypeStr = "VTK_LINE";
947 jgs 154 break;
948    
949     case 2:
950     cellType = VTK_QUAD;
951     numVTKNodesPerElement = 4;
952 woo409 757 elemTypeStr = "VTK_QUAD";
953 jgs 154 break;
954    
955     case 3:
956     cellType = VTK_HEXAHEDRON;
957     numVTKNodesPerElement = 8;
958 woo409 757 elemTypeStr = "VTK_HEXAHEDRON";
959 jgs 154 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 woo409 757 Data& d = boost::python::extract<Data&>(dataDict[keys[i_data]]);
1014 phornby 870 const std::string &
1015     full_name = boost::python::extract<std::string>(keys[i_data]);
1016     std::string name(full_name, 0, MAX_namelength);
1017 woo409 757
1018     if (!d.isEmpty() && d.getFunctionSpace().getTypeCode()==Function) {
1019    
1020     switch(d.getDataPointRank()) {
1021 jgs 154
1022     case 0:
1023     if (!set_scalar) {
1024 woo409 757 archiveFile << " Scalars=" << name;
1025 jgs 154 set_scalar=true;
1026     }
1027     break;
1028    
1029     case 1:
1030     if (!set_vector) {
1031 woo409 757 archiveFile << " Vectors=" << name;
1032 jgs 154 set_vector=true;
1033     }
1034     break;
1035    
1036     case 2:
1037     if (!set_tensor) {
1038 woo409 757 archiveFile << " Tensors=" << name;
1039 jgs 154 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 woo409 757 Data& d = boost::python::extract<Data&>(dataDict[keys[i_data]]);
1055 phornby 870 const std::string &
1056     full_name = boost::python::extract<std::string>(keys[i_data]);
1057     std::string name(full_name, 0, MAX_namelength);
1058 woo409 757
1059     if (!d.isEmpty() && d.getFunctionSpace().getTypeCode()==Function) {
1060 jgs 154 int numPointsPerSample=1;
1061 woo409 757 int rank = d.getDataPointRank();
1062     int nComp = d.getDataPointSize();
1063 jgs 154 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 woo409 757 shape = d.getDataPointShape()[0];
1074 jgs 154 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 woo409 757 shape=d.getDataPointShape()[0];
1083     if (shape>3 || shape!=d.getDataPointShape()[1]) {
1084 jgs 154 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 woo409 757 archiveFile << "\t\t\t\t\t<DataArray Name=" << name << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1093 jgs 154 /*
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 woo409 757 Data& d = boost::python::extract<Data&>(dataDict[keys[i_data]]);
1162 phornby 870 const std::string &
1163     full_name = boost::python::extract<std::string>(keys[i_data]);
1164     std::string name(full_name, 0, MAX_namelength);
1165 jgs 154
1166 woo409 757 if (!d.isEmpty() && d.getFunctionSpace().getTypeCode()==ContinuousFunction) {
1167 jgs 154
1168 woo409 757 switch(d.getDataPointRank()) {
1169    
1170 jgs 154 case 0:
1171     if (!set_scalar) {
1172 woo409 757 archiveFile << " Scalars=" << name;
1173 jgs 154 set_scalar=true;
1174     }
1175     break;
1176    
1177     case 1:
1178     if (!set_vector) {
1179 woo409 757 archiveFile << " Vectors=" << name;
1180 jgs 154 set_vector=true;
1181     }
1182     break;
1183    
1184     case 2:
1185     if (!set_tensor) {
1186 woo409 757 archiveFile << " Tensors=" << name;
1187 jgs 154 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 woo409 757 Data& d = boost::python::extract<Data&>(dataDict[keys[i_data]]);
1203 phornby 870 const std::string &
1204     full_name = boost::python::extract<std::string>(keys[i_data]);
1205     std::string name(full_name, 0, MAX_namelength);
1206 jgs 154
1207 phornby 870 if (!d.isEmpty() &&
1208     d.getFunctionSpace().getTypeCode()==ContinuousFunction) {
1209    
1210 jgs 154 int numPointsPerSample=1;
1211 woo409 757 int rank = d.getDataPointRank();
1212     int nComp = d.getDataPointSize();
1213 jgs 154 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 woo409 757 break;
1221 jgs 154 case 1:
1222 woo409 757 shape=d.getDataPointShape()[0];
1223 jgs 154 if (shape>3) {
1224     archiveFile.close();
1225     throw BruceException("saveVTK: rank 1 object must have less than 4 components");
1226     }
1227     nCompReqd=3;
1228 woo409 757 break;
1229 jgs 154 case 2:
1230 woo409 757 shape=d.getDataPointShape()[0];
1231     if (shape>3 || shape!=d.getDataPointShape()[1]) {
1232 jgs 154 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 woo409 757 break;
1237 jgs 154 }
1238    
1239 woo409 757 archiveFile << "\t\t\t\t\t<DataArray Name=" << name << " type=\"Float32\" NumberOfComponents=" << nCompReqd << " format=\"ascii\">" << endl;
1240 jgs 154 /*
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 jgs 153 }
1311    
1312 jgs 154 void
1313     Bruce::interpolateOnDomain(escript::Data& target,
1314     const escript::Data& source) const
1315     {
1316     }
1317    
1318 jgs 150 bool
1319 jgs 154 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 jgs 151 Bruce::isZero(DimVec vec)
1341 jgs 149 {
1342 jgs 150 for (int i=0; i<vec.size(); i++) {
1343     if (vec[i] != 0) {
1344     return false;
1345     }
1346     }
1347     return true;
1348 jgs 149 }
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