SCIRun  5.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VFDataT.h
Go to the documentation of this file.
1 /*
2  For more information, please see: http://software.sci.utah.edu
3 
4  The MIT License
5 
6  Copyright (c) 2009 Scientific Computing and Imaging Institute,
7  University of Utah.
8 
9 
10  Permission is hereby granted, free of charge, to any person obtaining a
11  copy of this software and associated documentation files (the "Software"),
12  to deal in the Software without restriction, including without limitation
13  the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  and/or sell copies of the Software, and to permit persons to whom the
15  Software is furnished to do so, subject to the following conditions:
16 
17  The above copyright notice and this permission notice shall be included
18  in all copies or substantial portions of the Software.
19 
20  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  DEALINGS IN THE SOFTWARE.
27 */
28 /// @todo Documentation Core/Datatypes/Legacy/Field/VFDataT.h
29 
31 
33 #include <sci_debug.h>
34 
35 #ifdef DEBUG
36 #define TESTRANGE(c, l, h) \
37  if(c < static_cast<index_type>(l) || c >= static_cast<index_type>(h)){ \
38  std::ostringstream msg; \
39  msg << #l "(value=" << l << ") <= " #c << "(value=" << c << ") < " << #h << "(value=" << h << ")"; \
40  throw(SCIRun::AssertionFailed(msg.str().c_str(), __FILE__, __LINE__)); \
41  }
42 #else
43 #define TESTRANGE(c, l, h)
44 #endif
45 
46 
47 #define VFDATAT_ACCESS_DEFINITION(type) \
48 template<class FDATA, class EFDATA, class HFDATA> \
49 void VFDataT<FDATA,EFDATA,HFDATA>::get_value(type &val, VMesh::index_type idx) const \
50 { \
51  TESTRANGE(idx,0,fdata_.size())\
52  val = CastFData<type>(fdata_[idx]); \
53 } \
54 \
55 template<class FDATA, class EFDATA, class HFDATA> \
56 void VFDataT<FDATA,EFDATA,HFDATA>::set_value(const type &val, VMesh::index_type idx) \
57 { \
58  TESTRANGE(idx,0,fdata_.size())\
59  fdata_[idx] = CastFData<typename FDATA::value_type>(val); \
60 } \
61 \
62 template<class FDATA, class EFDATA, class HFDATA> \
63 void VFDataT<FDATA,EFDATA,HFDATA>::get_evalue(type &val, VMesh::index_type idx) const \
64 {\
65  TESTRANGE(idx,0,efdata_.size())\
66  val = CastFData<type>(efdata_[idx]);\
67 } \
68 \
69 template<class FDATA, class EFDATA, class HFDATA> \
70 void VFDataT<FDATA,EFDATA,HFDATA>::set_evalue(const type &val, VMesh::index_type idx) \
71 {\
72  TESTRANGE(idx,0,efdata_.size()) \
73  efdata_[idx] = CastFData<typename FDATA::value_type>(val);\
74 } \
75 \
76 template<class FDATA, class EFDATA, class HFDATA> \
77 void VFDataT<FDATA,EFDATA,HFDATA>::get_values(type *ptr, VMesh::size_type sz, VMesh::size_type offset) const \
78 { if (static_cast<size_type>(fdata_.size()) < sz+offset) sz = static_cast<size_type>(fdata_.size())-offset; for (size_type i=0; i< sz; i++) ptr[i] = CastFData<type>(fdata_[i+offset]); } \
79 \
80 template<class FDATA, class EFDATA, class HFDATA> \
81 void VFDataT<FDATA,EFDATA,HFDATA>::set_values(const type *ptr, VMesh::size_type sz, VMesh::size_type offset) \
82 { if (static_cast<size_type>(fdata_.size()) < sz+offset) sz = static_cast<size_type>(fdata_.size())-offset; for (size_type i=0; i< sz; i++) fdata_[i+offset] = CastFData<typename FDATA::value_type>(ptr[i]); } \
83 \
84 template<class FDATA, class EFDATA, class HFDATA> \
85 void VFDataT<FDATA,EFDATA,HFDATA>::get_evalues(type *ptr, VMesh::size_type sz, VMesh::size_type offset) const \
86 { if (static_cast<size_type>(efdata_.size()) < sz+offset) sz = static_cast<size_type>(efdata_.size())-offset; for (size_type i=0; i< sz; i++) ptr[i] = CastFData<type>(efdata_[i+offset]); } \
87 \
88 template<class FDATA, class EFDATA, class HFDATA> \
89 void VFDataT<FDATA,EFDATA,HFDATA>::set_evalues(const type *ptr, VMesh::size_type sz, VMesh::size_type offset) \
90 { if (static_cast<size_type>(efdata_.size()) < sz+offset) sz = static_cast<size_type>(efdata_.size())-offset; for (size_type i=0; i< sz; i++) efdata_[i+offset] = CastFData<typename FDATA::value_type>(ptr[i]); } \
91 \
92 template<class FDATA, class EFDATA, class HFDATA> \
93 void VFDataT<FDATA,EFDATA,HFDATA>::set_all_values(const type &val) \
94 { typename FDATA::value_type tval = CastFData<typename FDATA::value_type>(val); \
95  size_type sz1 = static_cast<size_type>(fdata_.size()); \
96  size_type sz2 = static_cast<size_type>(efdata_.size()); \
97  for (size_type i=0; i < sz1; i++) fdata_[i] = tval; \
98  for (size_type i=0; i < sz2; i++) efdata_[i] = tval; \
99 } \
100 template<class FDATA, class EFDATA, class HFDATA> \
101 void VFDataT<FDATA,EFDATA,HFDATA>::get_weighted_value(type &val, const VMesh::index_type* idx, const VMesh::weight_type* w, VMesh::size_type sz ) const \
102 { typename FDATA::value_type tval = typename FDATA::value_type(0); for(size_type i=0; i<sz; i++) { TESTRANGE(idx[i],0,fdata_.size()) tval = tval + static_cast<typename FDATA::value_type>(w[i]*fdata_[idx[i]]); } val = CastFData<type>(tval); } \
103 \
104 template<class FDATA, class EFDATA, class HFDATA> \
105 void VFDataT<FDATA,EFDATA,HFDATA>::get_weighted_evalue(type &val, const VMesh::index_type* idx, const VMesh::weight_type* w, VMesh::size_type sz ) const \
106 { typename EFDATA::value_type tval = typename EFDATA::value_type(0); for(size_type i=0; i<sz; i++) { TESTRANGE(idx[i],0,efdata_.size()) tval = tval + static_cast<typename EFDATA::value_type>(w[i]*efdata_[idx[i]]); } val = CastFData<type>(tval); } \
107 \
108 template<class FDATA, class EFDATA, class HFDATA> \
109 void VFDataT<FDATA,EFDATA,HFDATA>::get_values(type *ptr, VMesh::Node::array_type& nodes) const \
110 { for(size_t j=0; j<nodes.size(); j++) { TESTRANGE(nodes[j],0,fdata_.size()) ptr[j] = CastFData<type>(fdata_[nodes[j]]); } } \
111 \
112 template<class FDATA, class EFDATA, class HFDATA> \
113 void VFDataT<FDATA,EFDATA,HFDATA>::get_values(type *ptr, VMesh::Elem::array_type& elems) const\
114 { for(size_t j=0; j<elems.size(); j++) { TESTRANGE(elems[j],0,fdata_.size()) ptr[j] = CastFData<type>(fdata_[elems[j]]); } } \
115 \
116 template<class FDATA, class EFDATA, class HFDATA> \
117 void VFDataT<FDATA,EFDATA,HFDATA>::set_values(const type *ptr, VMesh::Node::array_type& nodes) \
118 { for(size_t j=0; j<nodes.size(); j++) { TESTRANGE(nodes[j],0,fdata_.size()) fdata_[nodes[j]] = CastFData<typename FDATA::value_type>(ptr[j]); } } \
119 \
120 template<class FDATA, class EFDATA, class HFDATA> \
121 void VFDataT<FDATA,EFDATA,HFDATA>::set_values(const type *ptr, VMesh::Elem::array_type& elems) \
122 { for(size_t j=0; j<elems.size(); j++) { TESTRANGE(elems[j],0,fdata_.size()) fdata_[elems[j]] = CastFData<typename FDATA::value_type>(ptr[j]);} } \
123 \
124 template<class FDATA, class EFDATA, class HFDATA> \
125 void VFDataT<FDATA,EFDATA,HFDATA>::get_values(type *ptr, index_type* idx, size_type size) const\
126 { for(index_type j=0; j<size; j++) { TESTRANGE(idx[j],0,fdata_.size()) ptr[j] = CastFData<type>(fdata_[idx[j]]); } } \
127 \
128 template<class FDATA, class EFDATA, class HFDATA> \
129 void VFDataT<FDATA,EFDATA,HFDATA>::set_values(const type *ptr, index_type* idx, size_type size) \
130 { for(index_type j=0; j<size; j++) { TESTRANGE(idx[j],0,fdata_.size()) fdata_[idx[j]] = CastFData<typename FDATA::value_type>(ptr[j]);} } \
131 
132 
133 #define VFDATAT_ACCESS_DEFINITION2(type) \
134 template<class FDATA, class EFDATA, class HFDATA> \
135 void VFDataT<FDATA,EFDATA,HFDATA>::interpolate(type &val, VMesh::ElemInterpolate &interp, type defval) const \
136 { interpolateT<type>(val,interp,defval); } \
137 \
138 template<class FDATA, class EFDATA, class HFDATA> \
139 void VFDataT<FDATA,EFDATA,HFDATA>::minterpolate(std::vector<type> &vals, VMesh::MultiElemInterpolate &interp, type defval) const \
140 { minterpolateT<type>(vals,interp,defval); } \
141 \
142 template<class FDATA, class EFDATA, class HFDATA> \
143 void VFDataT<FDATA,EFDATA,HFDATA>::gradient(StackVector<type,3> &val, VMesh::ElemGradient &interp, type defval) const \
144 { gradientT<type>(val,interp,defval); } \
145 \
146 template<class FDATA, class EFDATA, class HFDATA> \
147 void VFDataT<FDATA,EFDATA,HFDATA>::mgradient(std::vector<StackVector<type,3> > &vals, VMesh::MultiElemGradient &interp, type defval) const \
148 { mgradientT<type>(vals,interp,defval); }
149 
150 #define VFDATA_FUNCTION_SCALAR_DEFINITION(type) \
151 VFData* CreateVFData(std::vector<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
152 { return new VFDataScalarT<std::vector<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
153 \
154 VFData* CreateVFData(Array2<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
155 { return new VFDataScalarT<Array2<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
156 \
157 VFData* CreateVFData(Array3<type>& fdata,std::vector<type>& efdata, std::vector<std::vector<type> >& hfdata) \
158 { return new VFDataScalarT<Array3<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); }
159 
160 #define VFDATA_FUNCTION_VECTOR_DEFINITION(type) \
161 VFData* CreateVFData(std::vector<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
162 { return new VFDataVectorT<std::vector<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
163 \
164 VFData* CreateVFData(Array2<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
165 { return new VFDataVectorT<Array2<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
166 \
167 VFData* CreateVFData(Array3<type>& fdata,std::vector<type>& efdata, std::vector<std::vector<type> >& hfdata) \
168 { return new VFDataVectorT<Array3<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); }
169 
170 #define VFDATA_FUNCTION_TENSOR_DEFINITION(type) \
171 VFData* CreateVFData(std::vector<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
172 { return new VFDataTensorT<std::vector<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
173 \
174 VFData* CreateVFData(Array2<type>& fdata,std::vector<type>& efdata,std::vector<std::vector<type> >& hfdata) \
175 { return new VFDataTensorT<Array2<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); } \
176 \
177 VFData* CreateVFData(Array3<type>& fdata,std::vector<type>& efdata, std::vector<std::vector<type> >& hfdata) \
178 { return new VFDataTensorT<Array3<type>,std::vector<type>,std::vector<std::vector<type> > >(fdata,efdata,hfdata); }
179 
180 
181 namespace SCIRun {
182 
183 /// Implementation layer of the functions we actually need. These classes are
184 /// defined here to prevent overload of the header files:
185 
186 template<class FDATA, class EFDATA, class HFDATA>
187 class VFDataT : public VFData {
188 
189 private:
190  template<class T>
191  void resize(std::vector<T>& fdata,VMesh::dimension_type dim)
192  {
193  VMesh::size_type sz = 1;
194  if (dim.size() > 0) sz = dim[0];
195  if (dim.size() > 1) sz *= dim[1];
196  if (dim.size() > 2) sz *= dim[2];
197  fdata.resize(sz);
198  }
199 
200  template<class T>
201  void resize(Array2<T>& fdata,VMesh::dimension_type dim)
202  {
203  VMesh::size_type sz1 = 1;
204  VMesh::size_type sz2 = 1;
205  if (dim.size() > 0) sz1 = dim[0];
206  if (dim.size() > 1) sz2 = dim[1];
207  fdata.resize(sz2,sz1);
208  }
209 
210  template<class T>
211  void resize(Array3<T>& fdata,VMesh::dimension_type dim)
212  {
213  VMesh::size_type sz1 = 1;
214  VMesh::size_type sz2 = 1;
215  VMesh::size_type sz3 = 1;
216  if (dim.size() > 0) sz1 = dim[0];
217  if (dim.size() > 1) sz2 = dim[1];
218  if (dim.size() > 2) sz3 = dim[2];
219  fdata.resize(sz3,sz2,sz1);
220  }
221 
222 public:
223  // constructor
224  VFDataT(FDATA& fdata, EFDATA& efdata, HFDATA& hfdata) :
225  fdata_(fdata), efdata_(efdata), hfdata_(hfdata)
226  { }
227 
228  // destructor
229  virtual ~VFDataT() {}
230 
231 
233  {
234  resize(fdata_,dim);
235  }
236 
238  {
239  resize(efdata_,dim);
240  }
241 
242  virtual VMesh::size_type fdata_size() const
243  { return (fdata_.size()); }
245  { return (efdata_.size()); }
246 
247  virtual void* fdata_pointer() const
248  {
249  if (fdata_.size() == 0) return (0);
250  return (&(fdata_[0]));
251  }
252  virtual void* efdata_pointer() const
253  {
254  if (efdata_.size() == 0) return (0);
255  return (&(efdata_[0]));
256  }
257 
259  VFDATA_ACCESS_DECLARATION(unsigned char)
261  VFDATA_ACCESS_DECLARATION(unsigned short)
263  VFDATA_ACCESS_DECLARATION(unsigned int)
264  VFDATA_ACCESS_DECLARATION(long long)
265  VFDATA_ACCESS_DECLARATION(unsigned long long)
268  VFDATA_ACCESS_DECLARATION(Core::Geometry::Vector)
269  VFDATA_ACCESS_DECLARATION(Core::Geometry::Tensor)
270 
271 
275  VFDATA_ACCESS_DECLARATION2(Core::Geometry::Vector)
276  VFDATA_ACCESS_DECLARATION2(Core::Geometry::Tensor)
277 
278  virtual void copy_value(VFData* fdata,
279  VMesh::index_type vidx,
280  VMesh::index_type idx)
281  {
282  fdata->get_value(fdata_[idx],vidx);
283  }
284 
285  virtual void copy_values(VFData* fdata,
286  VMesh::index_type vidx,
287  VMesh::index_type idx,
288  VMesh::size_type num)
289  {
290  fdata->get_values(&(fdata_[idx]),num,vidx);
291  }
292 
293 
294  virtual void copy_weighted_value(VFData* fdata,
295  const VMesh::index_type* vidx,
296  const VMesh::weight_type* vw,
297  VMesh::size_type sz,
298  VMesh::index_type idx)
299  {
300  fdata->get_weighted_value(fdata_[idx],vidx,vw,sz);
301  }
302 
303  virtual void copy_evalue(VFData* fdata,
304  VMesh::index_type vidx,
305  VMesh::index_type idx)
306  {
307  typename EFDATA::value_type val;
308  fdata->get_value(val,vidx);
309  efdata_[idx] = val;
310  }
311 
312  virtual void copy_evalues(VFData* fdata,
313  VMesh::index_type vidx,
314  VMesh::index_type idx,
315  VMesh::size_type num)
316  {
317  fdata->get_evalues(&(fdata_[idx]),num,vidx);
318  }
319 
320 
321  virtual void copy_weighted_evalue(VFData* fdata,
322  VMesh::index_type* vidx,
323  VMesh::weight_type* vw,
324  VMesh::size_type sz,
325  VMesh::index_type idx)
326  {
327  typename FDATA::value_type val;
328  fdata->get_weighted_evalue(val,vidx,vw,sz);
329  efdata_[idx] = val;
330  }
331 
332  virtual void copy_values(VFData* fdata)
333  {
334  VMesh::size_type sz1 = fdata_.size();
335  VMesh::size_type sz2 = fdata->fdata_size();
336  if (sz1 == 0) return;
337  if (sz2 < sz1) sz1 = sz2;
338  fdata->get_values(&(fdata_[0]),sz1,0);
339  }
340 
341  virtual void copy_evalues(VFData* fdata)
342  {
343  VMesh::size_type sz1 = efdata_.size();
344  VMesh::size_type sz2 = fdata->efdata_size();
345  if (sz2 < sz1) sz1 = sz2;
346  fdata->get_evalues(&(efdata_[0]),sz1,0);
347  }
348 
349  template<class T>
350  inline void interpolateT(T& val, VMesh::ElemInterpolate& ei, T defval) const
351  {
352  switch(ei.basis_order)
353  {
354  case 0:
355  {
356  if (ei.elem_index >= 0)
357  {
358  TESTRANGE(ei.elem_index,0,fdata_.size())
359  val = CastFData<T>(fdata_[ei.elem_index]);
360  }
361  else
362  val = defval;
363  return;
364  }
365  case 1:
366  {
367  if (ei.elem_index >= 0)
368  {
369  val = static_cast<T>(0.0);
370  for (size_t p=0;p<ei.node_index.size(); p++)
371  {
372  TESTRANGE(ei.node_index[p],0,fdata_.size())
373  val += CastFData<T>(fdata_[ei.node_index[p]]*ei.weights[p]);
374  }
375  }
376  else
377  {
378  val = defval;
379  }
380  return;
381  }
382  case 2:
383  {
384  if (ei.elem_index >= 0)
385  {
386  index_type k = 0;
387  val = static_cast<T>(0.0);
388  for (size_t p=0;p<ei.node_index.size(); p++)
389  {
390  TESTRANGE(ei.node_index[p],0,fdata_.size())
391  val += CastFData<T>(fdata_[ei.node_index[p]]*ei.weights[k]); k++;
392  }
393  for (size_t p=0;p<ei.edge_index.size(); p++)
394  {
395  TESTRANGE(ei.edge_index[p],0,efdata_.size())
396  val += CastFData<T>(efdata_[ei.edge_index[p]]*ei.weights[k]); k++;
397  }
398  }
399  else
400  {
401  val = defval;
402  }
403  return;
404  }
405  case 3:
406  {
407  if (ei.elem_index >= 0)
408  {
409  val = static_cast<T>(0.0);
410  index_type k = 0;
411  for (size_t p=0;p<ei.node_index.size();p++)
412  {
413  TESTRANGE(ei.node_index[p],0,fdata_.size())
414  val += CastFData<T>(fdata_[ei.node_index[p]]*ei.weights[k]); k++;
415  for (index_type q=0;q<ei.num_hderivs;q++)
416  {
417  TESTRANGE(ei.node_index[p],0,hfdata_.size())
418  val += CastFData<T>(hfdata_[ei.node_index[p]][q]*ei.weights[k]); k++;
419  }
420  }
421  }
422  else
423  {
424  val = defval;
425  }
426  return;
427  }
428  }
429  ASSERTFAIL("Interpolate encountered an unknown basis_order");
430  }
431 
432  template<class T>
433  inline void minterpolateT(std::vector<T>& vals, VMesh::MultiElemInterpolate& ei, T defval) const
434  {
435  vals.resize(ei.size());
436 
437  if (ei.size() == 0) return;
438 
439  switch(ei[0].basis_order)
440  {
441  case 0:
442  {
443  for (size_t j=0; j<ei.size(); j++)
444  if (ei[j].elem_index >= 0)
445  vals[j] = CastFData<T>(fdata_[ei[j].elem_index]);
446  else
447  vals[j] = defval;
448  return;
449  }
450  case 1:
451  {
452  for (size_t j=0; j<ei.size(); j++)
453  {
454  if (ei[j].elem_index >= 0)
455  {
456  vals[j] = static_cast<T>(0.0);
457  for (size_t p=0;p<ei[j].node_index.size(); p++)
458  vals[j] += CastFData<T>(fdata_[ei[j].node_index[p]]*ei[j].weights[p]);
459  }
460  else
461  {
462  vals[j] = defval;
463  }
464  }
465  return;
466  }
467  case 2:
468  {
469  for (size_t j=0; j<ei.size(); j++)
470  {
471  if (ei[j].elem_index >= 0)
472  {
473  vals[j] = static_cast<T>(0.0);
474  index_type k = 0;
475  for (size_t p=0;p<ei[j].node_index.size(); p++)
476  { vals[j] += CastFData<T>(fdata_[ei[j].node_index[p]]*ei[j].weights[k]); k++; }
477  for (size_t p=0;p<ei[j].edge_index.size(); p++)
478  { vals[j] += CastFData<T>(efdata_[ei[j].edge_index[p]]*ei[j].weights[k]); k++; }
479  }
480  else
481  {
482  vals[j] = defval;
483  }
484  }
485  return;
486  }
487  case 3:
488  {
489  for (size_t j=0; j<ei.size(); j++)
490  {
491  if (ei[j].elem_index >= 0)
492  {
493  index_type k = 0;
494  vals[j] = static_cast<T>(0.0);
495  for (size_t p=0;p<ei[j].node_index.size();p++)
496  {
497  vals[j] += CastFData<T>(fdata_[ei[j].node_index[p]]*ei[j].weights[k]); k++;
498  for (index_type q=0;q<ei[j].num_hderivs;q++)
499  { vals[j] += CastFData<T>(hfdata_[ei[j].node_index[p]][q]*ei[j].weights[k]); k++; }
500  }
501  }
502  else
503  {
504  vals[j] = defval;
505  }
506  }
507  return;
508  }
509  }
510  ASSERTFAIL("Interpolate encountered an unknown basis_order");
511  }
512 
513 
514  template<class T>
515  inline void mgradientT(std::vector<StackVector<T,3> >& vals, VMesh::MultiElemGradient& eg, T defval) const
516  {
517  vals.resize(eg.size());
518 
519  if (eg.size() == 0) return;
520  T grad;
521 
522  switch(eg[0].basis_order)
523  {
524  case 0:
525  for (size_t j=0; j<eg.size(); j++)
526  {
527  vals[j].resize(3);
528  if (eg[j].elem_index >= 0)
529  {
530  vals[j][0] = T(0); vals[j][1] = T(0); vals[j][2] = T(0);
531  }
532  else
533  {
534  vals[j][0] = defval; vals[j][1] = defval; vals[j][2] = defval;
535  }
536 
537  }
538  return;
539  case 1:
540  {
541  for (size_t j=0; j<eg.size(); j++)
542  {
543  vals[j].resize(3);
544  vals[j][0] = T(0); vals[j][1] = T(0); vals[j][2] = T(0);
545  if (eg[j].elem_index >= 0)
546  {
547  for (index_type k=0, q = 0; k<eg[j].num_derivs; k++)
548  {
549  grad = T(0);
550  for (size_t p=0;p<eg[j].node_index.size();p++)
551  { grad += CastFData<T>(fdata_[eg[j].node_index[p]]*eg[j].weights[q]); q++; }
552 
553  vals[j][0] += CastFData<T>(grad*eg[j].inverse_jacobian[k]);
554  vals[j][1] += CastFData<T>(grad*eg[j].inverse_jacobian[k+3]);
555  vals[j][2] += CastFData<T>(grad*eg[j].inverse_jacobian[k+6]);
556  }
557  }
558  else
559  {
560  vals[j][0] = defval; vals[j][1] = defval; vals[j][2] = defval;
561  }
562  }
563  return;
564  }
565  case 2:
566  {
567  for (size_t j=0; j<eg.size(); j++)
568  {
569  vals[j].resize(3);
570  vals[j][0] = T(0); vals[j][1] = T(0); vals[j][2] = T(0);
571 
572  if (eg[j].elem_index >= 0)
573  {
574  index_type q = 0;
575  for (index_type k=0; k<eg[j].num_derivs; k++)
576  {
577  grad = T(0);
578  for (size_t p=0;p<eg[j].node_index.size();p++)
579  { grad += CastFData<T>(fdata_[eg[j].node_index[p]]*eg[j].weights[q]); q++; }
580  for (size_t p=0;p<eg[j].edge_index.size();p++)
581  { grad += CastFData<T>(fdata_[eg[j].edge_index[p]]*eg[j].weights[q]); q++; }
582 
583  vals[j][0] += CastFData<T>(grad*eg[j].inverse_jacobian[k]);
584  vals[j][1] += CastFData<T>(grad*eg[j].inverse_jacobian[k+3]);
585  vals[j][2] += CastFData<T>(grad*eg[j].inverse_jacobian[k+6]);
586  }
587  }
588  else
589  {
590  vals[j][0] = defval; vals[j][1] = defval; vals[j][2] = defval;
591  }
592 
593  }
594  return;
595  }
596  case 3:
597  {
598  for (size_t j=0; j<eg.size(); j++)
599  {
600  vals[j].resize(3);
601  vals[j][0] = T(0); vals[j][1] = T(0); vals[j][2] = T(0);
602 
603  if (eg[j].elem_index >= 0)
604  {
605  index_type q = 0;
606  for (index_type k=0; k<eg[j].num_derivs; k++)
607  {
608  grad = T(0);
609 
610  for (size_t p=0;p<eg[j].node_index.size();p++)
611  {
612  grad += CastFData<T>(fdata_[eg[j].node_index[p]]*eg[j].weights[q]); q++;
613  for (index_type r=1;r<eg[j].num_hderivs;r++)
614  { grad += CastFData<T>(hfdata_[eg[j].node_index[p]][r]*eg[j].weights[q]); q++; }
615  }
616 
617  vals[j][0] += CastFData<T>(grad*eg[j].inverse_jacobian[k]);
618  vals[j][1] += CastFData<T>(grad*eg[j].inverse_jacobian[k+3]);
619  vals[j][2] += CastFData<T>(grad*eg[j].inverse_jacobian[k+6]);
620  }
621  }
622  else
623  {
624  vals[j][0] = defval; vals[j][1] = defval; vals[j][2] = defval;
625  }
626 
627  }
628  return;
629  }
630  }
631  ASSERTFAIL("Gradient encountered an unknown basis_order");
632  }
633 
634  template<class T>
635  inline void gradientT(StackVector<T,3>& val, VMesh::ElemGradient& eg, T defval) const
636  {
637  T grad;
638 
639  switch(eg.basis_order)
640  {
641  case 0:
642  val.resize(3);
643  if (eg.elem_index >= 0)
644  {
645  val[0] = T(0); val[1] = T(0); val[2] = T(0);
646  }
647  else
648  {
649  val[0] = defval; val[1] = defval; val[2] = defval;
650  }
651 
652  return;
653  case 1:
654  {
655  val.resize(3);
656  val[0] = T(0); val[1] = T(0); val[2] = T(0);
657 
658  if (eg.elem_index >= 0)
659  {
660  for (index_type k=0, q = 0; k<eg.num_derivs; k++)
661  {
662  grad = T(0);
663  for (size_t p=0;p<eg.node_index.size();p++)
664  { grad += CastFData<T>(fdata_[eg.node_index[p]]*eg.weights[q]); q++; }
665 
666  val[0] += CastFData<T>(grad*eg.inverse_jacobian[k]);
667  val[1] += CastFData<T>(grad*eg.inverse_jacobian[k+3]);
668  val[2] += CastFData<T>(grad*eg.inverse_jacobian[k+6]);
669  }
670  }
671  else
672  {
673  val[0] = defval; val[1] = defval; val[2] = defval;
674  }
675  return;
676  }
677  case 2:
678  {
679  val.resize(3);
680  val[0] = T(0); val[1] = T(0); val[2] = T(0);
681 
682  if (eg.elem_index >= 0)
683  {
684  index_type q = 0;
685  for (index_type k=0; k<eg.num_derivs; k++)
686  {
687  grad = T(0);
688  for (size_t p=0;p<eg.node_index.size();p++)
689  { grad += CastFData<T>(fdata_[eg.node_index[p]]*eg.weights[q]); q++; }
690  for (size_t p=0;p<eg.edge_index.size();p++)
691  { grad += CastFData<T>(fdata_[eg.edge_index[p]]*eg.weights[q]); q++; }
692 
693  val[0] += CastFData<T>(grad*eg.inverse_jacobian[k]);
694  val[1] += CastFData<T>(grad*eg.inverse_jacobian[k+3]);
695  val[2] += CastFData<T>(grad*eg.inverse_jacobian[k+6]);
696  }
697  }
698  else
699  {
700  val[0] = defval; val[1] = defval; val[2] = defval;
701  }
702 
703  return;
704  }
705  case 3:
706  {
707  val.resize(3);
708  val[0] = T(0); val[1] = T(0); val[2] = T(0);
709 
710  if (eg.elem_index >= 0)
711  {
712  index_type q = 0;
713  for (index_type k=0; k<eg.num_derivs; k++)
714  {
715  grad = T(0);
716 
717  for (size_t p=0;p<eg.node_index.size();p++)
718  {
719  grad += CastFData<T>(fdata_[eg.node_index[p]]*eg.weights[q]); q++;
720  for (index_type r=1;r<eg.num_hderivs;r++)
721  { grad += CastFData<T>(hfdata_[eg.node_index[p]][r]*eg.weights[q]); q++; }
722  }
723 
724  val[0] += CastFData<T>(grad*eg.inverse_jacobian[k]);
725  val[1] += CastFData<T>(grad*eg.inverse_jacobian[k+3]);
726  val[2] += CastFData<T>(grad*eg.inverse_jacobian[k+6]);
727  }
728  }
729  else
730  {
731  val[0] = defval; val[1] = defval; val[2] = defval;
732  }
733 
734  return;
735  }
736  }
737  ASSERTFAIL("Gradient encountered an unknown basis_order");
738  }
739 
740  virtual VMesh::size_type size() { return (VMesh::size_type(fdata_.size())); }
741 
742 protected:
743  FDATA& fdata_;
744  EFDATA& efdata_; // Additional data for lagrangian interpolation data
745  HFDATA& hfdata_; // Additional data for hermitian interpolation data
746 };
747 
748 
749 
750 
751 
752 template<class FDATA, class EFDATA, class HFDATA>
753 class VFDataScalarT : public VFDataT<FDATA,EFDATA,HFDATA> {
754 public:
755  // constructor
756  VFDataScalarT(FDATA& fdata, EFDATA& efdata, HFDATA& hfdata) :
757  VFDataT<FDATA,EFDATA,HFDATA>(fdata,efdata,hfdata)
758  {}
759 
760  // destructor
761  virtual ~VFDataScalarT() {}
762 
763  virtual bool min(double& val, VMesh::index_type& idx) const
764  {
765  typename FDATA::value_type tval(0);
766  idx = 0;
767  size_type sz = static_cast<size_type>(this->fdata_.size());
768  if (sz > 0)
769  tval = this->fdata_[0];
770  else
771  return (false);
772 
773  for (size_type p=1; p<sz; p++)
774  {
775  if (this->fdata_[p] < tval)
776  {
777  tval = this->fdata_[p];
778  idx = VMesh::index_type(p);
779  }
780  }
781  val = CastFData<double>(tval);
782  return (true);
783  }
784 
785  virtual bool max(double& val, VMesh::index_type& idx) const
786  {
787  typename FDATA::value_type tval(0);
788  idx = 0;
789  size_type sz = static_cast<size_type>(this->fdata_.size());
790  if (sz > 0) tval = this->fdata_[0]; else return (false);
791  for (size_type p=1; p<sz; p++)
792  if (this->fdata_[p] > tval) { tval = this->fdata_[p]; idx = VMesh::index_type(p); }
793  val = CastFData<double>(tval);
794  return (true);
795  }
796 
797  virtual bool minmax(double& min, VMesh::index_type& idxmin,
798  double& max, VMesh::index_type& idxmax) const
799  {
800  typename FDATA::value_type tval(0);
801  typename FDATA::value_type tval2(0);
802  idxmin = 0;
803  idxmax = 0;
804  size_type sz = static_cast<size_type>(this->fdata_.size());
805  if (sz > 0) { tval = this->fdata_[0]; tval2 = tval; } else return (false);
806  for (size_type p=1; p<sz; p++)
807  {
808  if (this->fdata_[p] < tval) { tval = this->fdata_[p]; idxmin = VMesh::index_type(p); }
809  if (this->fdata_[p] > tval2) { tval2 = this->fdata_[p]; idxmax = VMesh::index_type(p); }
810  }
811  min = CastFData<double>(tval);
812  max = CastFData<double>(tval2);
813  return (true);
814  }
815 };
816 
817 template<class FDATA, class EFDATA, class HFDATA>
818 class VFDataVectorT : public VFDataT<FDATA,EFDATA,HFDATA> {
819 public:
820  // constructor
821  VFDataVectorT(FDATA& fdata, EFDATA& efdata, HFDATA& hfdata) :
822  VFDataT<FDATA,EFDATA,HFDATA>(fdata,efdata,hfdata)
823  {}
824 
825  // destructor
826  virtual ~VFDataVectorT() {}
827 
828  virtual bool min(double& val, VMesh::index_type& idx) const
829  {
830  double tval = 0;
831  idx = 0;
832  size_type sz = static_cast<size_type>(this->fdata_.size());
833  if (sz > 0) tval = this->fdata_[0].length();
834 
835  for (size_type p=1; p<sz; p++)
836  {
837  double len = this->fdata_[p].length();
838  if (len < tval) { tval = len; idx = VMesh::index_type(p); }
839  }
840 
841  val = CastFData<double>(tval);
842  return (true);
843  }
844 
845  virtual bool max(double& val, VMesh::index_type& idx) const
846  {
847  double tval = 0;
848  idx = 0;
849  size_type sz = static_cast<size_type>(this->fdata_.size());
850  if (sz > 0) tval = this->fdata_[0].length();
851  for (size_type p=1; p<sz; p++)
852  {
853  double len = this->fdata_[p].length();
854  if (len > tval) { tval = len; idx = VMesh::index_type(p); }
855  }
856 
857  val = CastFData<double>(tval);
858  return (true);
859  }
860 
861  virtual bool minmax(double& min, VMesh::index_type& idxmin,
862  double& max, VMesh::index_type& idxmax) const
863  {
864  double tval = 0;
865  double tval2 = 0;
866  idxmin = 0;
867  idxmax = 0;
868  size_type sz = static_cast<size_type>(this->fdata_.size());
869  if (sz > 0) { tval = this->fdata_[0].length(); tval2 = tval; }
870  for (size_type p=1; p<sz; p++)
871  {
872  if (this->fdata_[p].length() < tval) { tval = this->fdata_[p].length(); idxmin = VMesh::index_type(p); }
873  if (this->fdata_[p].length() > tval2) { tval2 = this->fdata_[p].length(); idxmax = VMesh::index_type(p); }
874  }
875  min = CastFData<double>(tval);
876  max = CastFData<double>(tval2);
877  return (true);
878  }
879 };
880 
881 
882 template<class FDATA, class EFDATA, class HFDATA>
883 class VFDataTensorT : public VFDataT<FDATA,EFDATA,HFDATA> {
884 public:
885  // constructor
886  VFDataTensorT(FDATA& fdata, EFDATA& efdata, HFDATA& hfdata) :
887  VFDataT<FDATA,EFDATA,HFDATA>(fdata,efdata,hfdata)
888  {}
889 
890  // destructor
891  virtual ~VFDataTensorT() {}
892 
893  virtual bool min(double& val, VMesh::index_type& idx) const
894  {
895  double tval = 0;
896  idx = 0;
897  size_type sz = static_cast<size_type>(this->fdata_.size());
898  if (sz > 0) tval = this->fdata_[0].norm();
899 
900  for (size_type p=1; p<sz; p++)
901  {
902  double len = this->fdata_[p].norm();
903  if (len < tval) { tval = len; idx = VMesh::index_type(p); }
904  }
905 
906  val = CastFData<double>(tval);
907  return (true);
908  }
909 
910  virtual bool max(double& val, VMesh::index_type& idx) const
911  {
912  double tval = 0;
913  idx = 0;
914  size_type sz = static_cast<size_type>(this->fdata_.size());
915  if (sz > 0) tval = this->fdata_[0].norm();
916  for (size_type p=1; p<sz; p++)
917  {
918  double len = this->fdata_[p].norm();
919  if (len > tval) { tval = len; idx = VMesh::index_type(p); }
920  }
921 
922  val = CastFData<double>(tval);
923  return (true); }
924 
925  virtual bool minmax(double& min, VMesh::index_type& idxmin,
926  double& max, VMesh::index_type& idxmax) const
927  {
928  double tval = 0;
929  double tval2 = 0;
930  idxmin = 0;
931  idxmax = 0;
932  size_type sz = static_cast<size_type>(this->fdata_.size());
933  if (sz > 0) { tval = this->fdata_[0].norm(); tval2 = tval; }
934  for (size_type p=1; p<sz; p++)
935  {
936  if (this->fdata_[p].norm() < tval) { tval = this->fdata_[p].norm(); idxmin = VMesh::index_type(p); }
937  if (this->fdata_[p].norm() > tval2) { tval2 = this->fdata_[p].norm(); idxmax = VMesh::index_type(p); }
938  }
939  min = CastFData<double>(tval);
940  max = CastFData<double>(tval2);
941  return (true);
942  }
943 };
944 
945 
947 VFDATAT_ACCESS_DEFINITION(Core::Geometry::Vector)
948 VFDATAT_ACCESS_DEFINITION(Core::Geometry::Tensor)
950 VFDATAT_ACCESS_DEFINITION(unsigned char)
952 VFDATAT_ACCESS_DEFINITION(unsigned short)
954 VFDATAT_ACCESS_DEFINITION(unsigned int)
956 VFDATAT_ACCESS_DEFINITION(long long)
957 VFDATAT_ACCESS_DEFINITION(unsigned long long)
958 
962 VFDATAT_ACCESS_DEFINITION2(Core::Geometry::Vector)
963 VFDATAT_ACCESS_DEFINITION2(Core::Geometry::Tensor)
964 
965 
966 }
967 
968 
969 
Definition: VFDataT.h:753
std::vector< ElemInterpolate > MultiElemInterpolate
Definition: VMesh.h:208
virtual void copy_values(VFData *fdata, VMesh::index_type vidx, VMesh::index_type idx, VMesh::size_type num)
Definition: VFDataT.h:285
StackVector< index_type, 8 > node_index
Definition: VMesh.h:199
#define VFDATA_ACCESS_DECLARATION(type)
Definition: VFData.h:47
virtual ~VFDataScalarT()
Definition: VFDataT.h:761
virtual VMesh::size_type efdata_size() const
Definition: VFDataT.h:244
void gradientT(StackVector< T, 3 > &val, VMesh::ElemGradient &eg, T defval) const
Definition: VFDataT.h:635
Definition: VMesh.h:212
virtual void resize_efdata(VMesh::dimension_type dim)
Definition: VFDataT.h:237
StackVector< double, 9 > inverse_jacobian
Definition: VMesh.h:224
virtual bool minmax(double &min, VMesh::index_type &idxmin, double &max, VMesh::index_type &idxmax) const
Definition: VFDataT.h:797
Generic exception for internal errors.
virtual bool min(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:763
index_type elem_index
Definition: VMesh.h:198
virtual void copy_values(VFData *fdata)
Copy values from one FData array to another FData array.
Definition: VFDataT.h:332
virtual void copy_evalues(VFData *fdata, VMesh::index_type vidx, VMesh::index_type idx, VMesh::size_type num)
Definition: VFDataT.h:312
StackBasedVector< double, 64 > weights
Definition: VMesh.h:202
virtual void * efdata_pointer() const
Definition: VFDataT.h:252
virtual VMesh::size_type fdata_size() const
Definition: VFData.cc:105
VFDataTensorT(FDATA &fdata, EFDATA &efdata, HFDATA &hfdata)
Definition: VFDataT.h:886
void minterpolateT(std::vector< T > &vals, VMesh::MultiElemInterpolate &ei, T defval) const
Definition: VFDataT.h:433
virtual void copy_weighted_evalue(VFData *fdata, VMesh::index_type *vidx, VMesh::weight_type *vw, VMesh::size_type sz, VMesh::index_type idx)
Copy a weighted edge value without needing to know the type.
Definition: VFDataT.h:321
HFDATA & hfdata_
Definition: VFDataT.h:745
Definition: VMesh.h:193
virtual ~VFDataT()
Definition: VFDataT.h:229
Definition: Vector.h:63
virtual bool max(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:785
virtual bool max(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:845
#define VFDATAT_ACCESS_DEFINITION(type)
Definition: VFDataT.h:47
virtual ~VFDataVectorT()
Definition: VFDataT.h:826
#define VFDATA_ACCESS_DECLARATION2(type)
Definition: VFData.h:67
size_type num_hderivs
Definition: VMesh.h:220
Definition: VFDataT.h:883
Definition: ParallelLinearAlgebraTests.cc:358
void resize(size_t size1, size_t size2, size_t size3)
Definition: Array3.h:77
#define TESTRANGE(c, l, h)
Definition: VFDataT.h:43
virtual VMesh::size_type size()
Definition: VFDataT.h:740
long long size_type
Definition: Types.h:40
Definition: StackVector.h:50
Definition: VFDataT.h:818
StackBasedVector< double, 64 > weights
Definition: VMesh.h:221
virtual void copy_value(VFData *fdata, VMesh::index_type vidx, VMesh::index_type idx)
Copy a value without needing to know the type.
Definition: VFDataT.h:278
Definition: VFData.h:86
virtual void copy_weighted_value(VFData *fdata, const VMesh::index_type *vidx, const VMesh::weight_type *vw, VMesh::size_type sz, VMesh::index_type idx)
Definition: VFDataT.h:294
double weight_type
Weights for interpolation.
Definition: VMesh.h:65
virtual bool min(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:828
virtual VMesh::size_type fdata_size() const
Definition: VFDataT.h:242
virtual VMesh::size_type efdata_size() const
Definition: VFData.cc:111
FDATA & fdata_
Definition: VFDataT.h:743
Definition: Tensor.h:65
virtual bool min(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:893
VFDataVectorT(FDATA &fdata, EFDATA &efdata, HFDATA &hfdata)
Definition: VFDataT.h:821
Definition: Array2.h:51
#define ASSERTFAIL(string)
Definition: Assert.h:52
StackVector< index_type, 12 > edge_index
Definition: VMesh.h:200
virtual void copy_evalue(VFData *fdata, VMesh::index_type vidx, VMesh::index_type idx)
Copy a edge value without needing to know the type.
Definition: VFDataT.h:303
VFDataScalarT(FDATA &fdata, EFDATA &efdata, HFDATA &hfdata)
Definition: VFDataT.h:756
EFDATA & efdata_
Definition: VFDataT.h:744
virtual ~VFDataTensorT()
Definition: VFDataT.h:891
long long index_type
Definition: Types.h:39
Definition: ParallelLinearAlgebraTests.cc:751
virtual void copy_evalues(VFData *fdata)
Definition: VFDataT.h:341
void resize(size_t size1, size_t size2)
Definition: Array2.h:64
void resize(size_t size, const value_type &val=value_type())
Definition: StackVector.h:61
void mgradientT(std::vector< StackVector< T, 3 > > &vals, VMesh::MultiElemGradient &eg, T defval) const
Definition: VFDataT.h:515
virtual bool minmax(double &min, VMesh::index_type &idxmin, double &max, VMesh::index_type &idxmax) const
Definition: VFDataT.h:861
virtual void * fdata_pointer() const
Definition: VFDataT.h:247
virtual void resize_fdata(VMesh::dimension_type dim)
Definition: VFDataT.h:232
size_type num_derivs
Definition: VMesh.h:223
#define VFDATAT_ACCESS_DEFINITION2(type)
Definition: VFDataT.h:133
virtual bool minmax(double &min, VMesh::index_type &idxmin, double &max, VMesh::index_type &idxmax) const
Definition: VFDataT.h:925
void interpolateT(T &val, VMesh::ElemInterpolate &ei, T defval) const
Definition: VFDataT.h:350
StackBasedVector< index_type, 8 > node_index
Definition: VMesh.h:218
VFDataT(FDATA &fdata, EFDATA &efdata, HFDATA &hfdata)
Definition: VFDataT.h:224
size_type num_hderivs
Definition: VMesh.h:201
std::vector< size_type > dimension_type
Definition: VMesh.h:76
Definition: Array3.h:64
Mesh::size_type size_type
Definition: VMesh.h:68
int basis_order
Definition: VMesh.h:216
StackVector< index_type, 12 > edge_index
Definition: VMesh.h:219
Definition: VMesh.h:53
virtual bool max(double &val, VMesh::index_type &idx) const
Definition: VFDataT.h:910
int basis_order
Definition: VMesh.h:197
index_type elem_index
Definition: VMesh.h:217
Mesh::index_type index_type
VIRTUAL INTERFACE.
Definition: VMesh.h:63
std::vector< ElemGradient > MultiElemGradient
Multiple gradient calculations in parallel.
Definition: VMesh.h:229
Definition: VFDataT.h:187