← Back to branch summary

~ubuntu-branches/ubuntu/saucy/ifrit/saucy

« back to all changes in this revision

Viewing changes to core/igenericparticlesdatasubject.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Mark Hymers
  • Date: 2006-10-28 15:06:32 UTC
  • mfrom: (1.1.4 upstream) (2.1.1 etch)
  • Revision ID: james.westby@ubuntu.com-20061028150632-hyvuhvsv6zpmf5ev
Tags: 3.0.5-1
New upstream version. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
core/igenericparticlesdatasubject.cpp
29
29
#include "igenericparticlesdatasubject.h"
30
30
 
31
31
 
 
32
#include "ibuffer.h"
32
33
#include "idata.h"
33
34
#include "idatalimits.h"
34
35
#include "ierror.h"
35
36
#include "ierrorstatus.h"
36
37
#include "ieventobserver.h"
 
38
#include "ihistogrammaker.h"
37
39
#include "iparallel.h"
 
40
#include "iparallelmanager.h"
38
41
#include "iparallelworker.h"
 
42
#include "iprobefilter.h"
39
43
#include "iviewmodule.h"
40
44
 
 
45
#include <vtkCellArray.h>
 
46
#include <vtkDoubleArray.h>
41
47
#include <vtkFloatArray.h>
42
48
#include <vtkMath.h>
43
49
#include <vtkPointData.h>
44
50
#include <vtkPoints.h>
45
51
#include <vtkPolyData.h>
46
52
 
47
 
 
48
 
//
49
 
//  Helper class for parallel executions
 
53
//
 
54
//  Templates
 
55
//
 
56
#include "ibuffertemplate.h"
 
57
 
 
58
 
 
59
//
 
60
//  Helper classes for parallel execution & templates
50
61
//
51
62
class iGenericParticlesDataHelper : protected iParallelWorker
52
63
{
56
67
        iGenericParticlesDataHelper(iGenericParticlesDataSubject *subject);
57
68
 
58
69
        void ShiftData(bool paf, int dim, long n, bool per[3], float dr, float *pf, double *pd);
 
70
        void FindRange(int dim, long n, float *pf, float *amin, float *amax);
59
71
 
60
72
protected:
61
73
 
63
75
 
64
76
        iGenericParticlesDataSubject *mSubject;
65
77
 
 
78
        iBuffer<float> mMin, mMax;
 
79
 
66
80
        bool *mPeriodic;
67
 
        float *mFarr1;
68
 
        double *mDarr1;
 
81
        float *mFarr;
 
82
        double *mDarr;
69
83
 
70
 
        int mItmp1;
71
 
        bool mBtmp1;
72
 
        long mLtmp1;
73
 
        float mFtmp1;
 
84
        int mDim, mNumProcs;
 
85
        bool mPaf;
 
86
        long mTot;
 
87
        float mDr;
74
88
};
75
89
 
76
90
 
77
91
//
78
92
//  Main class
79
93
//
80
 
iGenericParticlesDataSubject::iGenericParticlesDataSubject(iDataReader *r, const iDataType &type) : iDataSubject(r,type)
 
94
iGenericParticlesDataSubject::iGenericParticlesDataSubject(iDataReader *r, const iDataType &type) : iPrimaryDataSubject(r,type)
81
95
{
82
96
        mHaveNormals = false;
83
97
 
84
98
        mDownsampleMode = 0;
85
99
        mDownsampleFactor = 1;
86
100
 
87
 
        mHelper = new iGenericParticlesDataHelper(this); IERROR_ASSERT_NULL_POINTER(mHelper);
 
101
        mHelper = new iGenericParticlesDataHelper(this); IERROR_ASSERT(mHelper);
88
102
}
89
103
 
90
104
 
159
173
        if(mHaveNormals)
160
174
        {
161
175
                vtkFloatArray *newNormals;
162
 
                newNormals = vtkFloatArray::New(); IERROR_ASSERT_NULL_POINTER(newNormals);
 
176
                newNormals = vtkFloatArray::New(); IERROR_ASSERT(newNormals);
163
177
                newNormals->SetNumberOfComponents(3);
164
178
                // Allocates and Sets MaxId
165
179
                newNormals->SetNumberOfTuples(ntot);
222
236
}
223
237
 
224
238
 
225
 
vtkIdType iGenericParticlesDataSubject::PrepareMask(vtkIdType ntot, int df, vtkIdType offm, vtkIdType offt)
 
239
vtkIdType iGenericParticlesDataSubject::PrepareMask(long ntot0, int df, vtkIdType offm, long offt)
226
240
{
227
241
        if(df < 1) df = 1;
228
242
 
229
 
        wMask.NumTotal = ntot;
230
 
        wMask.NumMasked = (ntot+df-1)/df;
 
243
        wMask.NumTotal = ntot0;
 
244
        wMask.NumMasked = (ntot0+df-1)/df;
231
245
        wMask.NumUnmasked = wMask.NumTotal - wMask.NumMasked;
232
246
        wMask.OffsetTotal = offt;
233
247
        wMask.OffsetMasked = offm;
238
252
}
239
253
 
240
254
 
 
255
bool iGenericParticlesDataSubject::ReadFortranRecordWithMask(iFile &F, long ntot0, int comp, vtkFloatArray *array, float updateStart, float updateDuration, float scale, float offset)
 
256
{
 
257
        return this->ReadFortranRecordWithMaskTemplate(F,ntot0,comp,array,updateStart,updateDuration,scale,offset);
 
258
}
 
259
 
 
260
 
 
261
bool iGenericParticlesDataSubject::ReadFortranRecordWithMask(iFile &F, long ntot0, int comp, vtkDoubleArray *array, float updateStart, float updateDuration, double scale, double offset)
 
262
{
 
263
        return this->ReadFortranRecordWithMaskTemplate(F,ntot0,comp,array,updateStart,updateDuration,scale,offset);
 
264
}
 
265
 
 
266
 
 
267
template<class T>
 
268
bool iGenericParticlesDataSubject::ReadFortranRecordWithMaskTemplate(iFile &F, long ntot0, int comp, vtkDataArray *array, float updateStart, float updateDuration, T scale, T offset)
 
269
{
 
270
        int j;
 
271
        long lrec1, lrec2;
 
272
        vtkIdType l, idm, idt;
 
273
 
 
274
        //
 
275
        //  Is component valid?
 
276
        //
 
277
        int nComp;
 
278
        if(array==0 || comp<0 || comp>=(nComp = array->GetNumberOfComponents()))
 
279
        {
 
280
                this->GetErrorStatus()->Set("Invalid call to iGenericParticlesDataSubject::ReadFortranRecordWithMask.");
 
281
                return false;
 
282
        }
 
283
 
 
284
        //
 
285
        //  Read the header
 
286
        //
 
287
        if(!this->ReadFortranHeaderFooter(F,lrec1) || lrec1!=sizeof(T)*ntot0)
 
288
        {
 
289
                this->GetErrorStatus()->Set("Corrupted data.");
 
290
                return false;
 
291
        }
 
292
 
 
293
        //
 
294
        //  Direct access to data
 
295
        //
 
296
        T *ptr = (T *)array->GetVoidPointer(0) + comp;
 
297
        if(ptr == 0)
 
298
        {
 
299
                this->GetErrorStatus()->Set("Not enough memory to create the data.");
 
300
                return false;
 
301
        }
 
302
 
 
303
        //
 
304
        // Decide how to read
 
305
        //
 
306
        vtkIdType ntot = array->GetNumberOfTuples();
 
307
        vtkIdType lpiece1, lpiece = ntot0/1000;
 
308
        if(lpiece < 1000) lpiece = 1000;
 
309
        int npieces = (ntot0+lpiece-1)/lpiece;
 
310
 
 
311
        //
 
312
        //  Create tmp array
 
313
        //
 
314
        T *d = new T[lpiece];
 
315
        if(d == 0) 
 
316
        { 
 
317
                this->GetErrorStatus()->Set("Not enough memory to create the data.");
 
318
                return false;
 
319
        }
 
320
 
 
321
        //
 
322
        //  parameters for the Progress Bar
 
323
        //
 
324
        updateDuration /= npieces;
 
325
 
 
326
        //
 
327
        //  Read piece by piece
 
328
        //
 
329
        idm = idt = 0;
 
330
        for(j=0; j<npieces; j++)
 
331
        {
 
332
                if(j < npieces-1)
 
333
                {
 
334
                        lpiece1 = lpiece;
 
335
                }
 
336
                else
 
337
                {
 
338
                        //
 
339
                        //  Correct for the last record
 
340
                        //
 
341
                        lpiece1 = ntot0 - j*lpiece;
 
342
                }
 
343
                if(!this->ReadBlock(F,d,lpiece1,updateStart+j*updateDuration,updateDuration))
 
344
                {
 
345
                        this->GetErrorStatus()->Set("Corrupted data.");
 
346
                        return false;
 
347
                }
 
348
                if(this->GetObserver()->IsAborted())
 
349
                {
 
350
                        this->GetErrorStatus()->SetAbort();
 
351
                        return false;
 
352
                }
 
353
                for(l=0; l<lpiece1; l++, idt++) if(this->IsMasked(idm,idt))
 
354
                {
 
355
#ifdef I_CHECK1
 
356
                        if(idm>=ntot || idt>=ntot0)
 
357
                        {
 
358
                                IERROR_REPORT_BUG;
 
359
                        }
 
360
#endif
 
361
                        ptr[nComp*idm] = d[l];
 
362
                        idm++;
 
363
                }
 
364
        }
 
365
 
 
366
        //
 
367
        //  Read the footer
 
368
        //
 
369
        if(!this->ReadFortranHeaderFooter(F,lrec2) || lrec1!=lrec2)
 
370
        {
 
371
                this->GetErrorStatus()->Set("Corrupted data.");
 
372
                return false;
 
373
        }
 
374
 
 
375
        delete [] d;
 
376
 
 
377
        //
 
378
        //  Do we need to scale?
 
379
        //
 
380
        if(fabs(scale) > iMath::_DoubleMin)
 
381
        {
 
382
                if(fabs(offset) > iMath::_DoubleMin)
 
383
                {
 
384
                        //
 
385
                        //  Scale as positions
 
386
                        //
 
387
                        for(l=0; l<ntot; l++)
 
388
                        {
 
389
                                ptr[nComp*l] = -1.0 + scale*(ptr[nComp*l]-offset);
 
390
                        }
 
391
                }
 
392
                else
 
393
                {
 
394
                        //
 
395
                        //  Just scale
 
396
                        //
 
397
                        for(l=0; l<ntot; l++)
 
398
                        {
 
399
                                ptr[nComp*l] *= scale;
 
400
                        }
 
401
                }
 
402
        }
 
403
 
 
404
        if(this->GetObserver()->IsAborted())
 
405
        {
 
406
                this->GetErrorStatus()->SetAbort();
 
407
                return false;
 
408
        }
 
409
 
 
410
        return true;
 
411
}
 
412
 
 
413
 
 
414
iProbeFilter* iGenericParticlesDataSubject::CreateProbeFilter(iViewSubject *vo) const
 
415
{
 
416
        return iProbeFilter::New(vo);
 
417
}
 
418
 
 
419
 
 
420
iHistogramMaker* iGenericParticlesDataSubject::CreateHistogramMaker() const
 
421
{
 
422
        return iHistogramMaker::New(this->GetViewModule());
 
423
}
 
424
 
 
425
 
 
426
void iGenericParticlesDataSubject::CreateArrays(vtkIdType ntot, int natt, bool paf, vtkPolyData *data, vtkPoints* &newPoints, vtkCellArray* &newVerts, vtkFloatArray* &newScalars)
 
427
{
 
428
        vtkPoints *oldPoints = data->GetPoints();
 
429
        vtkCellArray *oldVerts = data->GetVerts();
 
430
        vtkFloatArray *oldScalars = iRequiredCast<vtkFloatArray>(data->GetPointData()->GetScalars());
 
431
 
 
432
        int type = paf ? VTK_FLOAT : VTK_DOUBLE;
 
433
 
 
434
        if(oldPoints==0 || oldPoints->GetNumberOfPoints()!=ntot || oldPoints->GetDataType()!=type)
 
435
        {
 
436
                data->SetPoints(0);
 
437
                newPoints = vtkPoints::New(type); IERROR_ASSERT(newPoints);
 
438
                // Allocates and Sets MaxId
 
439
                newPoints->SetNumberOfPoints(ntot);
 
440
        }
 
441
        else
 
442
        {
 
443
                newPoints = oldPoints;
 
444
                newPoints->Register(0);
 
445
        }
 
446
 
 
447
        if(oldVerts==0 || oldVerts->GetNumberOfCells()!=ntot)
 
448
        {
 
449
                data->SetVerts(0);
 
450
                newVerts = vtkCellArray::New(); IERROR_ASSERT(newVerts);
 
451
                // This allocates but does not Set Max Id
 
452
                newVerts->Allocate(newVerts->EstimateSize(ntot,1));
 
453
                vtkIdType l;
 
454
                for(l=0; l<ntot; l++)
 
455
                {
 
456
                        newVerts->InsertNextCell(1);
 
457
                        newVerts->InsertCellPoint(l);
 
458
                }
 
459
        }
 
460
        else
 
461
        {
 
462
                newVerts = oldVerts;
 
463
                newVerts->Register(0);
 
464
        }
 
465
 
 
466
        if(oldScalars==0 || oldScalars->GetNumberOfTuples()!=ntot || oldScalars->GetNumberOfComponents()!=natt)
 
467
        {
 
468
                data->GetPointData()->SetScalars(0);
 
469
                newScalars = vtkFloatArray::New(); IERROR_ASSERT(newScalars);
 
470
                newScalars->SetNumberOfComponents(natt);
 
471
                newScalars->SetNumberOfTuples(ntot);
 
472
        }
 
473
        else
 
474
        {
 
475
                newScalars = oldScalars;
 
476
                newScalars->Register(0);
 
477
        }
 
478
}
 
479
 
 
480
        
 
481
void iGenericParticlesDataSubject::AttachArrays(vtkPolyData *data, vtkPoints *newPoints, vtkCellArray *newVerts, vtkFloatArray *newScalars)
 
482
{
 
483
        data->SetPoints(newPoints);
 
484
        newPoints->Delete();
 
485
 
 
486
        data->SetVerts(newVerts);
 
487
        newVerts->Delete();
 
488
 
 
489
        data->GetPointData()->SetScalars(newScalars);
 
490
        newScalars->Delete();
 
491
 
 
492
        //
 
493
        //  Attribute limits.
 
494
        //
 
495
        if(mAdjustableLimits)
 
496
        {
 
497
                vtkFloatArray *att = iRequiredCast<vtkFloatArray>(data->GetPointData()->GetScalars());
 
498
 
 
499
                if(att->GetNumberOfComponents() == this->GetLimits()->GetNumVars()) // just in case
 
500
                {
 
501
                        iBuffer<float> aMin, aMax;
 
502
                        aMin.Extend(att->GetNumberOfComponents());
 
503
                        aMax.Extend(att->GetNumberOfComponents());
 
504
 
 
505
                        mHelper->FindRange(att->GetNumberOfComponents(),att->GetNumberOfTuples(),att->GetPointer(0),aMin,aMax);
 
506
 
 
507
                        int i;
 
508
                        for(i=0; i<att->GetNumberOfComponents(); i++)
 
509
                        {
 
510
                                this->GetLimits()->SetMin(i,aMin[i]);
 
511
                                this->GetLimits()->SetMax(i,aMax[i]);
 
512
                        }
 
513
                }
 
514
        }
 
515
}
 
516
 
 
517
 
241
518
//
242
519
//  Helper class
243
520
//
249
526
 
250
527
void iGenericParticlesDataHelper::ShiftData(bool paf, int dim, long n, bool per[3], float dr, float *pf, double *pd)
251
528
{
252
 
        mBtmp1 = paf;
253
 
        mItmp1 = dim;
254
 
        mLtmp1 = n;
255
 
        mFtmp1 = dr;
256
 
        mFarr1 = pf;
257
 
        mDarr1 = pd;
 
529
        mPaf = paf;
 
530
        mDim = dim;
 
531
        mTot = n;
 
532
        mDr = dr;
 
533
        mFarr = pf;
 
534
        mDarr = pd;
258
535
        mPeriodic = per;
259
536
 
260
537
        this->ParallelExecute(1);
261
538
}
262
539
 
263
540
 
264
 
int iGenericParticlesDataHelper::ExecuteStep(int /*step*/, iParallel::ProcessorInfo &p)
 
541
void iGenericParticlesDataHelper::FindRange(int dim, long tot, float *pf, float *fMin, float *fMax)
 
542
{
 
543
        mDim = dim;
 
544
        mTot = tot;
 
545
        mFarr = pf;
 
546
 
 
547
        mNumProcs = this->GetManager()->GetNumberOfProcessors();
 
548
        mMin.Extend(mDim*mNumProcs);
 
549
        mMax.Extend(mDim*mNumProcs);
 
550
        
 
551
        this->ParallelExecute(2);
 
552
 
 
553
        int i, n;
 
554
        for(n=0; n<dim; n++)
 
555
        {
 
556
                fMin[n] = mMin[n];
 
557
                fMax[n] = mMax[n];
 
558
                for(i=1; i<mNumProcs; i++)
 
559
                {
 
560
                        if(fMin[n] > mMin[n+dim*i]) fMin[n] = mMin[n+dim*i];
 
561
                        if(fMax[n] < mMax[n+dim*i]) fMax[n] = mMax[n+dim*i];
 
562
                }
 
563
        }
 
564
}
 
565
 
 
566
 
 
567
int iGenericParticlesDataHelper::ExecuteStep(int step, iParallel::ProcessorInfo &p)
265
568
{
266
569
        long l, kstp, kbeg, kend;
267
 
 
268
 
        iParallel::SplitRange(p,mLtmp1,kbeg,kend,kstp);
269
 
 
270
 
        int d = mItmp1;
271
 
        float dr = mFtmp1;
272
 
 
273
 
        if(mBtmp1)
274
 
        {
275
 
                float *x = mFarr1 + 3*kbeg;
276
 
                for(l=kbeg; l<kend; l++)
277
 
                {
278
 
                        if(l%1000 == 0)
279
 
                        {
280
 
                                if(this->IsMaster(p)) mSubject->GetObserver()->SetProgress((d+(float)(l-kbeg)/(kend-kbeg))/3.0);
281
 
                                if(mSubject->GetObserver()->IsAborted()) return 2;
282
 
                        }
283
 
                        x[d] += 2.0*dr;
284
 
                        if(mPeriodic[d])
285
 
                        {
286
 
                                if(x[d] >  1.0) x[d] -= 2.0;
287
 
                                if(x[d] < -1.0) x[d] += 2.0;
288
 
                        }
289
 
                        x += 3;
290
 
                }
291
 
        }
292
 
        else
293
 
        {
294
 
                double *x = mDarr1 + 3*kbeg;
295
 
                for(l=kbeg; l<kend; l++)
296
 
                {
297
 
                        if(l%1000 == 0)
298
 
                        {
299
 
                                if(this->IsMaster(p)) mSubject->GetObserver()->SetProgress((float)(l-kbeg)/(kend-kbeg));
300
 
                                if(mSubject->GetObserver()->IsAborted()) return 2;
301
 
                        }
302
 
                        x[d] += 2.0*dr;
303
 
                        if(mPeriodic[d])
304
 
                        {
305
 
                                if(x[d] >  1.0) x[d] -= 2.0;
306
 
                                if(x[d] < -1.0) x[d] += 2.0;
307
 
                        }
308
 
                        x += 3;
309
 
                }
310
 
        }
311
 
 
312
 
        return 0;
 
570
        int d = mDim;
 
571
 
 
572
        iParallel::SplitRange(p,mTot,kbeg,kend,kstp);
 
573
 
 
574
        switch(step)
 
575
        {
 
576
        case 1:
 
577
                {
 
578
                        float dr = mDr;
 
579
 
 
580
                        if(mPaf)
 
581
                        {
 
582
                                float *x = mFarr + 3*kbeg;
 
583
                                for(l=kbeg; l<kend; l++)
 
584
                                {
 
585
                                        if(l%1000 == 0)
 
586
                                        {
 
587
                                                if(this->IsMaster(p)) mSubject->GetObserver()->SetProgress((d+(float)(l-kbeg)/(kend-kbeg))/3.0);
 
588
                                                if(mSubject->GetObserver()->IsAborted()) return 2;
 
589
                                        }
 
590
                                        x[d] += 2.0*dr;
 
591
                                        if(mPeriodic[d])
 
592
                                        {
 
593
                                                if(x[d] >  1.0) x[d] -= 2.0;
 
594
                                                if(x[d] < -1.0) x[d] += 2.0;
 
595
                                        }
 
596
                                        x += 3;
 
597
                                }
 
598
                        }
 
599
                        else
 
600
                        {
 
601
                                double *x = mDarr + 3*kbeg;
 
602
                                for(l=kbeg; l<kend; l++)
 
603
                                {
 
604
                                        if(l%1000 == 0)
 
605
                                        {
 
606
                                                if(this->IsMaster(p)) mSubject->GetObserver()->SetProgress((float)(l-kbeg)/(kend-kbeg));
 
607
                                                if(mSubject->GetObserver()->IsAborted()) return 2;
 
608
                                        }
 
609
                                        x[d] += 2.0*dr;
 
610
                                        if(mPeriodic[d])
 
611
                                        {
 
612
                                                if(x[d] >  1.0) x[d] -= 2.0;
 
613
                                                if(x[d] < -1.0) x[d] += 2.0;
 
614
                                        }
 
615
                                        x += 3;
 
616
                                }
 
617
                        }
 
618
                        return 0;
 
619
                }
 
620
        case 2:
 
621
                {
 
622
                        int j;
 
623
                        float *f = mFarr + d*kbeg, *fmin = mMin + d*p.ThisProc, *fmax = mMax + d*p.ThisProc;
 
624
                        for(j=0; j<d; j++) fmin[j] = fmax[j] = f[j];
 
625
                        for(l=kbeg+1; l<kend; l++)
 
626
                        {
 
627
                                f += d;
 
628
                                for(j=0; j<d; j++)
 
629
                                {
 
630
                                        if(fmin[j] > f[j]) fmin[j] = f[j];
 
631
                                        if(fmax[j] < f[j]) fmax[j] = f[j];
 
632
                                }
 
633
                        }
 
634
                        return 0;
 
635
                }
 
636
        default:
 
637
                {
 
638
                        return 1;
 
639
                }
 
640
        }
313
641
}
314
642