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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/bruce/src/Bruce/Bruce.cpp revision 149 by jgs, Thu Sep 1 03:31:39 2005 UTC trunk/bruce/src/Bruce/Bruce.cpp revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC
# Line 14  Line 14 
14  */  */
15    
16  #include "bruce/Bruce/Bruce.h"  #include "bruce/Bruce/Bruce.h"
17    #include "bruce/Bruce/BruceException.h"
18    
19    #include "vtkCellType.h"
20    #include <boost/python/extract.hpp>
21    
22  using namespace std;  using namespace std;
23  using namespace escript;  using namespace escript;
24    
25  namespace bruce {  namespace bruce {
26    
27  const int Bruce::Nodes=0;  const int Bruce::ContinuousFunction=0;
28  const int Bruce::Elements=1;  const int Bruce::Function=1;
29    
30    Bruce::FunctionSpaceNamesMapType Bruce::m_functionSpaceTypeNames;
31    
32  Bruce::Bruce()  Bruce::Bruce()
33  {  {
34      setFunctionSpaceTypeNames();
35    }
36    
37    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    {
44      if (!checkParameters()) {
45        stringstream temp;
46        temp << "Error - Invalid parameters supplied to constructor.";
47        throw BruceException(temp.str());
48      }
49      setFunctionSpaceTypeNames();
50  }  }
51    
52  Bruce::Bruce(const Bruce& other)  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  Bruce::~Bruce()  Bruce::~Bruce()
61  {  {
62      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    }
70    
71    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    bool
81    Bruce::checkParameters()
82    {
83      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      // reorder vectors and point counts according to point counts
98    
99      //
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  }  }
274    
275  bool  bool
276  Bruce::isValidFunctionSpaceType(int functionSpaceType) const  Bruce::isValidFunctionSpaceType(int functionSpaceCode) const
277  {  {
278    return (true);    FunctionSpaceNamesMapType::iterator loc;
279      loc=m_functionSpaceTypeNames.find(functionSpaceCode);
280      return (loc!=m_functionSpaceTypeNames.end());
281  }  }
282    
283  int  std::string
284  Bruce::getDim() const  Bruce::functionSpaceTypeAsString(int functionSpaceCode) const
285  {  {
286    return 0;    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  pair<int,int>  pair<int,int>
298  Bruce::getDataShape(int functionSpaceCode) const  Bruce::getDataShape(int functionSpaceCode) const
299  {  {
300    int numDataPointsPerSample=0;    int numDataPointsPerSample=1;
301    int numSamples=0;    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    return pair<int,int>(numDataPointsPerSample,numSamples);    return pair<int,int>(numDataPointsPerSample,numSamples);
368  }  }
369    
370    int
371    Bruce::getNumSamples(int functionSpaceCode) const
372    {
373      std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
374      return domainShape.second;
375    }
376    
377    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      if (dim>0 && (dataShape.size()!=1 || dataShape[0]!=dim)) {
412        stringstream temp;
413        temp << "Error - Incompatible shape Data object supplied to Bruce::setToX";
414        throw BruceException(temp.str());
415      } 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      }
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            int sampleNo=i;
440            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              int sampleNo=(m_n1*i)+j;
451              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                int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
466                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        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              int sampleNo=i;
504              DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
505              sampleData[0] = m_origin[0] + m_v0[0]*(i + 0.5);
506            }
507          }
508    
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              int sampleNo=i;
523              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                int sampleNo=(n1max*i)+j;
532                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              int sampleNo=i;
555              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                int sampleNo=(n1max*i)+j;
564                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                  int sampleNo=(n1max*n2max*i)+(n2max*j)+k;
575                  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        }
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      //
607      // determine functionSpace type expected by supplied Data object
608      int functionSpaceCode = out.getFunctionSpace().getTypeCode();
609    
610      //
611      // ensure supplied Data object has a matching number of data-points
612      // for this Bruce domain object
613      std::pair<int,int> domainShape = getDataShape(functionSpaceCode);
614      if(domainShape.first!=out.getNumDataPointsPerSample() ||
615         domainShape.second!=out.getNumSamples()) {
616        stringstream temp;
617        temp << "Error - Incompatible number of data-points Data object supplied to Bruce::setToSize";
618        throw BruceException(temp.str());
619      }
620    
621      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      std::vector<int> dataShape = out.getDataPointShape();
628      // 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      if (dataShape.size()!=1 || dataShape[0]!=1) {
633        stringstream temp;
634        temp << "Error - Incompatible shape Data object supplied to Bruce::setToSize";
635        throw BruceException(temp.str());
636      }
637    */
638    
639      double dp_size;
640    
641      if (functionSpaceCode==ContinuousFunction) {
642    
643        //
644        // calculate the size of data-points
645        // located on the nodes of this Bruce domain
646    
647        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            int sampleNo=i;
665            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              int sampleNo=(m_n1*i)+j;
680              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                int sampleNo=(m_n1*m_n2*i)+(m_n2*j)+k;
698                DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
699                sampleData[0] = dp_size;
700              }
701            }
702          }
703    
704        }
705    
706      } else if (functionSpaceCode==Function) {
707    
708        //
709        // calculate the sizes of data-points
710        // located on the cell centres of this Bruce domain
711    
712        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              int sampleNo=i;
736              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              int sampleNo=i;
759              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                int sampleNo=(n1max*i)+j;
766                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              int sampleNo=i;
792              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                int sampleNo=(n1max*i)+j;
799                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                  int sampleNo=(n2max*n1max*i)+(n1max*j)+k;
808                  DataAbstract::ValueType::value_type* sampleData = out.getSampleData(sampleNo);
809                  sampleData[0] = dp_size;
810                }
811              }
812            }
813          }
814    
815        }
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    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  bool  bool
838  Bruce::operator==(const AbstractDomain& other) const  Bruce::operator==(const AbstractDomain& other) const
839  {  {
840    const Bruce* temp=dynamic_cast<const Bruce*>(&other);    const Bruce* temp=dynamic_cast<const Bruce*>(&other);
841    if (temp!=0) {    if (temp!=0) {
842      return (true);      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    } else {    } else {
853      return false;      return false;
854    }    }
# Line 72  Bruce::operator!=(const AbstractDomain& Line 860  Bruce::operator!=(const AbstractDomain&
860    return !(operator==(other));    return !(operator==(other));
861  }  }
862    
863  Data  int
864  Bruce::getX() const  Bruce::getReferenceNoFromSampleNo(int functionSpaceCode,
865                                      int sampleNo) const
866  {  {
867    return continuousFunction(asAbstractContinuousDomain()).getX();    // 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  }  }
878    
879  Data  void
880  Bruce::getSize() const  Bruce::saveVTK(const std::string& filename,
881                   const boost::python::dict& dataDict) const
882    {
883    
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    }
1320    
1321    void
1322    Bruce::interpolateOnDomain(escript::Data& target,
1323                               const escript::Data& source) const
1324    {
1325    }
1326    
1327    bool
1328    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 function(asAbstractContinuousDomain()).getSize();    return true;
1346    }
1347    
1348    bool
1349    Bruce::isZero(DimVec vec)
1350    {
1351      for (int i=0; i<vec.size(); i++) {
1352        if (vec[i] != 0) {
1353          return false;
1354        }
1355      }
1356      return true;
1357  }  }
1358    
1359  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.149  
changed lines
  Added in v.155

  ViewVC Help
Powered by ViewVC 1.1.26