← Back to branch summary

~esys-p-dev/esys-particle/trunk

« back to all changes in this revision

Viewing changes to Tools/StressCalculator/InteractionToStressConverter.cpp

  • Committer: slatham
  • Date: 2004-12-08 07:29:43 UTC
  • Revision ID: svn-v4:1ddc92f8-1f06-0410-a909-b11019f1b28a:lsmtrunk:553
Added copy-constructor and assignment operator for VertexBox. Removed std-deviation parameter and removed redundant computations from the raw2tostress converter.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Tools/StressCalculator/InteractionToStressConverter.cpp
28
28
{
29
29
  namespace lsm
30
30
  {
 
31
    template <typename TmplSphere, typename TmplBox>
 
32
    std::string getDetailsString(const TmplSphere &sphere, const TmplBox &box)
 
33
    {
 
34
      std::stringstream msg;
 
35
      msg
 
36
        << "box = ("
 
37
        << box.getMinPt() << ", " << box.getMaxPt() << ")"
 
38
        << ", sphere = (" << sphere.getCentre() << ", " << sphere.getRadius() << ")";
 
39
      return msg.str();
 
40
    }
 
41
 
 
42
    template <typename TmplSphere, typename TmplBox>
 
43
    void checkIntersectionVolume(double vol, const TmplSphere &sphere, const TmplBox &box)
 
44
    {
 
45
      if (isnan(vol))
 
46
      {
 
47
        std::stringstream msg;
 
48
        msg 
 
49
          << "nan encountered during volume calculation: "
 
50
          << getDetailsString(sphere, box);
 
51
 
 
52
        throw std::runtime_error(msg.str());
 
53
      }
 
54
      else if ((vol < 0.0) && (fabs(vol) > 1.0e-6))
 
55
      {
 
56
        std::stringstream msg;
 
57
        msg 
 
58
          << "Negative intersection volume " << vol 
 
59
          << ". " << getDetailsString(sphere, box);
 
60
        throw std::runtime_error(msg.str());
 
61
      }
 
62
      else if (vol > (box.getVolume() + (box.getVolume()*1.0e-6)))
 
63
      {
 
64
        std::stringstream msg;
 
65
        msg 
 
66
          << "Volume " << vol 
 
67
          << " larger than box volume " << box.getVolume()
 
68
          << ". " << getDetailsString(sphere, box);
 
69
        throw std::runtime_error(msg.str());
 
70
      }
 
71
      else if (vol > (sphere.getVolume() + (sphere.getVolume()*1.0e-6)))
 
72
      {
 
73
        std::stringstream msg;
 
74
        msg 
 
75
          << "Volume " << vol
 
76
          << " larger than sphere volume " << sphere.getVolume()
 
77
          << ". " << getDetailsString(sphere, box);
 
78
        throw std::runtime_error(msg.str());
 
79
      }
 
80
    }
 
81
 
31
82
    class ThreeDIntersectionCalker : public SphereBoxVolCalculator
32
83
    {
33
84
    public:
34
 
      typedef SphereBoxVolCalculator::Sphere Sphere;
35
85
      ThreeDIntersectionCalker(const BoundingBox &box)
36
86
        : SphereBoxVolCalculator(Box(box.getMinPt(), box.getMaxPt()))
37
87
      {
38
88
      }
39
 
      
40
 
      std::string getDetailsString(const Sphere &sphere) const
41
 
      {
42
 
        std::stringstream msg;
43
 
        msg
44
 
          << "box = ("
45
 
          << getBox().getMinPt() << ", " << getBox().getMaxPt() << ")"
46
 
          << ", sphere = (" << sphere.getCentre() << ", " << sphere.getRadius() << ")";
47
 
        return msg.str();
48
 
      }
49
 
      
 
89
 
50
90
      double getVolume(const Sphere &sphere)
51
91
      {
52
92
        const double vol = SphereBoxVolCalculator::getVolume(sphere);
53
 
        if (isnan(vol))
54
 
        {
55
 
          std::stringstream msg;
56
 
          msg 
57
 
            << "nan encountered during volume calculation: " << getDetailsString(sphere);
58
 
          throw std::runtime_error(msg.str());
59
 
        }
60
 
        else if (vol > (getBox().getVolume() + (getBox().getVolume()*1.0e-6)))
61
 
        {
62
 
          std::stringstream msg;
63
 
          msg 
64
 
            << "Volume " << vol 
65
 
            << " larger than box volume " << getBox().getVolume()
66
 
            << ". " << getDetailsString(sphere);
67
 
          throw std::runtime_error(msg.str());
68
 
        }
69
 
        else if (vol > (sphere.getVolume() + (sphere.getVolume()*1.0e-6)))
70
 
        {
71
 
          std::stringstream msg;
72
 
          msg 
73
 
            << "Volume " << vol
74
 
            << " larger than sphere volume " << sphere.getVolume()
75
 
            << ". " << getDetailsString(sphere);
76
 
          throw std::runtime_error(msg.str());
77
 
        }
 
93
        checkIntersectionVolume(vol, sphere, getBox());
78
94
        return vol;
79
95
      }
80
96
    };
82
98
    class TwoDIntersectionCalker  : public CircleBoxVolCalculator
83
99
    {
84
100
    public:
85
 
      typedef CircleBoxVolCalculator::Sphere Sphere;
86
101
      TwoDIntersectionCalker(const BoundingBox &box)
87
102
        :  CircleBoxVolCalculator(Box(box.getMinPt(), box.getMaxPt()))
88
103
      {
89
104
      }
 
105
 
 
106
      double getVolume(const Sphere &sphere)
 
107
      {
 
108
        const double vol = CircleBoxVolCalculator::getVolume(sphere);
 
109
        checkIntersectionVolume(vol, sphere, getBox());
 
110
        return vol;
 
111
      }
90
112
    };
91
113
 
92
114
    class PointDataType
120
142
    typedef vtk::Vec3Type                                    PointType;
121
143
    typedef vtk::UnstructuredPiece<PointType, PointDataType> Piece;
122
144
    typedef EigenvalueCalculator::ComplexRealImagComparer    RealImagComparer;
 
145
    typedef EigenvalueCalculator::ComplexAbsRealImagComparer AbsRealImagComparer;
123
146
    typedef EigenvalueCalculator::ComplexNormComparer        NormComparer;
124
147
 
125
 
    InteractionToStressConverter::InteractionToStressConverter(double stdDeviation)
126
 
      : m_stressCalculator(),
 
148
    InteractionToStressConverter::InteractionToStressConverter(const BoundingBox &bBox, double gridSpacing)
 
149
      : m_gridSpacing(gridSpacing),
 
150
        m_bBox(bBox),
 
151
        m_stressCalculator(),
127
152
        m_stressTensorCollection(m_stressCalculator),
128
 
        m_stdDeviation(stdDeviation)
 
153
        m_regTensorGrid(bBox, gridSpacing),
 
154
        m_regDevStressGrid(bBox, gridSpacing),
 
155
        m_irrStressTensorGrid(bBox, gridSpacing)
129
156
    {
130
157
    }
131
 
    
 
158
 
132
159
    double InteractionToStressConverter::getRealDevStress(const Tensor &stressTensor) const
133
160
    {
134
161
      StressTensor::ComplexVector eigenvals = stressTensor.getEigenvalues();
135
 
      std::sort(eigenvals.begin(), eigenvals.end(), RealImagComparer());
136
 
      
137
162
      double devStress = 0.0;
 
163
 
138
164
      if (is3d())
139
165
      {
 
166
        std::sort(eigenvals.begin(), eigenvals.end(), RealImagComparer());
140
167
        devStress = (eigenvals[2].real() - eigenvals[0].real());
141
168
      }
142
169
      else
143
170
      {
144
 
        if (abs(eigenvals[1].real()) > abs(eigenvals[0].real()))
145
 
        {
146
 
          devStress = (eigenvals[2].real() - eigenvals[1].real());
147
 
        }
148
 
        else
149
 
        {
150
 
          devStress = (eigenvals[2].real() - eigenvals[0].real());
151
 
        }
 
171
        //
 
172
        // In 2d, make sure we only ever subtract \sigma_1 and \sigma_2,
 
173
        // so we sort by absolute value. This ensures that the zero-eigenvalue
 
174
        // (corresponding to the z-dimension) will be sorted into the
 
175
        // eigenvals[0] element.
 
176
        //
 
177
        std::sort(eigenvals.begin(), eigenvals.end(), AbsRealImagComparer());
 
178
        devStress = fabs(eigenvals[2].real() - eigenvals[1].real());
152
179
      }
153
 
      
 
180
 
154
181
      return devStress;
155
182
    }
156
183
 
269
296
    typedef vtk::Vec3Type                                            StrctPointType;
270
297
    typedef vtk::StructuredPiece<StrctPointType, StrctPointDataType> StrctPiece;
271
298
 
272
 
#if 0
273
 
    DoubleGrid InteractionToStressConverter::getDevRegularGrid
274
 
    (
275
 
      const BoundingBox &bBox,
276
 
      double gridSpacing,
277
 
      double stdDeviation
278
 
    )
279
 
    {
280
 
      DoubleGrid irregularGrid = DoubleGrid(bBox, gridSpacing);
 
299
    StressTensorPtrGrid &InteractionToStressConverter::getTensorIrregularGrid()
 
300
    {
 
301
      if (m_irrStressTensorGrid.size() <= 0)
 
302
      {
 
303
        calcTensorIrregularGrid();
 
304
      }
 
305
      return m_irrStressTensorGrid;
 
306
    }
 
307
 
 
308
    double InteractionToStressConverter::getMaxRadius()
 
309
    {
 
310
      StressTensorPtrGrid::ValueIterator it = getTensorIrregularGrid().getValueIterator();
 
311
      double maxRadius = -1.0;
 
312
 
 
313
      while (it.hasNext())
 
314
      {
 
315
        StressTensorPtrGrid::const_reference stressTensor = it.next();
 
316
        if (stressTensor->getRad() > maxRadius)
 
317
        {
 
318
          maxRadius = stressTensor->getRad();
 
319
        }
 
320
      }
 
321
      return maxRadius;
 
322
    }
 
323
 
 
324
    void InteractionToStressConverter::calcTensorIrregularGrid()
 
325
    {
 
326
      m_irrStressTensorGrid = StressTensorPtrGrid(m_bBox, m_gridSpacing);
281
327
      StressTensorCollection::StressTensorIterator it = m_stressTensorCollection.getIterator();
282
 
      
 
328
 
283
329
      while (it.hasNext()) {
284
330
        StressTensorCollection::StressTensorIterator::reference stressTensor = it.next();
285
 
        const double val = getRealDevStress(stressTensor);
286
 
        irregularGrid.insert(stressTensor.getPos(), val);
287
 
      }
288
 
      
289
 
      DoubleGrid regularGrid = DoubleGrid(bBox, gridSpacing);
290
 
      
291
 
      GaussianGridder gridder = GaussianGridder(stdDeviation);
292
 
      gridder.transform(irregularGrid, regularGrid);
293
 
      
294
 
      return regularGrid;
295
 
    }
296
 
#endif
 
331
        m_irrStressTensorGrid.insert(stressTensor.getPos(), &stressTensor);
 
332
      }
 
333
    }
 
334
 
 
335
    TensorGrid &InteractionToStressConverter::getTensorRegularGrid()
 
336
    {
 
337
      if (m_regTensorGrid.size() <= 0)
 
338
      {
 
339
        calcTensorRegularGrid();
 
340
      }
 
341
      return m_regTensorGrid;
 
342
    }
297
343
 
298
344
    template <typename TmplCellIterator, typename TmplIntsectVolCalker>
299
345
    Matrix3 getBoxTensor(
320
366
      return tensor*(1.0/intersectCalker.getBox().getVolume());
321
367
    }
322
368
 
323
 
    TensorGrid InteractionToStressConverter::getTensorRegularGrid
324
 
    (
325
 
      const BoundingBox &bBox,
326
 
      double gridSpacing
327
 
    )
 
369
    void InteractionToStressConverter::calcTensorRegularGrid()
328
370
    {
329
 
      StressTensorPtrGrid irregularGrid = StressTensorPtrGrid(bBox, gridSpacing);
330
 
      StressTensorCollection::StressTensorIterator it = m_stressTensorCollection.getIterator();
331
 
 
332
 
      double maxRadius = -1.0;
333
 
      while (it.hasNext()) {
334
 
        StressTensorCollection::StressTensorIterator::reference stressTensor = it.next();
335
 
        irregularGrid.insert(stressTensor.getPos(), &stressTensor);
336
 
        if (stressTensor.getRad() > maxRadius)
337
 
        {
338
 
          maxRadius = stressTensor.getRad();
339
 
        }
340
 
      }
341
 
 
342
 
      TensorGrid regularGrid = TensorGrid(bBox, gridSpacing);
343
 
 
344
 
      maxRadius += gridSpacing;
345
 
 
346
 
      TensorGrid::CellIterator regCellIt = regularGrid.getCellIterator();
 
371
      StressTensorPtrGrid &irregularGrid = getTensorIrregularGrid();
 
372
 
 
373
      m_regTensorGrid = TensorGrid(m_bBox, m_gridSpacing);
 
374
 
 
375
      const double maxRadius = (getMaxRadius() + m_gridSpacing);
 
376
 
 
377
      TensorGrid::CellIterator regCellIt = m_regTensorGrid.getCellIterator();
347
378
      while (regCellIt.hasNext())
348
379
      {
349
380
        TensorGrid::Cell &regCell = regCellIt.next();
351
382
 
352
383
        if (is3d())
353
384
        {
354
 
          regularGrid.insert(
 
385
          m_regTensorGrid.insert(
355
386
            regCell.getPos(),
356
387
            Tensor(
357
388
              regCell.getPos(),
364
395
        }
365
396
        else
366
397
        {
367
 
          regularGrid.insert(
 
398
          m_regTensorGrid.insert(
368
399
            regCell.getPos(),
369
400
            Tensor(
370
401
              regCell.getPos(),
376
407
          );
377
408
        }
378
409
      }
379
 
 
380
 
      return regularGrid;
381
 
    }
382
 
 
383
 
    DoubleGrid InteractionToStressConverter::getDevRegularGrid
384
 
    (
385
 
      const BoundingBox &bBox,
386
 
      double gridSpacing
387
 
    )
388
 
    {
389
 
      const TensorGrid regularTensorGrid = getTensorRegularGrid(bBox, gridSpacing);
390
 
      DoubleGrid regularGrid = DoubleGrid(bBox, gridSpacing);
 
410
    }
 
411
 
 
412
    DoubleGrid &InteractionToStressConverter::getDevRegularGrid()
 
413
    {
 
414
      if (m_regDevStressGrid.size() <= 0)
 
415
      {
 
416
        calcDevRegularGrid();
 
417
      }
 
418
      return m_regDevStressGrid;
 
419
    }
 
420
 
 
421
    void InteractionToStressConverter::calcDevRegularGrid()
 
422
    {
 
423
      const TensorGrid &regularTensorGrid = getTensorRegularGrid();
 
424
      m_regDevStressGrid = DoubleGrid(m_bBox, m_gridSpacing);
391
425
 
392
426
      TensorGrid::CellConstIterator cellIt = regularTensorGrid.getCellIterator();
393
427
      while (cellIt.hasNext())
396
430
        while (pairIt.hasNext())
397
431
        {
398
432
          const Tensor &tensor = pairIt.next().getValue();
399
 
          regularGrid.insert(
 
433
          m_regDevStressGrid.insert(
400
434
            tensor.getPos(),
401
435
            getRealDevStress(tensor)
402
436
          );
403
437
        }
404
438
      }
405
 
 
406
 
      return regularGrid;
407
439
    }
408
440
 
409
 
    void InteractionToStressConverter::writeVtkStructuredXml
410
 
    (
411
 
      std::ostream &oStream,
412
 
      const BoundingBox &bBox,
413
 
      double gridSpacing
414
 
    )
 
441
    void InteractionToStressConverter::writeVtkStructuredXml(std::ostream &oStream)
415
442
    {
416
 
      TensorGrid regularGrid = getTensorRegularGrid(bBox, gridSpacing);
 
443
      TensorGrid &regularGrid = getTensorRegularGrid();
417
444
      StrctPiece piece(StrctPointType("points"), StrctPointDataType());
418
445
      piece.setExtent(regularGrid.getMinVecIndex(), regularGrid.getMaxVecIndex());
419
446
      TensorGrid::CellIterator cellIt = regularGrid.getCellIterator();
440
467
      grid.writeXml(oStream);
441
468
    }
442
469
    
443
 
    void InteractionToStressConverter::writeFlatStructured
444
 
    (
445
 
      std::ostream &oStream,
446
 
      const BoundingBox &bBox,
447
 
      double gridSpacing
448
 
    )
 
470
    void InteractionToStressConverter::writeFlatStructured(std::ostream &oStream)
449
471
    {
450
 
      DoubleGrid regular = getDevRegularGrid(bBox, gridSpacing);
 
472
      DoubleGrid regular = getDevRegularGrid();
451
473
      DoubleGrid::CellIterator cellIt = regular.getCellIterator();
452
474
      while (cellIt.hasNext())
453
475
      {