/[escript]/branches/RW_WIN32/bruce/src/Bruce.cpp
ViewVC logotype

Annotation of /branches/RW_WIN32/bruce/src/Bruce.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 210 - (hide annotations)
Wed Nov 23 09:54:02 2005 UTC (15 years ago) by robwdcock
File size: 36072 byte(s)
PARTIAL WIN32 BUILD SYSTEM AND PORT
+ All libraries now build under new build system. No unit test routines yet
+ Reversed some incorrect refactorings in Bruce.cpp

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26