Free Electron
Continuum.h
Go to the documentation of this file.
1 /* Copyright (C) 2003-2021 Free Electron Organization
2  Any use of this software requires a license. If a valid license
3  was not distributed with this file, visit freeelectron.org. */
4 
5 /** @file */
6 
7 #ifndef __spatial_Continuum_h__
8 #define __spatial_Continuum_h__
9 
10 #include "SpaceI.h"
11 #include "UnitContinuum.h"
12 
13 namespace fe
14 {
15 namespace ext
16 {
17 
18 template<typename T>
19 class Continuum
20 {
21  public:
22  Continuum(void)
23  {
24  }
25 virtual ~Continuum(void)
26  {
27  }
28 
29  Continuum<T> &operator=(const Continuum<T> &a_other)
30  {
31  m_unitContinuum = a_other.m_unitContinuum;
32  m_space = a_other.m_space;
33  return *this;
34  }
35 
36  bool operator==(const Continuum<T> &a_other)
37  {
38  if(m_unitContinuum != a_other.m_unitContinuum)
39  {
40  return false;
41  }
42  if(m_space != a_other.m_space)
43  {
44  return false;
45  }
46  return true;
47  }
48 
49  typedef typename UnitContinuum<T,3>::t_index t_index;
50  typedef typename UnitContinuum<T,3>::t_cells t_cells;
51  typedef SpaceI t_space;
52 
53  void create( const t_index &a_count,
54  const T &a_init,
55  const sp<t_space> &a_space);
56 
57 
58  // unsafe location to index methods (indices may be out of range)
59  void index( t_index &a_index,
60  const SpatialVector &a_location) const;
61  void lower( t_index &a_index,
62  const SpatialVector &a_location) const;
63 
64  // index to location methods
65  void location( SpatialVector &a_location,
66  const t_index &a_index) const;
67 
68  // out-of-range safety
69  SpatialVector &safe( SpatialVector &a_location) const;
70  bool isSafe( const SpatialVector &a_location) const;
71 
72  // unit continuum access
73  const UnitContinuum<T,3> &unit(void) const;
74  UnitContinuum<T,3> &accessUnit(void);
75 
76  sp<t_space> &space(void) { return m_space; }
77 
78  void cellsize(SpatialVector &a_cellsize);
79 
80 
81  private:
82  UnitContinuum<T,3> m_unitContinuum;
83  sp<t_space> m_space;
84 };
85 
86 template<typename T>
87 inline void Continuum<T>::create(const t_index &a_count, const T &a_init,
88  const sp<t_space> &a_space)
89 {
90  m_unitContinuum.create(a_count, a_init);
91  m_space = a_space;
92 }
93 
94 
95 template<typename T>
96 inline void Continuum<T>::index(t_index &a_index,
97  const SpatialVector &a_location) const
98 {
99  SpatialVector unit_location;
100  m_space->from(unit_location, a_location);
101  typename UnitContinuum<T,3>::t_span location;
102  location = unit_location;
103  m_unitContinuum.index(a_index, location);
104 }
105 
106 template<typename T>
107 inline void Continuum<T>::lower(t_index &a_index,
108  const SpatialVector &a_location) const
109 {
110  SpatialVector unit_location;
111  m_space->from(unit_location, a_location);
112  m_unitContinuum.lower(a_index, unit_location);
113 }
114 
115 template<typename T>
116 inline void Continuum<T>::location(SpatialVector &a_location,
117  const t_index &a_index) const
118 {
119  typename UnitContinuum<T,3>::t_span unit_location;
120  m_unitContinuum.location(unit_location, a_index);
121  SpatialVector location = unit_location;
122  m_space->to(a_location, location);
123 }
124 
125 template<typename T>
126 inline SpatialVector &Continuum<T>::safe(SpatialVector &a_span) const
127 {
128  SpatialVector unit_span;
129  m_space->from(unit_span, a_span);
130  typename UnitContinuum<T,3>::t_span uspan;
131  uspan = unit_span;
132  m_unitContinuum.safe(uspan);
133  unit_span = uspan;
134  m_space->to(a_span, unit_span);
135  return a_span;
136 }
137 
138 template<typename T>
139 inline bool Continuum<T>::isSafe(const SpatialVector &a_span) const
140 {
141  SpatialVector unit_span;
142  m_space->from(unit_span, a_span);
143  typename UnitContinuum<T,3>::t_index index, safe_index;
144  typename UnitContinuum<T,3>::t_span uspan;
145  uspan = unit_span;
146  m_unitContinuum.index(index, uspan);
147  safe_index = index;
148  m_unitContinuum.safe(safe_index);
149  return (index == safe_index);
150 }
151 
152 template<typename T>
153 inline const UnitContinuum<T,3> &Continuum<T>::unit(void) const
154 {
155  return m_unitContinuum;
156 }
157 
158 template<typename T>
159 inline UnitContinuum<T,3> &Continuum<T>::accessUnit(void)
160 {
161  return m_unitContinuum;
162 }
163 
164 template<typename T>
165 inline void Continuum<T>::cellsize(SpatialVector &a_cellsize)
166 {
167  typename UnitContinuum<T,3>::t_span unitsz = unit().cellsize();
168  SpatialVector origin(0.0,0.0,0.0);
169  SpatialVector sz;
170  sz = unitsz;
171  m_space->to(origin, origin);
172  m_space->to(a_cellsize, sz);
173  a_cellsize -= origin;
174 }
175 
176 // ============================================================================
177 // ============================================================================
178 
179 
180 /** sample a Continuum by returning the direct cell value
181  @relates Continuum */
182 template<typename T>
183 const T &directSample(const Continuum<T> &a_continuum,
184  const SpatialVector &a_location)
185 {
186  typename Continuum<T>::t_index idx;
187  a_continuum.index(idx, a_location);
188  return a_continuum.unit()[idx];
189 }
190 
191 
192 /** Calculate the gradient of a UnitContinuum at a point using midpoint
193  differencing and linear interpolation.
194  @relates Continuum */
195 void gradient(Continuum<Real> &a_continuum, SpatialVector &a_gradient,
196  const SpatialVector &a_location, Real a_h = 0.1);
197 
198 struct t_n_radial_set
199 {
200  SpatialVector *m_location;
201  SpatialVector *m_radius;
202  Continuum<Real> *m_continuum;
203  Real m_value;
204 };
205 
206 void rad_func(UnitContinuum<Real,3> &a_unit,
207  UnitContinuum<Real,3>::t_index &a_index, t_n_radial_set *a_data);
208 
209 /** Radial set function.
210  @relates Continuum */
211 FE_DL_EXPORT void radialSet(Continuum<Real> &a_continuum,
212  const SpatialVector &a_location,
213  const SpatialVector &a_radius, Real a_value);
214 
215 /** Calculate a gradient vector field for a scalar field.
216  @relates Continuum */
217 template<typename T>
218 void derive(Continuum<T> &a_vector, Continuum<Real> &a_scalar)
219 {
220  derive<T,3>(a_vector.accessUnit(), a_scalar.accessUnit());
221  a_vector.space() = a_scalar.space();
222 }
223 
224 
225 // ============================================================================
226 // old implementation. left here for a while for easy reference
227 // ============================================================================
228 #if 0
229 
230 /** N-Dimensional Spatial Continuum using a evenly spaced grid descretization
231  */
232 template<typename T, unsigned int D>
233 class Continuum
234 {
235  public:
236  Continuum(void)
237  {
238  }
239 virtual ~Continuum(void)
240  {
241  }
242 
243  Continuum<T,D> &operator=(const Continuum<T,D> &other)
244  {
245  m_origin = other.m_origin;
246  m_count = other.m_count;
247  for(unsigned int i = 0; i < D; i++)
248  {
249  m_boundary[i] = other.m_boundary[i];
250  m_size[i] = other.m_size[i];
251  }
252  m_cells = other.m_cells;
253  m_totalCount = other.m_totalCount;
254  return *this;
255  }
256 
257  bool operator==(const Continuum<T,D> &other)
258  {
259  if(m_origin != other.m_origin) { return false; }
260  if(m_count != other.m_count) { return false; }
261  if(m_boundary != other.m_boundary) { return false; }
262  if(m_size != other.m_size) { return false; }
263  if(m_totalCount != other.m_totalCount) { return false; }
264  if(m_cells == other.m_cells) { return true; }
265  return false;
266  }
267 
268  template <typename P>
269  void copyTopology(const Continuum<P,D> &a_other,
270  const T &a_init)
271  {
272  create(a_other.m_count, a_other.m_boundary, a_init);
273  }
274 
275  typedef Vector<D,int> t_index;
276  typedef Vector<D,Real> t_span;
277  typedef std::vector<T> t_cells;
278 
279  void create( const t_index &a_count,
280  const t_span &a_boundary,
281  const T &a_init);
282 
283  void create( const t_span &a_origin,
284  const t_index &a_count,
285  const t_span &a_boundary,
286  const T &a_init);
287 
288  const T &operator[](const t_index &a_index) const;
289 
290  void set(const t_index &a_index, const T &value);
291  void set(const T &value);
292  T &access(const t_index &a_index);
293 
294  void direct( t_index &a_direct,
295  const t_span &a_location) const;
296  void closest( t_index &a_closest,
297  const t_span &a_location) const;
298  void space( t_span &a_space,
299  const t_index &a_index) const;
300  const t_span &cellsize(void) const { return m_size; }
301  const t_index &count(void) const { return m_count; }
302  const t_span &origin(void) const { return m_origin; }
303  const t_span &boundary(void) const { return m_boundary; }
304  unsigned int totalCount(void) const { return m_totalCount; }
305  //t_span locate(const t_index &a_index) const;
306 
307  template<typename A>
308  void iterate( void (*a_function)(Continuum<T,D> &,
309  t_index &,
310  A),
311  A a_userarg);
312 
313  void addNeighbors(
314  std::list<t_index> &a_neighbors,
315  const t_index &a_index) const;
316  void addNeighborsDiagonal(
317  std::list<t_index> &a_neighbors,
318  const t_index &a_index) const;
319  t_index &safe( t_index &a_index) const;
320  t_span &safe( t_span &a_span) const;
321  private:
322  bool safe( unsigned int a_i,
323  t_index &a_index) const;
324  void addNeighborsD(unsigned int a_d,
325  std::list<t_index> &a_neighbors,
326  const t_index &a_index,
327  const t_index &a_base) const;
328  unsigned int calc_array_index( const t_index &a_index) const;
329  template<typename A>
330  void rIterate( unsigned int a_d,
331  t_index &a_index,
332  void (*a_function)(Continuum<T,D> &,
333  t_index &,
334  A),
335  A a_userarg);
336 
337  private:
338  t_span m_origin;
339  t_index m_count;
340  t_span m_boundary;
341  t_span m_size;
342  t_cells m_cells;
343  unsigned int m_totalCount;
344 };
345 
346 template<typename T, unsigned int D>
347 inline bool Continuum<T,D>::safe(unsigned int a_i, t_index &a_index) const
348 {
349  bool rv = true;
350  if(m_count[a_i] == 1)
351  {
352  if(a_index[a_i] != 0)
353  {
354  if(m_boundary[a_i] >= 0.0)
355  {
356  rv = false;
357  }
358  a_index[a_i] = 0;
359  }
360  }
361  else if(m_count[a_i] == 2)
362  {
363  if(a_index[a_i] < 0)
364  {
365  if(m_boundary[a_i] >= 0.0)
366  {
367  rv = false;
368  }
369  a_index[a_i] = 0;
370  }
371  else if(a_index[a_i] > 1)
372  {
373  if(m_boundary[a_i] >= 0.0)
374  {
375  rv = false;
376  }
377  a_index[a_i] = 1;
378  }
379  }
380  else
381  {
382  while(a_index[a_i] < 0)
383  {
384  if(m_boundary[a_i] < 0.0)
385  {
386  a_index[a_i] += m_count[a_i] - 1;
387  }
388  else
389  {
390  rv = false;
391  a_index[a_i] = 0;
392  }
393  }
394  while(a_index[a_i] >= m_count[a_i])
395  {
396  if(m_boundary[a_i] < 0.0)
397  {
398  a_index[a_i] -= m_count[a_i] - 1;
399  }
400  else
401  {
402  rv = false;
403  a_index[a_i] = m_count[a_i]-1;
404  }
405  }
406  }
407 
408  return rv;
409 }
410 
411 template<typename T, unsigned int D>
412 inline typename Continuum<T,D>::t_index &Continuum<T,D>::safe(t_index &a_index) const
413 {
414  for(unsigned int i = 0; i < D; i++)
415  {
416  if(m_count[i] == 1)
417  {
418  a_index[i] = 0;
419  }
420  else if(m_count[i] == 2)
421  {
422  if(a_index[i] < 0)
423  {
424  a_index[i] = 0;
425  }
426  else if(a_index[i] > 1)
427  {
428  a_index[i] = 1;
429  }
430  }
431  else
432  {
433  while(a_index[i] < 0)
434  {
435  if(m_boundary[i] < 0.0)
436  {
437  a_index[i] += m_count[i] - 1;
438  }
439  else
440  {
441  a_index[i] = 0;
442  }
443  }
444  while(a_index[i] >= m_count[i])
445  {
446  if(m_boundary[i] < 0.0)
447  {
448  a_index[i] -= m_count[i] - 1;
449  }
450  else
451  {
452  a_index[i] = m_count[i]-1;
453  }
454  }
455  }
456  }
457 
458  return a_index;
459 }
460 
461 template<typename T, unsigned int D>
462 inline typename Continuum<T,D>::t_span &Continuum<T,D>::safe(
463  t_span &a_span) const
464 {
465  a_span -= m_origin;
466  for(unsigned int i = 0; i < D; i++)
467  {
468  if(m_count[i] == 1)
469  {
470  a_span[i] = 0.0;
471  }
472  else
473  {
474  Real bnd = fabs(m_boundary[i]);
475  while(a_span[i] < 0.0)
476  {
477  if(m_boundary[i] < 0.0)
478  {
479  a_span[i] += bnd;
480  }
481  else
482  {
483  a_span[i] = 0.0;
484  }
485  }
486  while(a_span[i] > bnd)
487  {
488  if(m_boundary[i] < 0.0)
489  {
490  a_span[i] -= bnd;
491  }
492  else
493  {
494  a_span[i] = bnd;
495  }
496  }
497  }
498  }
499  a_span += m_origin;
500 
501  return a_span;
502 }
503 
504 template<typename T, unsigned int D>
505 template<typename A>
506 inline void Continuum<T,D>::iterate(void (*a_function)(Continuum<T,D> &,
507  t_index &, A), A a_userarg)
508 {
509  t_index index;
510  rIterate(D-1, index, a_function, a_userarg);
511 }
512 
513 template<typename T, unsigned int D>
514 template<typename A>
515 inline void Continuum<T,D>::rIterate(unsigned int a_d, t_index &a_index,
516  void (*a_function)(Continuum<T,D> &, t_index &, A), A a_userarg)
517 {
518  for(int i = 0; i < count()[a_d]; i++)
519  {
520  a_index[a_d] = i;
521  if(a_d == 0)
522  {
523  a_function(*this, a_index, a_userarg);
524  }
525  else
526  {
527  rIterate(a_d-1, a_index, a_function, a_userarg);
528  }
529  }
530 }
531 
532 
533 template<typename T, unsigned int D>
534 inline void Continuum<T,D>::direct(t_index &a_index,
535  const t_span &a_location) const
536 {
537  t_span location = a_location - m_origin;
538  for(int i = D-1; i >= 0; i--)
539  {
540  if(m_size[i] > 0.0)
541  {
542  a_index[i] = (int)(location[i] / m_size[i]);
543  }
544  else
545  {
546  a_index[i] = 0;
547  continue;
548  }
549  if(a_index[i] >= m_count[i] - 1)
550  {
551  a_index[i] = m_count[i] - 2;
552  }
553  if(a_index[i] < 0)
554  {
555  a_index[i] = 0;
556  }
557  }
558 }
559 
560 template<typename T, unsigned int D>
561 inline void Continuum<T,D>::closest(t_index &a_index,
562  const t_span &a_location) const
563 {
564  t_span location = a_location - m_origin;
565  for(int i = D-1; i >= 0; i--)
566  {
567  if(m_size[i] > 0.0)
568  {
569  a_index[i] = (int)((location[i]+(0.5*m_size[i])) / m_size[i]);
570  }
571  else
572  {
573  a_index[i] = 0;
574  continue;
575  }
576  if(a_index[i] >= m_count[i])
577  {
578  a_index[i] = m_count[i] - 1;
579  }
580  if(a_index[i] < 0)
581  {
582  a_index[i] = 0;
583  }
584  }
585 }
586 
587 template<typename T, unsigned int D>
588 inline unsigned int Continuum<T,D>::calc_array_index(
589  const t_index &a_index) const
590 {
591  unsigned int index = 0;
592  unsigned int step = 1;
593  for(int d = D-1; d >= 0; d--)
594  {
595  index += a_index[d] * step;
596  step *= m_count[d];
597  }
598  return index;
599 }
600 
601 
602 template<typename T, unsigned int D>
603 inline void Continuum<T,D>::space(t_span &a_space, const t_index &a_index) const
604 {
605  for(unsigned int i = 0; i < D; i++)
606  {
607  a_space[i] = (Real)a_index[i] * cellsize()[i] + m_origin[i];
608  }
609 }
610 
611 template<typename T, unsigned int D>
612 inline void Continuum<T,D>::create(const t_index &a_count,
613  const t_span &a_boundary, const T &a_init)
614 {
615  m_totalCount = 1;
616  for(unsigned int i = 0; i < D; i++)
617  {
618  m_boundary[i] = a_boundary[i];
619  m_count[i] = a_count[i];
620  Real bnd = fabs(m_boundary[i]);
621  if(m_count[i] == 1)
622  {
623  m_size[i] = 0.0;
624  }
625  else
626  {
627  m_size[i] = bnd / (m_count[i] - 1);
628  }
629  m_totalCount *= m_count[i];
630  }
631 
632  m_cells.resize(m_totalCount, a_init);
633 }
634 
635 template<typename T, unsigned int D>
636 inline void Continuum<T,D>::create(const t_span &a_origin,
637  const t_index &a_count,
638  const t_span &a_boundary, const T &a_init)
639 {
640  m_origin = a_origin;
641  create(a_count, a_boundary, a_init);
642 }
643 
644 
645 template<typename T, unsigned int D>
646 inline const T &Continuum<T,D>::operator[](const t_index &a_index) const
647 {
648  return m_cells[calc_array_index(a_index)];
649 }
650 
651 template<typename T, unsigned int D>
652 inline T &Continuum<T,D>::access(const t_index &a_index)
653 {
654  return m_cells[calc_array_index(a_index)];
655 }
656 
657 template<typename T, unsigned int D>
658 inline void Continuum<T,D>::set(const t_index &a_index, const T &a_value)
659 {
660  m_cells[calc_array_index(a_index)] = a_value;
661 }
662 
663 template<typename T, unsigned int D>
664 inline void Continuum<T,D>::set(const T &a_value)
665 {
666  for(unsigned int i = 0; i < m_cells.size(); i++)
667  {
668  m_cells[i] = a_value;
669  }
670 }
671 
672 template<typename T, unsigned int D>
673 inline void Continuum<T,D>::addNeighborsD(unsigned int a_d,
674  std::list<t_index> &a_neighbors, const t_index &a_index,
675  const t_index &a_base) const
676 {
677  t_index candidate = a_index;
678  if(a_d == 0)
679  {
680  candidate[a_d] = a_index[a_d] + 1;
681  if(safe(a_d, candidate))
682  {
683  if(!(candidate == a_base))
684  {
685  a_neighbors.push_back(candidate);
686  }
687  }
688 
689  candidate[a_d] = a_index[a_d] - 1;
690  if(safe(a_d, candidate))
691  {
692  if(!(candidate == a_base))
693  {
694  a_neighbors.push_back(candidate);
695  }
696  }
697 
698  candidate[a_d] = a_index[a_d];
699  if(safe(a_d, candidate))
700  {
701  if(!(candidate == a_base))
702  {
703  a_neighbors.push_back(candidate);
704  }
705  }
706  }
707  else
708  {
709  candidate[a_d] = a_index[a_d] + 1;
710  safe(a_d, candidate);
711  if(candidate[a_d] != a_index[a_d])
712  {
713  addNeighborsD(a_d-1, a_neighbors, candidate, a_base);
714  }
715 
716  candidate[a_d] = a_index[a_d] - 1;
717  safe(a_d, candidate);
718  if(candidate[a_d] != a_index[a_d])
719  {
720  addNeighborsD(a_d-1, a_neighbors, candidate, a_base);
721  }
722 
723  candidate[a_d] = a_index[a_d];
724  safe(a_d, candidate);
725  addNeighborsD(a_d-1, a_neighbors, candidate, a_base);
726 
727 
728  }
729 }
730 
731 template<typename T, unsigned int D>
732 inline void Continuum<T,D>::addNeighborsDiagonal(
733  std::list<t_index> &a_neighbors, const t_index &a_index) const
734 {
735  t_index candidate = a_index;
736  safe(candidate);
737  if(!(candidate == a_index))
738  {
739  feX("Continuum<T,D>::addNeighbors",
740  "invalid reference index");
741  }
742  addNeighborsD(D-1, a_neighbors, candidate, a_index);
743 }
744 
745 template<typename T, unsigned int D>
746 inline void Continuum<T,D>::addNeighbors(std::list<t_index> &a_neighbors,
747  const t_index &a_index) const
748 {
749  t_index candidate = a_index;
750  safe(candidate);
751  if(!(candidate == a_index))
752  {
753  feX("Continuum<T,D>::addNeighbors",
754  "invalid reference index");
755  }
756  for(unsigned int i = 0; i < D; i++)
757  {
758  candidate[i] = a_index[i] + 1;
759  safe(i, candidate);
760  if(candidate[i] != a_index[i])
761  {
762  a_neighbors.push_back(candidate);
763  }
764 
765  candidate[i] = a_index[i] - 1;
766  safe(i, candidate);
767  if(candidate[i] != a_index[i])
768  {
769  a_neighbors.push_back(candidate);
770  }
771 
772  candidate[i] = a_index[i];
773  }
774 }
775 
776 #if 0
777 template<typename T, unsigned int D>
778 inline typename Continuum<T,D>::t_span Continuum<T,D>::locate(
779  const t_index &a_index) const
780 {
781  typename Continuum<T,D>::t_span location;
782  space(location, a_index);
783  return location;
784 }
785 #endif
786 
787 
788 
789 // ============================================================================
790 // ============================================================================
791 
792 
793 
794 /** sample a Continuum by returning the value at the nearest grid point with all
795  dimension values lower than the sample point (basicly, the cell origin)
796  @relates Continuum */
797 template<typename T, unsigned int D>
798 const T &directSample(const Continuum<T,D> &a_continuum,
799  const typename Continuum<T,D>::t_span &a_location)
800 {
801  typename Continuum<T,D>::t_index direct;
802  a_continuum.direct(direct, a_location);
803  return a_continuum[direct];
804 }
805 
806 /** sample a Continuum by returning the value at the nearest grid point.
807  @relates Continuum */
808 template<typename T, unsigned int D>
809 const T &closestSample(const Continuum<T,D> &a_continuum,
810  const typename Continuum<T,D>::t_span &a_location)
811 {
812  typename Continuum<T,D>::t_index closest;
813  a_continuum.closest(closest, a_location);
814  a_continuum[closest];
815  return a_continuum[closest];
816 }
817 
818 struct d_bnd
819 {
820  Real m_lo;
821  Real m_hi;
822 };
823 
824 template<unsigned int D>
825 Real rLinearSample(unsigned int a_d, const Continuum<Real,3> &a_continuum,
826  typename Continuum<Real,D>::t_index &a_index,
827  typename Continuum<Real,D>::t_index &a_direct, float a_mult,
828  d_bnd *a_bnd)
829 {
830  Real accum = 0.0;
831  Real m;
832  a_index[a_d] = a_direct[a_d];
833  a_continuum.safe(a_index);
834  m = a_bnd[a_d].m_hi * a_mult;
835  if(m > 0.0)
836  {
837  if(a_d == 0)
838  {
839  accum += m * a_continuum[a_index];
840  }
841  else
842  {
843  accum += rLinearSample<D>(a_d-1, a_continuum, a_index, a_direct,
844  m, a_bnd);
845  }
846  }
847 
848  a_index[a_d] = a_direct[a_d] + 1;
849  a_continuum.safe(a_index);
850  m = a_bnd[a_d].m_lo * a_mult;
851  if(m > 0.0)
852  {
853  if(a_d == 0)
854  {
855  accum += m * a_continuum[a_index];
856  }
857  else
858  {
859  accum += rLinearSample<D>(a_d-1, a_continuum, a_index, a_direct,
860  m, a_bnd);
861  }
862  }
863 
864  return accum;
865 }
866 
867 
868 
869 /** sample a Continuum by returning the a linear interpolated value.
870  @relates Continuum */
871 template<unsigned int D>
872 Real linearSample(const Continuum<Real,D> &a_continuum,
873  const typename Continuum<Real,D>::t_span &a_location)
874 {
875  d_bnd bnd[D];
876 
877  typename Continuum<Real,D>::t_index direct;
878  a_continuum.direct(direct, a_location);
879 
880  typename Continuum<Real,D>::t_span location =
881  a_location - a_continuum.origin();
882 
883  for(unsigned int i = 0; i < D; i++)
884  {
885  if(a_continuum.cellsize()[i] == 0.0)
886  {
887  bnd[i].m_lo = 0.0;
888  bnd[i].m_hi = 1.0;
889  }
890  else
891  {
892  bnd[i].m_lo = location[i]/a_continuum.cellsize()[i] -
893  (Real)direct[i];
894 
895  bnd[i].m_hi = 1.0 - bnd[i].m_lo;
896  }
897  }
898 
899  Real r = 1.0;
900 
901  typename Continuum<Real,D>::t_index index;
902  return rLinearSample<D>(D-1, a_continuum, index, direct, r, bnd);
903 }
904 
905 
906 template<unsigned int D>
907 Real rPolySample( unsigned int a_d,
908  const Continuum<Real,3> &a_continuum,
909  const typename Continuum<Real,D>::t_span &a_location,
910  typename Continuum<Real,D>::t_index &a_index,
911  typename Continuum<Real,D>::t_index &a_direct,
912  const typename Continuum<Real,D>::t_index &a_block)
913 {
914 
915  int size = a_block[a_d];
916  int minus = (size / 2) - 1;
917  if(minus < 0) { minus = 0; }
918  int pre = 0;
919  int post = 0;
920  if((int)a_direct[a_d] < minus)
921  {
922  pre = minus - a_direct[a_d];
923  size -= pre;
924  minus = a_direct[a_d];
925  }
926  int root = (int)a_direct[a_d] - minus;
927  if(root + size > (int)a_continuum.count()[a_d])
928  {
929  post = size - (a_continuum.count()[a_d] - root);
930  size = a_continuum.count()[a_d] - root;
931  }
932 
933  if(size < 1) { size = 1; }
934 
935  std::vector<Vector<2,Real> > xy(size+pre+post);
936  for(int i = 0; i < size; i++)
937  {
938  a_index[a_d] = root + i;
939 
940  xy[i+pre][0] = a_continuum.cellsize()[a_d] * (Real)(root + i);
941  if(a_d == 0)
942  {
943  xy[i+pre][1] = a_continuum[a_index];
944  }
945  else
946  {
947  xy[i+pre][1] = rPolySample<D>( a_d-1,
948  a_continuum,
949  a_location,
950  a_index,
951  a_direct,
952  a_block);
953  }
954  }
955  for(int i = 0; i < pre; i++)
956  {
957  xy[i][0] = a_continuum.cellsize()[a_d] * (Real)(i - pre);
958  xy[i][1] = xy[pre][1];
959  }
960  for(int i = 0; i < post; i++)
961  {
962  xy[size+pre+i][0] = a_continuum.cellsize()[a_d] * (Real)(i + size + root);
963  xy[size+pre+i][1] = xy[size+pre-1][1];
964  }
965 
966  Real y,dy;
967  polyInterp<Real>(xy,a_location[a_d],y,dy);
968  return y;
969 }
970 
971 /** sample a Continuum by returning the a polynomial interpolated value.
972  @relates Continuum */
973 template<unsigned int D>
974 Real polySample( const Continuum<Real,D> &a_continuum,
975  const typename Continuum<Real,D>::t_span &a_location,
976  const typename Continuum<Real,D>::t_index &a_block)
977 {
978  typename Continuum<Real,D>::t_index index;
979  typename Continuum<Real,D>::t_index direct;
980  a_continuum.direct(direct, a_location);
981  typename Continuum<Real,D>::t_span location =
982  a_location - a_continuum.origin();
983  return rPolySample<D>( D-1,
984  a_continuum,
985  location,
986  index,
987  direct,
988  a_block);
989 }
990 
991 /** sample a Continuum by returning the a polynomial interpolated value.
992  @relates Continuum */
993 template<unsigned int D>
994 Real polySample( const Continuum<Real,D> &a_continuum,
995  const typename Continuum<Real,D>::t_span &a_location)
996 {
997  typename Continuum<Real,D>::t_index block;
998  for(unsigned int i = 0; i < D; i++)
999  {
1000  if(a_continuum.count()[i] < 4)
1001  {
1002  block[i] = a_continuum.count()[i];
1003  }
1004  else
1005  {
1006  block[i] = 4;
1007  }
1008  }
1009  return polySample(a_continuum, a_location, block);
1010 }
1011 
1012 template<unsigned int D>
1013 struct t_radial_set
1014 {
1015  typename Continuum<Real,D>::t_span *m_location;
1016  typename Continuum<Real,D>::t_span *m_radius;
1017  Real m_value;
1018 };
1019 
1020 template<unsigned int D>
1021 void rad_func(Continuum<Real,D> &a_continuum, typename Continuum<Real,D>::t_index &a_index, t_radial_set<D> *a_data)
1022 {
1023  Real value = 0.0;
1024  typename Continuum<Real,D>::t_span index_loc;
1025  a_continuum.space(index_loc, a_index);
1026  for(unsigned int i = 0; i < D; i++)
1027  {
1028  Real dist = fabs(index_loc[i] - (*(a_data->m_location))[i]);
1029  if(a_continuum.boundary()[i] < 0.0)
1030  {
1031  if(dist > fabs(a_continuum.boundary()[i]) / 2.0)
1032  {
1033  dist = fabs(a_continuum.boundary()[i]) - dist;
1034  }
1035  }
1036  dist /= (*(a_data->m_radius))[i];
1037  dist *= dist;
1038  if(dist > 1.0) return;
1039  value += dist;
1040  }
1041  value = 1.0-sqrt(value);
1042  if(value < 0.0) return;
1043  value = a_data->m_value*value;
1044  a_continuum.set(a_index, value);
1045 }
1046 
1047 /** Radial set function.
1048  @relates Continuum */
1049 template<unsigned int D>
1050 void radialSet(Continuum<Real,D> &a_continuum,
1051  const typename Continuum<Real,D>::t_span &a_location,
1052  const typename Continuum<Real,D>::t_span &a_radius, Real a_value)
1053 {
1054  t_radial_set<D> radial_data;
1055  radial_data.m_location =
1056  const_cast<typename Continuum<Real,D>::t_span *>(&a_location);
1057  radial_data.m_radius =
1058  const_cast<typename Continuum<Real,D>::t_span *>(&a_radius);
1059  radial_data.m_value = a_value;
1060  a_continuum.iterate(rad_func<D>, &radial_data);
1061 }
1062 
1063 
1064 /** Calculate the gradient of a Continuum at a point using midpoint
1065  differencing and polynomial interpolation.
1066  @relates Continuum */
1067 template<unsigned int D>
1068 void gradientMidPoly(Continuum<Real,D> &a_continuum,
1069  typename Continuum<Real,D>::t_span &a_gradient,
1070  const typename Continuum<Real,D>::t_span &a_location, Real a_h = 0.1)
1071 {
1072  for(unsigned int i = 0; i < D; i++)
1073  {
1074  typename Continuum<Real,D>::t_span lo(a_location);
1075  typename Continuum<Real,D>::t_span hi(a_location);
1076  lo[i] -= a_h;
1077  hi[i] += a_h;
1078  a_continuum.safe(lo);
1079  a_continuum.safe(hi);
1080  Real lo_smpl = polySample(a_continuum, lo);
1081  Real hi_smpl = polySample(a_continuum, hi);
1082  a_gradient[i] = (hi_smpl - lo_smpl)/(2.0*a_h);
1083  }
1084 }
1085 
1086 /** Calculate the gradient of a Continuum at a point using forward
1087  differencing and linear interpolation.
1088  @relates Continuum */
1089 template<unsigned int D>
1090 void gradientFwdLinear(Continuum<Real,D> &a_continuum,
1091  typename Continuum<Real,D>::t_span &a_gradient,
1092  const typename Continuum<Real,D>::t_span &a_location, Real a_h = 0.1)
1093 {
1094  for(unsigned int i = 0; i < D; i++)
1095  {
1096  typename Continuum<Real,D>::t_span lo(a_location);
1097  typename Continuum<Real,D>::t_span hi(a_location);
1098  hi[i] += a_h;
1099  a_continuum.safe(lo);
1100  a_continuum.safe(hi);
1101  Real lo_smpl = linearSample(a_continuum, lo);
1102  Real hi_smpl = linearSample(a_continuum, hi);
1103  a_gradient[i] = (hi_smpl - lo_smpl)/(a_h);
1104  }
1105 }
1106 
1107 /** Calculate the gradient of a Continuum at a point using midpoint
1108  differencing and linear interpolation.
1109  @relates Continuum */
1110 template<unsigned int D>
1111 void gradient(Continuum<Real,D> &a_continuum,
1112  typename Continuum<Real,D>::t_span &a_gradient,
1113  const typename Continuum<Real,D>::t_span &a_location, Real a_h = 0.1)
1114 {
1115  for(unsigned int i = 0; i < D; i++)
1116  {
1117  typename Continuum<Real,D>::t_span lo(a_location);
1118  typename Continuum<Real,D>::t_span hi(a_location);
1119  lo[i] -= a_h;
1120  hi[i] += a_h;
1121  a_continuum.safe(lo);
1122  a_continuum.safe(hi);
1123  Real lo_smpl = linearSample(a_continuum, lo);
1124  Real hi_smpl = linearSample(a_continuum, hi);
1125  a_gradient[i] = (hi_smpl - lo_smpl)/(2.0*a_h);
1126  }
1127 }
1128 
1129 
1130 template<typename T, unsigned int D>
1131 struct t_derive
1132 {
1133  Continuum<T,D> *m_vectors;
1134 };
1135 
1136 template<typename T, unsigned int D>
1137 void derive_fun(Continuum<Real,D> &a_continuum, typename Continuum<Real,D>::t_index &a_index, t_derive<T,D> *a_data)
1138 {
1139  for(unsigned int d = 0; d < D; d++)
1140  {
1141  if(a_continuum.cellsize()[d] > 0.0)
1142  {
1143  typename Continuum<Real,D>::t_index dn(a_index);
1144  dn[d] -= 1;
1145  a_continuum.safe(dn);
1146  typename Continuum<Real,D>::t_index up(a_index);
1147  up[d] += 1;
1148  a_continuum.safe(up);
1149  Real value = (a_continuum[up] - a_continuum[dn]) /
1150  (2.0*a_continuum.cellsize()[d]);
1151  setAt(a_data->m_vectors->access(a_index), d, value);
1152  }
1153  else
1154  {
1155  setAt(a_data->m_vectors->access(a_index), d, 0.0f);
1156  }
1157  }
1158 }
1159 
1160 /** Calculate a gradient vector field for a scalar field.
1161  @relates Continuum */
1162 template<typename T, unsigned int D>
1163 void derive(Continuum<T,D> &a_vector, Continuum<Real,D> &a_scalar)
1164 {
1165  T dummy;
1166  a_vector.create(a_scalar.count(), a_scalar.boundary(), dummy);
1167  t_derive<T,D> derive_data;
1168  derive_data.m_vectors = &a_vector;
1169  a_scalar.iterate(derive_fun<T,D>, &derive_data);
1170 }
1171 #endif
1172 
1173 
1174 } /* namespace ext */
1175 } /* namespace fe */
1176 
1177 #endif /* __spatial_Continuum_h__ */
1178 
fe
kernel
Definition: namespace.dox:3
fe::operator==
BWORD operator==(const DualString &s1, const DualString &s2)
Compare two DualString&#39;s.
Definition: DualString.h:208
Generated by   doxygen 1.8.13