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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 11 months ago) by jgs
File size: 35271 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26