Free Electron
UnitContinuum.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_UnitContinuum_h__
8 #define __spatial_UnitContinuum_h__
9 
10 namespace fe
11 {
12 namespace ext
13 {
14 
15 
16 template<typename T, unsigned int D>
17 class UnitContinuum
18 {
19  public:
20  UnitContinuum(void) { }
21 virtual ~UnitContinuum(void) { }
22 
23  UnitContinuum<T,D> &operator=(const UnitContinuum<T,D> &other)
24  {
25  m_count = other.m_count;
26  for(unsigned int i = 0; i < D; i++)
27  {
28  m_periodic[i] = other.m_periodic[i];
29  m_size[i] = other.m_size[i];
30  m_halfsize[i] = other.m_halfsize[i];
31  }
32  m_cells = other.m_cells;
33  m_totalCount = other.m_totalCount;
34  return *this;
35  }
36 
37  bool operator==(const UnitContinuum<T,D> &other)
38  {
39  if(m_count != other.m_count) { return false; }
40  if(m_periodic != other.m_periodic) { return false; }
41  if(m_size != other.m_size) { return false; }
42  if(m_totalCount != other.m_totalCount) { return false; }
43  if(m_cells == other.m_cells) { return true; }
44  return false;
45  }
46 
47  bool operator!=(const UnitContinuum<T,D> &other)
48  {
49  return !operator==(other);
50  }
51 
52  typedef Vector<D,bool> t_periodic;
53  typedef Vector<D,int> t_index;
54  typedef Vector<D,Real> t_span;
55  typedef std::vector<T> t_cells;
56 
57  void create( const t_index &a_count,
58  const t_periodic &a_periodic,
59  const T &a_init);
60 
61  void create( const t_index &a_count,
62  const T &a_init);
63 
64  // safe index oriented access
65  const T &operator[](const t_index &a_index) const;
66  void set(const t_index &a_index, const T &value);
67  void set(const T &value);
68 
69  // unsafe index oriented access
70  T &access(const t_index &a_index);
71  const T &get(const t_index &a_index) const;
72 
73  // unsafe location to index methods (indices may be out of range)
74  void index( t_index &a_index,
75  const t_span &a_location) const;
76  void lower( t_index &a_index,
77  const t_span &a_location) const;
78 
79  // index to location methods
80  void location( t_span &a_location,
81  const t_index &a_index) const;
82 
83  // data structure informational
84  const t_span &cellsize(void) const { return m_size; }
85  const t_index &count(void) const { return m_count; }
86  const t_periodic &periodic(void) const { return m_periodic; }
87  unsigned int totalCount(void) const { return m_totalCount; }
88 
89  // data struction navigational
90  template<typename A>
91  void iterate( void (*a_function)(UnitContinuum<T,D> &,
92  t_index &,
93  A),
94  A a_userarg);
95  void addNeighbors(
96  std::list<t_index> &a_neighbors,
97  const t_index &a_index) const;
98 
99  // out-of-range safety
100  t_index &safe( t_index &a_index) const;
101  t_span &safe( t_span &a_span) const;
102 
103  private:
104  bool safe( unsigned int a_d,
105  t_index &a_index) const;
106  unsigned int calc_array_index( const t_index &a_index) const;
107  template<typename A>
108  void rIterate( unsigned int a_d,
109  t_index &a_index,
110  void (*a_function)(UnitContinuum<T,D> &,
111  t_index &,
112  A),
113  A a_userarg);
114 
115  private:
116  t_index m_count;
117  unsigned int m_totalCount;
118  t_cells m_cells;
119  t_periodic m_periodic;
120  t_span m_size;
121  t_span m_halfsize; // precomputed as optimization
122 };
123 
124 template<typename T, unsigned int D>
125 inline bool UnitContinuum<T,D>::safe(unsigned int a_d, t_index &a_index) const
126 {
127  // return true if the requested index is a valid index
128  bool rv = true;
129  if(m_count[a_d] == 1)
130  {
131  if(a_index[a_d] != 0)
132  {
133  if(!m_periodic[a_d])
134  {
135  rv = false;
136  }
137  a_index[a_d] = 0;
138  }
139  }
140  else if(m_count[a_d] == 2)
141  {
142  if(a_index[a_d] < 0)
143  {
144  if(!m_periodic[a_d])
145  {
146  rv = false;
147  }
148  a_index[a_d] = 0;
149  }
150  else if(a_index[a_d] > 1)
151  {
152  if(!m_periodic[a_d])
153  {
154  rv = false;
155  }
156  a_index[a_d] = 1;
157  }
158  }
159  else
160  {
161  if(a_index[a_d] < 0)
162  {
163  if(m_periodic[a_d])
164  {
165  a_index[a_d] = m_count[a_d] - (-a_index[a_d])%m_count[a_d];
166  }
167  else
168  {
169  if(!m_periodic[a_d])
170  {
171  rv = false;
172  }
173  a_index[a_d] = 0;
174  }
175  }
176  else if(a_index[a_d] >= m_count[a_d])
177  {
178  if(m_periodic[a_d])
179  {
180  a_index[a_d] = a_index[a_d]%m_count[a_d];
181  }
182  else
183  {
184  if(!m_periodic[a_d])
185  {
186  rv = false;
187  }
188  a_index[a_d] = m_count[a_d]-1;
189  }
190  }
191  }
192 
193  return rv;
194 }
195 
196 template<typename T, unsigned int D>
197 inline typename UnitContinuum<T,D>::t_index &UnitContinuum<T,D>::safe(t_index &a_index) const
198 {
199  for(unsigned int i = 0; i < D; i++)
200  {
201  safe(i, a_index);
202  }
203  return a_index;
204 }
205 
206 template<typename T, unsigned int D>
207 inline typename UnitContinuum<T,D>::t_span &UnitContinuum<T,D>::safe(
208  t_span &a_span) const
209 {
210  for(unsigned int i = 0; i < D; i++)
211  {
212  if(m_periodic[i])
213  {
214  float integer;
215  a_span[i] = modff(a_span[i], &integer);
216  if(integer < 0.0)
217  {
218  a_span[i] = 1.0 + a_span[i];
219  }
220  }
221  else
222  {
223  if(a_span[i] < 0.0)
224  {
225  a_span[i] = 0.0;
226  }
227  else if(a_span[i] > 1.0)
228  {
229  a_span[i] = 1.0;
230  }
231  }
232  }
233 
234  return a_span;
235 }
236 
237 template<typename T, unsigned int D>
238 template<typename A>
239 inline void UnitContinuum<T,D>::iterate(void (*a_function)(UnitContinuum<T,D> &,
240 t_index &, A), A a_userarg)
241 {
242  t_index index;
243  rIterate(D-1, index, a_function, a_userarg);
244 }
245 
246 template<typename T, unsigned int D>
247 template<typename A>
248 inline void UnitContinuum<T,D>::rIterate(unsigned int a_d, t_index &a_index,
249  void (*a_function)(UnitContinuum<T,D> &, t_index &, A), A a_userarg)
250 {
251  for(int i = 0; i < count()[a_d]; i++)
252  {
253  a_index[a_d] = i;
254  if(a_d == 0)
255  {
256  a_function(*this, a_index, a_userarg);
257  }
258  else
259  {
260  rIterate(a_d-1, a_index, a_function, a_userarg);
261  }
262  }
263 }
264 
265 template<typename T, unsigned int D>
266 inline void UnitContinuum<T,D>::index(t_index &a_index, const t_span &a_location) const
267 {
268  for(unsigned int i = 0; i < D; i++)
269  {
270  a_index[i] = (int)(a_location[i] / m_size[i]);
271  }
272 }
273 
274 template<typename T, unsigned int D>
275 inline void UnitContinuum<T,D>::lower(t_index &a_index, const t_span &a_location) const
276 {
277  for(unsigned int i = 0; i < D; i++)
278  {
279  a_index[i] = (int)((a_location[i] - m_halfsize[i]) / m_size[i]);
280  }
281 }
282 
283 template<typename T, unsigned int D>
284 inline unsigned int UnitContinuum<T,D>::calc_array_index(
285  const t_index &a_index) const
286 {
287  unsigned int index = 0;
288  unsigned int step = 1;
289  for(int d = D-1; d >= 0; d--)
290  {
291  index += a_index[d] * step;
292  step *= m_count[d];
293  }
294  return index;
295 }
296 
297 template<typename T, unsigned int D>
298 inline void UnitContinuum<T,D>::location(t_span &a_location, const t_index &a_index) const
299 {
300  for(unsigned int i = 0; i < D; i++)
301  {
302  a_location[i] = (Real)a_index[i] * cellsize()[i] + m_halfsize[i];
303  }
304 }
305 
306 template<typename T, unsigned int D>
307 inline void UnitContinuum<T,D>::create(const t_index &a_count,
308  const T &a_init)
309 {
310  m_totalCount = 1;
311  for(unsigned int i = 0; i < D; i++)
312  {
313  m_periodic[i] = false;
314  m_count[i] = a_count[i];
315  m_size[i] = 1.0/m_count[i];
316  m_totalCount *= m_count[i];
317  m_halfsize[i] = m_size[i]/2.0;
318  }
319 
320  m_cells.resize(m_totalCount, a_init);
321 }
322 
323 template<typename T, unsigned int D>
324 inline void UnitContinuum<T,D>::create(const t_index &a_count,
325  const t_periodic &a_periodic, const T &a_init)
326 {
327  m_totalCount = 1;
328  for(unsigned int i = 0; i < D; i++)
329  {
330  m_periodic[i] = a_periodic[i];
331  m_count[i] = a_count[i];
332  m_size[i] = 1.0/m_count[i];
333  m_totalCount *= m_count[i];
334  m_halfsize[i] = m_size[i]/2.0;
335  }
336 
337  m_cells.resize(m_totalCount, a_init);
338 }
339 
340 template<typename T, unsigned int D>
341 inline const T &UnitContinuum<T,D>::operator[](const t_index &a_index) const
342 {
343  t_index safe_index(a_index);
344  safe(safe_index);
345  return m_cells[calc_array_index(safe_index)];
346 }
347 
348 template<typename T, unsigned int D>
349 inline T &UnitContinuum<T,D>::access(const t_index &a_index)
350 {
351  return m_cells[calc_array_index(a_index)];
352 }
353 
354 template<typename T, unsigned int D>
355 inline const T &UnitContinuum<T,D>::get(const t_index &a_index) const
356 {
357  return m_cells[calc_array_index(a_index)];
358 }
359 
360 template<typename T, unsigned int D>
361 inline void UnitContinuum<T,D>::set(const t_index &a_index, const T &a_value)
362 {
363  t_index safe_index(a_index);
364  safe(safe_index);
365  m_cells[calc_array_index(safe_index)] = a_value;
366 }
367 
368 template<typename T, unsigned int D>
369 inline void UnitContinuum<T,D>::set(const T &a_value)
370 {
371  for(unsigned int i = 0; i < m_cells.size(); i++)
372  {
373  m_cells[i] = a_value;
374  }
375 }
376 
377 template<typename T, unsigned int D>
378 inline void UnitContinuum<T,D>::addNeighbors(std::list<t_index> &a_neighbors,
379  const t_index &a_index) const
380 {
381  t_index candidate = a_index;
382  safe(candidate);
383  if(!(candidate == a_index))
384  {
385  feX("UnitContinuum<T,D>::addNeighbors",
386  "invalid reference index");
387  }
388  for(unsigned int i = 0; i < D; i++)
389  {
390  candidate[i] = a_index[i] + 1;
391  safe(i, candidate);
392  if(candidate[i] != a_index[i])
393  {
394  a_neighbors.push_back(candidate);
395  }
396 
397  candidate[i] = a_index[i] - 1;
398  safe(i, candidate);
399  if(candidate[i] != a_index[i])
400  {
401  a_neighbors.push_back(candidate);
402  }
403 
404  candidate[i] = a_index[i];
405  }
406 }
407 
408 // ============================================================================
409 // ============================================================================
410 
411 
412 
413 /** sample a UnitContinuum by returning the direct cell value
414  @relates Continuum */
415 template<typename T, unsigned int D>
416 const T &directSample(const UnitContinuum<T,D> &a_continuum, const typename UnitContinuum<T,D>::t_span &a_location)
417 {
418  typename UnitContinuum<T,D>::t_index idx;
419  a_continuum.index(idx, a_location);
420  return a_continuum[idx];
421 }
422 
423 struct t_bnd
424 {
425  Real m_lo;
426  Real m_hi;
427 };
428 
429 template<typename T, unsigned int D>
430 const T &rLinearSample(unsigned int a_d, const UnitContinuum<T,3> &a_continuum, typename UnitContinuum<T,D>::t_index &a_index, typename UnitContinuum<T,D>::t_index &a_lower, float a_mult, t_bnd *a_bnd)
431 {
432  Real accum = 0.0;
433  Real m;
434  a_index[a_d] = a_lower[a_d];
435  a_continuum.safe(a_index);
436  m = a_bnd[a_d].m_hi * a_mult;
437  if(m > 0.0)
438  {
439  if(a_d == 0)
440  {
441  accum += m * a_continuum.access(a_index);
442  }
443  else
444  {
445  accum += rLinearSample<T,D>(a_d-1, a_continuum, a_index, a_lower,
446  m, a_bnd);
447  }
448  }
449 
450  a_index[a_d] = a_lower[a_d] + 1;
451  a_continuum.safe(a_index);
452  m = a_bnd[a_d].m_lo * a_mult;
453  if(m > 0.0)
454  {
455  if(a_d == 0)
456  {
457  accum += m * a_continuum.access(a_index);
458  }
459  else
460  {
461  accum += rLinearSample<T,D>(a_d-1, a_continuum, a_index, a_lower,
462  m, a_bnd);
463  }
464  }
465 
466  return accum;
467 }
468 
469 /** sample a Continuum by returning the a linear interpolated value.
470  @relates Continuum */
471 template<typename T, unsigned int D>
472 const T &linearSample(const UnitContinuum<T,D> &a_continuum, const typename UnitContinuum<T,D>::t_span &a_location)
473 {
474  t_bnd bnd[D];
475 
476  typename UnitContinuum<T,D>::t_index lower;
477  a_continuum.lower(lower, a_location);
478 
479 
480  for(unsigned int i = 0; i < D; i++)
481  {
482  bnd[i].m_lo = a_location[i]/a_continuum.cellsize()[i]
483  - (Real)lower[i] - 0.5;
484  bnd[i].m_hi = 1.0 - bnd[i].m_lo;
485  }
486 
487  Real r = 1.0;
488 
489  typename UnitContinuum<T,D>::t_index index;
490  return rLinearSample<T,D>(D-1, a_continuum, index, lower, r, bnd);
491 }
492 
493 
494 /** Calculate the gradient of a UnitContinuum at a point using midpoint
495  differencing and linear interpolation.
496  @relates Continuum */
497 template<unsigned int D>
498 void gradient(UnitContinuum<Real,D> &a_continuum, typename UnitContinuum<Real,D>::t_span &a_gradient, const typename UnitContinuum<Real,D>::t_span &a_location, Real a_h = 0.1)
499 {
500  for(unsigned int i = 0; i < D; i++)
501  {
502  typename UnitContinuum<Real,D>::t_span lo(a_location);
503  typename UnitContinuum<Real,D>::t_span hi(a_location);
504  lo[i] -= a_h;
505  hi[i] += a_h;
506  a_continuum.safe(lo);
507  a_continuum.safe(hi);
508  Real lo_smpl = linearSample(a_continuum, lo);
509  Real hi_smpl = linearSample(a_continuum, hi);
510  a_gradient[i] = (hi_smpl - lo_smpl)/(2.0*a_h);
511  }
512 }
513 
514 template<typename T, unsigned int D>
515 struct t_unit_derive
516 {
517  UnitContinuum<T,D> *m_vectors;
518 };
519 
520 template<typename T, unsigned int D>
521 void derive_fun(UnitContinuum<Real,D> &a_continuum, typename UnitContinuum<Real,D>::t_index &a_index, t_unit_derive<T,D> *a_data)
522 {
523  for(unsigned int d = 0; d < D; d++)
524  {
525  typename UnitContinuum<Real,D>::t_index dn(a_index);
526  dn[d] -= 1;
527  typename UnitContinuum<Real,D>::t_index up(a_index);
528  up[d] += 1;
529  Real value = (a_continuum[up] - a_continuum[dn]) /
530  (2.0*a_continuum.cellsize()[d]);
531  setAt(a_data->m_vectors->access(a_index), d, value);
532  }
533 }
534 
535 /** Calculate a gradient vector field for a scalar field.
536  @relates Continuum */
537 template<typename T, unsigned int D>
538 void derive(UnitContinuum<T,D> &a_vector, UnitContinuum<Real,D> &a_scalar)
539 {
540  T dummy;
541  a_vector.create(a_scalar.count(), a_scalar.periodic(), dummy);
542  t_unit_derive<T,D> derive_data;
543  derive_data.m_vectors = &a_vector;
544  a_scalar.iterate(derive_fun<T,D>, &derive_data);
545 }
546 
547 
548 } /* namespace ext */
549 } /* namespace fe */
550 
551 #endif /* __spatial_UnitContinuum_h__ */
fe
kernel
Definition: namespace.dox:3
fe::operator!=
BWORD operator!=(const DualString &s1, const DualString &s2)
Compare two DualString&#39;s (reverse logic)
Definition: DualString.h:229
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