Dendro  5.01
Dendro in Greek language means tree. The Dendro library is a large scale (262K cores on ORNL's Titan) distributed memory adaptive octree framework. The main goal of Dendro is to perform large scale multiphysics simulations efficeiently in mordern supercomputers. Dendro consists of efficient parallel data structures and algorithms to perform variational ( finite element) methods and finite difference mthods on 2:1 balanced arbitary adaptive octrees which enables the users to perform simulations raning from black holes (binary black hole mergers) to blood flow in human body, where applications ranging from relativity, astrophysics to biomedical engineering.
matvec.h
1 //
2 // Created by milinda on 11/20/17.
13 //
14 
15 #ifndef SFCSORTBENCH_FEMMATVEC_H
16 #define SFCSORTBENCH_FEMMATVEC_H
17 
18 #define NUM_VERTEX_BUCKETS 27
19 #define NUM_EDGE_BUCKETS 54
20 #define NUM_FACE_BUCKETS 36
21 #define NUM_INTERNAL_BUCKETS 8
22 #define NUM_TOTAL_BUCKETS 125
23 
24 #define IDX(i,j,k) i+nx*(j+k*ny)
25 #define IDX2D(i,j) j*ny+i
26 
27 #include "binUtils.h"
28 #include "dendro.h"
29 #include "refel.h"
30 #include "assert.h"
31 #include "TreeNode.h"
32 #include "mesh.h"
33 #include "operators.h"
34 
35 #define PT_EXTERNAL 0 // if the point is external to a given element
36 #define PT_INTERNAL 1 // if the point is internal but not a DG node.
37 #define PT_INTERNAL_DG_NODE 2 // if the point is internal and it is a DG node.
38 
39 #define SFC_MVEC_TABLE_DEFAULT UCHAR_MAX
40 #define SFC_MVEC_TABLE_Rz 9
41 #define SFC_MVEC_TABLE_Ry 3
42 
43 
44 namespace fem
45 {
46  namespace seq
47  {
49  //m_uiDim 3
50  static const unsigned char SFC_MATVEC_BUCKET_TABLE[m_uiDim*m_uiDim*m_uiDim][NUM_CHILDREN]={
51  {0,255,255,255,255,255,255,255},
52  {0,1,255,255,255,255,255,255},
53  {1,255,255,255,255,255,255,255},
54  {0,2,255,255,255,255,255,255},
55  {0,1,2,3,255,255,255,255},
56  {1,3,255,255,255,255,255,255},
57  {2,255,255,255,255,255,255,255},
58  {2,3,255,255,255,255,255,255},
59  {3,255,255,255,255,255,255,255},
60  {0,4,255,255,255,255,255,255},
61  {0,1,4,5,255,255,255,255},
62  {1,5,255,255,255,255,255,255},
63  {0,2,4,6,255,255,255,255},
64  {0,1,2,3,4,5,6,7},
65  {1,3,5,7,255,255,255,255},
66  {2,6,255,255,255,255,255,255},
67  {2,3,6,7,255,255,255,255},
68  {3,7,255,255,255,255,255,255},
69  {4,255,255,255,255,255,255,255},
70  {4,5,255,255,255,255,255,255},
71  {5,255,255,255,255,255,255,255},
72  {4,6,255,255,255,255,255,255},
73  {4,5,6,7,255,255,255,255},
74  {5,7,255,255,255,255,255,255},
75  {6,255,255,255,255,255,255,255},
76  {6,7,255,255,255,255,255,255},
77  {7,255,255,255,255,255,255,255}
78  };
79 
80  }
81 }
82 
83 
84 namespace fem
85 {
86 
87  // Note: if possible remove this sequenction namespace once the FEM mesh free matvec implemented for parallel case as well.
88  namespace seq
89  {
90 
91  /***
92  * @brief : computes i,j,k for a given nodal point relative to an given element.
93  *
94  * */
95  template<typename T>
96  inline unsigned int computeIJK(const T& pt,const ot::TreeNode elem,const unsigned int dx,const unsigned int logDx,unsigned int *ijk);
97 
98 
99  template<typename T>
100  inline bool isInternal(const T& pt,const ot::TreeNode elem,const unsigned int dx,const unsigned int logDx);
101 
108  template <typename T,typename da>
109  inline void matvec(const T *pt_coords,const da* in,unsigned int pBegin, unsigned int pEnd,const unsigned int lev,const unsigned int maxDepth,const unsigned int order,const ot::TreeNode& pOctant,const RefElement * refEl,da * out,unsigned int mVecType);
110 
111 
112 
113 
114  }
115 
116 
117 
118 
119 
120 }
121 
122 template<typename T>
123 inline unsigned int fem::seq::computeIJK(const T& pt,const ot::TreeNode elem,const unsigned int dx,const unsigned int logDx,unsigned int *ijk)
124 {
125  /*ijk[0]=UINT_MAX;
126  ijk[1]=UINT_MAX;
127  ijk[2]=UINT_MAX;*/
128 
129  //dendro::timer::sfcmatvec::t_computeIJK.start();
130  /*const unsigned int upperX=elem.maxX();
131  const unsigned int upperY=elem.maxY();
132  const unsigned int upperZ=elem.maxZ();
133 
134  const unsigned int lowerX=elem.minX();
135  const unsigned int lowerY=elem.minY();
136  const unsigned int lowerZ=elem.minZ();
137 
138  const unsigned int x=pt.xint();
139  const unsigned int y=pt.yint();
140  const unsigned int z=pt.zint();*/
141 
142  /*int p = (x | ~upper) & ((x ^ upper) | (~upper + x));
143  int q = (lower | ~x) & ((lower ^ x) | (~x + lower));
144  return 1 & ((p & q) >> 31);*/
145 
146  //if((((pt.xint() | ~elem.maxX()) & ((pt.xint() ^ elem.maxX()) | (~elem.maxX() + pt.xint())))>>31u) && (((elem.minX() | ~pt.xint()) & ((elem.minX() ^ pt.xint()) | (~pt.xint() + elem.minX())))>>31u) && (((pt.yint() | ~elem.maxY()) & ((pt.yint() ^ elem.maxY()) | (~elem.maxY() + pt.yint())))>>31u) && (((elem.minY() | ~pt.yint()) & ((elem.minY() ^ pt.yint()) | (~pt.yint() + elem.minY())))>>31u) && (((pt.zint() | ~elem.maxZ()) & ((pt.zint() ^ elem.maxZ()) | (~elem.maxZ() + pt.zint()))) && ((elem.minZ() | ~pt.zint()) & ((elem.minZ() ^ pt.zint()) | (~pt.zint() + elem.minZ())))>>31u) )
147  if(elem.minX()<=pt.xint() && elem.minY()<=pt.yint() && elem.minZ()<=pt.zint() && pt.xint()<=elem.maxX() && pt.yint()<=elem.maxY() && pt.zint() <=elem.maxZ())
148  {
149  // check if diff %dx ==0
150  if(!((pt.xint()-elem.minX())&((1u<<logDx)-1)) && !((pt.yint()-elem.minY())&((1u<<logDx)-1)) && !((pt.zint()-elem.minZ())&((1u<<(logDx))-1)) )
151  {
152  ijk[0]=(pt.xint()-elem.minX())>>logDx;
153  ijk[1]=(pt.yint()-elem.minY())>>logDx;
154  ijk[2]=(pt.zint()-elem.minZ())>>logDx;
155 
156  //dendro::timer::sfcmatvec::t_computeIJK.stop();
157 
158  return PT_INTERNAL_DG_NODE;
159  }else {
160  //dendro::timer::sfcmatvec::t_computeIJK.stop();
161  return PT_INTERNAL;
162  }
163 
164 
165  }else {
166  //dendro::timer::sfcmatvec::t_computeIJK.stop();
167  return PT_EXTERNAL;
168  }
169 
170 
171 }
172 
173 
174 template<typename T>
175 inline bool fem::seq::isInternal(const T& pt,const ot::TreeNode elem,const unsigned int dx,const unsigned int logDx)
176 {
177 
178  return (elem.minX()<=pt.xint() && elem.minY()<=pt.yint() && elem.minZ()<=pt.zint() && pt.xint()<=elem.maxX() && pt.yint()<=elem.maxY() && pt.zint() <=elem.maxZ()) ? true : false;
179 }
180 
181 
182 
183 template <typename T,typename da>
184 inline void fem::seq::matvec(const T *pt_coords,const da* in,unsigned int pBegin, unsigned int pEnd,const unsigned int lev,const unsigned int maxDepth,const unsigned int order,const ot::TreeNode& pOctant,const RefElement * refEl,da * out,unsigned int mVecType){
185 
186 
187  const unsigned int nPe=(order+1)*(order+1)*(order+1);
188  const unsigned int pMaxDepthBit=maxDepth-lev-1;
189  assert(lev==pOctant.getLevel());
190  assert(((1u<<(maxDepth-lev)) % order)==0);
191 
192  const unsigned int sz=(1u<<(maxDepth-pOctant.getLevel()));
193  const unsigned int szb2=sz>>1u;
194 
195  const unsigned int dx=(1u<<(maxDepth-lev))/order;
196  const unsigned int dxb2=(1u<<(maxDepth-lev-1))/order;
197  const unsigned int logDx=binOp::fastLog2(dx);
198  const unsigned int logDxb2=binOp::fastLog2(dxb2);
199 
200  const unsigned int nx=order+1;
201  const unsigned int ny=order+1;
202  const unsigned int nz=order+1;
203 
204  register unsigned int ijk[3];
205  unsigned int cnum;
206  const unsigned int pMin[]={pOctant.minX(),pOctant.minY(),pOctant.minZ()};
207  const unsigned int pMax[]={pOctant.maxX(),pOctant.maxY(),pOctant.maxZ()};
208  const unsigned int pMid[]={(pMax[0]+pMin[0])>>1u,(pMax[1]+pMin[1])>>1u,(pMax[2]+pMin[2])>>1u};
209 
210  unsigned int bCounts[NUM_CHILDREN];
211  for(unsigned int child=0;child<NUM_CHILDREN;child++)
212  bCounts[child]=0;
213 
214 
215 
216  for(unsigned int p=pBegin;p<pEnd;p++)
217  out[p]=(da)0; // initialize the out vector to zero.
218 
219 
220 
221  ot::TreeNode childElem[NUM_CHILDREN];
222 
223 
224  register unsigned int pt_status;
225  childElem[0]=ot::TreeNode(pOctant.minX(),pOctant.minY(),pOctant.minZ(),pOctant.getLevel()+1,m_uiDim,maxDepth);
226  childElem[1]=ot::TreeNode(pOctant.minX()+szb2,pOctant.minY(),pOctant.minZ(),pOctant.getLevel()+1,m_uiDim,maxDepth);
227  childElem[2]=ot::TreeNode(pOctant.minX(),pOctant.minY()+szb2,pOctant.minZ(),pOctant.getLevel()+1,m_uiDim,maxDepth);
228  childElem[3]=ot::TreeNode(pOctant.minX()+szb2,pOctant.minY()+szb2,pOctant.minZ(),pOctant.getLevel()+1,m_uiDim,maxDepth);
229 
230  childElem[4]=ot::TreeNode(pOctant.minX(),pOctant.minY(),pOctant.minZ()+szb2,pOctant.getLevel()+1,m_uiDim,maxDepth);
231  childElem[5]=ot::TreeNode(pOctant.minX()+szb2,pOctant.minY(),pOctant.minZ()+szb2,pOctant.getLevel()+1,m_uiDim,maxDepth);
232  childElem[6]=ot::TreeNode(pOctant.minX(),pOctant.minY()+szb2,pOctant.minZ()+szb2,pOctant.getLevel()+1,m_uiDim,maxDepth);
233  childElem[7]=ot::TreeNode(pOctant.minX()+szb2,pOctant.minY()+szb2,pOctant.minZ()+szb2,pOctant.getLevel()+1,m_uiDim,maxDepth);
234 
235  // keep track of the interpolations needed.
236 
237 
238 #ifdef MATVEC_PROFILE
239  dendro::timer::sfcmatvec::t_pts_p1_count.start();
240 #endif
241 
242  std::vector<unsigned int> duplicateCoordIndex;
243  //pass 1
244  for(unsigned int p=pBegin;p<pEnd;p++)
245  {
246 
247  ijk[2]=0;
248  ijk[1]=0;
249  ijk[0]=0;
250 
251  if(pt_coords[p].zint()>pMid[2]) ijk[2]=2;
252  if(pt_coords[p].yint()>pMid[1]) ijk[1]=2;
253  if(pt_coords[p].xint()>pMid[0]) ijk[0]=2;
254 
255  if(pt_coords[p].zint()<pMid[2]) ijk[2]=0;
256  if(pt_coords[p].yint()<pMid[1]) ijk[1]=0;
257  if(pt_coords[p].xint()<pMid[0]) ijk[0]=0;
258 
259  if(pt_coords[p].zint()==pMid[2]) ijk[2]=1;
260  if(pt_coords[p].yint()==pMid[1]) ijk[1]=1;
261  if(pt_coords[p].xint()==pMid[0]) ijk[0]=1;
262 
263  if(ijk[0]!=1 && ijk[1]!=1 && ijk[2]!=1)
264  {
265 
266  cnum=(((ijk[2]>>1u)<<2u) | ((ijk[1]>>1u)<<1u) | ((ijk[0]>>1u)));
267  bCounts[cnum]++;
268 
269  }else{
270  duplicateCoordIndex.push_back(p);
271  }
272 
273  }
274 
275 #ifdef MATVEC_PROFILE
276  dendro::timer::sfcmatvec::t_pts_p1_count.stop();
277 #endif
278 
279 #ifdef MATVEC_PROFILE
280  dendro::timer::sfcmatvec::t_pts_p2_count.start();
281 #endif
282 
283  //pass 2
284  for(unsigned int index=0;index<duplicateCoordIndex.size();index++)
285  {
286  ijk[2]=0;
287  ijk[1]=0;
288  ijk[0]=0;
289 
290  if(pt_coords[duplicateCoordIndex[index]].zint()>pMid[2]) ijk[2]=2;
291  if(pt_coords[duplicateCoordIndex[index]].yint()>pMid[1]) ijk[1]=2;
292  if(pt_coords[duplicateCoordIndex[index]].xint()>pMid[0]) ijk[0]=2;
293 
294  if(pt_coords[duplicateCoordIndex[index]].zint()<pMid[2]) ijk[2]=0;
295  if(pt_coords[duplicateCoordIndex[index]].yint()<pMid[1]) ijk[1]=0;
296  if(pt_coords[duplicateCoordIndex[index]].xint()<pMid[0]) ijk[0]=0;
297 
298  if(pt_coords[duplicateCoordIndex[index]].zint()==pMid[2]) ijk[2]=1;
299  if(pt_coords[duplicateCoordIndex[index]].yint()==pMid[1]) ijk[1]=1;
300  if(pt_coords[duplicateCoordIndex[index]].xint()==pMid[0]) ijk[0]=1;
301 
302  for(unsigned int child=0;child<NUM_CHILDREN;child++)
303  {
304 
305  cnum=SFC_MATVEC_BUCKET_TABLE[(ijk[2]*SFC_MVEC_TABLE_Rz+ijk[1]*SFC_MVEC_TABLE_Ry+ijk[0])][child];
306  if(cnum==UCHAR_MAX)
307  break;
308 
309  //std::cout<<"pt : "<<pt_coords[p].xint()<<" "<<pt_coords[p].yint()<<" "<<pt_coords[p].zint()<<" ijk: "<<ijk[0]<<" "<<ijk[1]<<" "<<ijk[2]<<" child: "<<child <<" cnum "<<cnum<<" pMid: "<<pMid[2]<<std::endl;
310  // std::cout<<"pt : "<<pt_coords[p].xint()<<" "<<pt_coords[p].yint()<<" "<<pt_coords[p].zint()<<" ijk: "<<ijk[0]<<" "<<ijk[1]<<" "<<ijk[2]<<" cnum "<<cnum<<" pMid: "<<pMid[2]<<"table 0 : "<<(int)SFC_MATVEC_BUCKET_TABLE[(ijk[2]*SFC_MVEC_TABLE_Rz+ijk[1]*SFC_MVEC_TABLE_Ry+ijk[0])][0]<<" table 1: "<<(int)SFC_MATVEC_BUCKET_TABLE[(ijk[2]*SFC_MVEC_TABLE_Rz+ijk[1]*SFC_MVEC_TABLE_Ry+ijk[0])][1]<<std::endl;
311  bCounts[cnum]++;
312  }
313 
314  }
315 
316 
317 #ifdef MATVEC_PROFILE
318  dendro::timer::sfcmatvec::t_pts_p2_count.stop();
319 #endif
320 
321 
322 
323 
324 #ifdef MATVEC_PROFILE
325  dendro::timer::sfcmatvec::t_malloc.start();
326 #endif
327 
328  da** u_internal=new da*[NUM_CHILDREN]; // in vector internal
329  da** v_internal=new da*[NUM_CHILDREN]; // in vector internal
330  T** pt_internal=new T*[NUM_CHILDREN];
331 
332 
333  for(unsigned int child=0;child<NUM_CHILDREN;child++)
334  {
335  u_internal[child]=NULL;
336  v_internal[child]=NULL;
337  pt_internal[child]=NULL;
338 
339  if(bCounts[child]!=0)
340  {
341  u_internal[child]=new da[binOp::getNextHighestPowerOfTwo(bCounts[child])];
342  v_internal[child]=new da[binOp::getNextHighestPowerOfTwo(bCounts[child])];
343  pt_internal[child]=new T[binOp::getNextHighestPowerOfTwo(bCounts[child])];
344  }
345 
346  }
347 
348 #ifdef MATVEC_PROFILE
349  dendro::timer::sfcmatvec::t_malloc.stop();
350 #endif
351 
352  DendroIntL counts[NUM_CHILDREN];
353  for(unsigned int child=0;child<NUM_CHILDREN;child++)
354  counts[child]=0;
355 
356 #ifdef MATVEC_PROFILE
357  dendro::timer::sfcmatvec::t_pts_p1_cpy.start();
358 #endif
359 
360  //pass 1
361  for(unsigned int p=pBegin;p<pEnd;p++)
362  {
363 
364  ijk[2]=0;
365  ijk[1]=0;
366  ijk[0]=0;
367 
368  if(pt_coords[p].zint()>pMid[2]) ijk[2]=2;
369  if(pt_coords[p].yint()>pMid[1]) ijk[1]=2;
370  if(pt_coords[p].xint()>pMid[0]) ijk[0]=2;
371 
372  if(pt_coords[p].zint()<pMid[2]) ijk[2]=0;
373  if(pt_coords[p].yint()<pMid[1]) ijk[1]=0;
374  if(pt_coords[p].xint()<pMid[0]) ijk[0]=0;
375 
376  if(pt_coords[p].zint()==pMid[2]) ijk[2]=1;
377  if(pt_coords[p].yint()==pMid[1]) ijk[1]=1;
378  if(pt_coords[p].xint()==pMid[0]) ijk[0]=1;
379 
380 
381  if(ijk[0]!=1 && ijk[1]!=1 && ijk[2]!=1)
382  {
383  cnum=(((ijk[2]>>1u)<<2u) | ((ijk[1]>>1u)<<1u) | ((ijk[0]>>1u)));
384  u_internal[cnum][counts[cnum]]=in[p];
385  pt_internal[cnum][counts[cnum]]=pt_coords[p];
386  counts[cnum]++;
387  }
388 
389 
390  }
391 #ifdef MATVEC_PROFILE
392  dendro::timer::sfcmatvec::t_pts_p1_cpy.stop();
393 #endif
394 
395 #ifdef MATVEC_PROFILE
396  dendro::timer::sfcmatvec::t_pts_p2_cpy.start();
397 #endif
398 
399  // pass 2;
400  for(unsigned int index=0;index<duplicateCoordIndex.size();index++) {
401 
402  ijk[2] = 0;
403  ijk[1] = 0;
404  ijk[0] = 0;
405 
406  if (pt_coords[duplicateCoordIndex[index]].zint() > pMid[2]) ijk[2] = 2;
407  if (pt_coords[duplicateCoordIndex[index]].yint() > pMid[1]) ijk[1] = 2;
408  if (pt_coords[duplicateCoordIndex[index]].xint() > pMid[0]) ijk[0] = 2;
409 
410  if (pt_coords[duplicateCoordIndex[index]].zint() < pMid[2]) ijk[2] = 0;
411  if (pt_coords[duplicateCoordIndex[index]].yint() < pMid[1]) ijk[1] = 0;
412  if (pt_coords[duplicateCoordIndex[index]].xint() < pMid[0]) ijk[0] = 0;
413 
414  if (pt_coords[duplicateCoordIndex[index]].zint() == pMid[2]) ijk[2] = 1;
415  if (pt_coords[duplicateCoordIndex[index]].yint() == pMid[1]) ijk[1] = 1;
416  if (pt_coords[duplicateCoordIndex[index]].xint() == pMid[0]) ijk[0] = 1;
417 
418  for(unsigned int child=0;child<NUM_CHILDREN;child++)
419  {
420  cnum=SFC_MATVEC_BUCKET_TABLE[(ijk[2]*SFC_MVEC_TABLE_Rz+ijk[1]*SFC_MVEC_TABLE_Ry+ijk[0])][child];
421  if(cnum==UCHAR_MAX)
422  break;
423 
424  u_internal[cnum][counts[cnum]]=in[duplicateCoordIndex[index]];
425  pt_internal[cnum][counts[cnum]]=pt_coords[duplicateCoordIndex[index]];
426  counts[cnum]++;
427  }
428 
429  }
430 
431 #ifdef MATVEC_PROFILE
432  dendro::timer::sfcmatvec::t_pts_p2_cpy.stop();
433 #endif
434 
435 
436 #ifdef MATVEC_PROFILE
437  dendro::timer::sfcmatvec::t_malloc.start();
438 #endif
439  da * parentVal=new da[nPe];
440 
441 #ifdef MATVEC_PROFILE
442  dendro::timer::sfcmatvec::t_malloc.stop();
443 #endif
444 
445  unsigned int pNodeCount=0;
446  bool computeParent=false;
447 
448  for(unsigned int child=0;child<NUM_CHILDREN;child++)
449  {
450  if(bCounts[child]<nPe)
451  {
452  computeParent=true;
453  break;
454  }
455 
456 
457  }
458 
459 
460 #ifdef MATVEC_PROFILE
461  dendro::timer::sfcmatvec::t_parent_bucket.start();
462 #endif
463 
464  if(computeParent)
465  {
466 
467 
468 
469  // find parent nodal index and corresponding values
470  for(unsigned int p=pBegin;p<pEnd;p++)
471  {
472  if(fem::seq::computeIJK(pt_coords[p],pOctant,dx,logDx,ijk)==PT_INTERNAL_DG_NODE)
473  {
474 
475  parentVal[IDX(ijk[0],ijk[1],ijk[2])]=in[p];
476  pNodeCount++;
477 
478  }
479 
480  }
481  }
482 
483 #ifdef MATVEC_PROFILE
484  dendro::timer::sfcmatvec::t_parent_bucket.stop();
485 #endif
486 
487 
488 
489  for(unsigned int child=0;child<NUM_CHILDREN;child++)
490  {
491 
492 
493  counts[child]=0;
494 
495 
496  if(bCounts[child]>nPe)
497  { // this is a non-leaf node. hence recurse.
498  fem::seq::matvec(pt_internal[child],u_internal[child],0,bCounts[child],lev+1,maxDepth,order,childElem[child],refEl,v_internal[child],mVecType);
499  }else{
500  // only reached by leaf nodes.
501 
502 
503 
504  da* interp3D_out=NULL;
505 
506  da* interp2D_in=NULL;
507  da* interp2D_out=NULL;
508 
509  da* interp1D_in=NULL;
510  da* interp1D_out=NULL;
511 
512  da* u_unzip=new da[nPe];
513  da* v_unzip=new da[nPe];
514 
515  bool* unzipStatus=new bool[nPe];
516  bool interpReq[OCT_DIR_TOTAL];
517  bool interpCpy[OCT_DIR_TOTAL];
518 
519  for(unsigned int i=0;i<nPe;i++)
520  unzipStatus[i]=false;
521 
522  for(unsigned int i=0;i<OCT_DIR_TOTAL;i++)
523  {
524  interpReq[i]=false;
525  interpCpy[i]=false;
526  }
527 
528 
529 #ifdef MATVEC_PROFILE
530  dendro::timer::sfcmatvec::t_pts_bucket.start();
531 #endif
532 
533  for(unsigned int p=0;p<bCounts[child];p++)
534  {
535  v_internal[child][p]=(da)0;
536  pt_status=fem::seq::computeIJK(pt_internal[child][p],childElem[child],dxb2,logDxb2,ijk);
537  if(pt_status==PT_INTERNAL_DG_NODE)
538  {
539  u_unzip[IDX(ijk[0],ijk[1],ijk[2])]=u_internal[child][p];
540  unzipStatus[IDX(ijk[0],ijk[1],ijk[2])]=true;
541  }
542 
543 
544  }
545 
546 #ifdef MATVEC_PROFILE
547  dendro::timer::sfcmatvec::t_pts_bucket.stop();
548 #endif
549 
550 
551 
552  if(bCounts[child]!=nPe) {
553  // we need to figure out the missing nodes and perform interpolations.
554 
555 #ifdef MATVEC_PROFILE
556  dendro::timer::sfcmatvec::t_p2cInterp.start();
557 #endif
558 
559  interp3D_out=new da[nPe];
560  interp2D_in=new da[(order+1)*(order+1)];
561  interp2D_out=new da[(order+1)*(order+1)];
562 
563  interp1D_in=new da[(order+1)];
564  interp1D_out=new da[(order+1)];
565 
566 
567 
568 
569 
570  const unsigned int intepCheckIndex=1;
571  //check for missing nodes
572  //OCT_DIR_LEFT_DOWN
573  if(!unzipStatus[IDX(0,0,intepCheckIndex)])
574  {
575  interpReq[OCT_DIR_LEFT_DOWN]=true;
576  }
577 
578 
579  //OCT_DIR_LEFT_UP
580  if(!unzipStatus[IDX(0,order,intepCheckIndex)])
581  {
582  interpReq[OCT_DIR_LEFT_UP]=true;
583  }
584 
585 
586  //OCT_DIR_LEFT_BACK
587  if(!unzipStatus[IDX(0,intepCheckIndex,0)])
588  {
589  interpReq[OCT_DIR_LEFT_BACK]=true;
590  }
591 
592 
593  //OCT_DIR_LEFT_FRONT
594  if(!unzipStatus[IDX(0,intepCheckIndex,order)])
595  {
596  interpReq[OCT_DIR_LEFT_FRONT]=true;
597  }
598 
599 
600  //OCT_DIR_RIGHT_DOWN
601  if(!unzipStatus[IDX(order,0,intepCheckIndex)])
602  {
603  interpReq[OCT_DIR_RIGHT_DOWN]=true;
604  }
605 
606 
607  //OCT_DIR_RIGHT_UP
608  if(!unzipStatus[IDX(order,order,intepCheckIndex)])
609  {
610  interpReq[OCT_DIR_RIGHT_UP]=true;
611  }
612 
613 
614  //OCT_DIR_RIGHT_BACK
615  if(!unzipStatus[IDX(order,intepCheckIndex,0)])
616  {
617  interpReq[OCT_DIR_RIGHT_BACK]=true;
618  }
619 
620 
621  //OCT_DIR_RIGHT_FRONT
622  if(!unzipStatus[IDX(order,intepCheckIndex,order)])
623  {
624  interpReq[OCT_DIR_RIGHT_FRONT]=true;
625  }
626 
627 
628  //OCT_DIR_DOWN_BACK
629  if(!unzipStatus[IDX(intepCheckIndex,0,0)])
630  {
631  interpReq[OCT_DIR_DOWN_BACK]=true;
632  }
633 
634 
635  //OCT_DIR_UP_BACK
636  if(!unzipStatus[IDX(intepCheckIndex,order,0)])
637  {
638  interpReq[OCT_DIR_UP_BACK]=true;
639  }
640 
641 
642  //OCT_DIR_DOWN_FRONT
643  if(!unzipStatus[IDX(intepCheckIndex,0,order)])
644  {
645  interpReq[OCT_DIR_DOWN_FRONT]=true;
646  }
647 
648 
649  //OCT_DIR_UP_FRONT
650  if(!unzipStatus[IDX(intepCheckIndex,order,order)])
651  {
652  interpReq[OCT_DIR_UP_FRONT]=true;
653  }
654 
655  // check for faces.
656 
657  //OCT_DIR_LEFT
658  if(!unzipStatus[IDX(0,intepCheckIndex,intepCheckIndex)])
659  {
660  interpReq[OCT_DIR_LEFT]=true;
661 
662  }
663 
664 
665  //OCT_DIR_RIGT
666  if(!unzipStatus[IDX(order,intepCheckIndex,intepCheckIndex)])
667  {
668  interpReq[OCT_DIR_RIGHT]=true;
669 
670 
671  }
672 
673 
674  //OCT_DIR_DOWN
675  if(!unzipStatus[IDX(intepCheckIndex,0,intepCheckIndex)])
676  {
677  interpReq[OCT_DIR_DOWN]=true;
678 
679  }
680 
681 
682  // OCT_DIR_UP
683  if(!unzipStatus[IDX(intepCheckIndex,order,intepCheckIndex)])
684  {
685  interpReq[OCT_DIR_UP]=true;
686 
687  }
688 
689 
690  //OCT_DIR_BACK
691  if(!unzipStatus[IDX(intepCheckIndex,intepCheckIndex,0)])
692  {
693  interpReq[OCT_DIR_BACK]=true;
694 
695  }
696 
697 
698  //OCT_DIR_FRONT
699  if(!unzipStatus[IDX(intepCheckIndex,intepCheckIndex,order)])
700  {
701  interpReq[OCT_DIR_FRONT]=true;
702  }
703 
704 
705 
706 
707  assert(pNodeCount==nPe);
708  if(pNodeCount!=nPe) std::cout<<"[Error]"<<__func__<<" pNodeCount: "<<pNodeCount<<std::endl;
709 
710 
711  //face copy.
712  if(interpReq[OCT_DIR_LEFT])
713  {
714 
715  assert(binOp::getBit(child,0)==0);
716  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,1));
717  //std::cout<<"cnumL: "<<cnum<<" binOp::getBit(child,2) "<<binOp::getBit(child,2)<<" binOp::getBit(child,1)"<<binOp::getBit(child,1)<<std::endl;
718  for(unsigned int k=0;k<(order+1);k++)
719  for(unsigned int j=0;j<(order+1);j++)
720  interp2D_in[IDX2D(j,k)]=parentVal[IDX(0,j,k)];
721 
722  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
723 
724 
725  for(unsigned int k=0;k<(order+1);k++)
726  for(unsigned int j=0;j<(order+1);j++)
727  u_unzip[IDX(0,j,k)]=interp2D_out[IDX2D(j,k)];
728 
729 
730  interpCpy[OCT_DIR_LEFT_DOWN]=true;
731  interpCpy[OCT_DIR_LEFT_UP]=true;
732  interpCpy[OCT_DIR_LEFT_BACK]=true;
733  interpCpy[OCT_DIR_LEFT_FRONT]=true;
734 
735  }
736 
737 
738  if(interpReq[OCT_DIR_RIGHT])
739  {
740 
741  assert(binOp::getBit(child,0)==1);
742  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,1));
743 
744  for(unsigned int k=0;k<(order+1);k++)
745  for(unsigned int j=0;j<(order+1);j++)
746  interp2D_in[IDX2D(j,k)]=parentVal[IDX(order,j,k)];
747 
748  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
749 
750  for(unsigned int k=0;k<(order+1);k++)
751  for(unsigned int j=0;j<(order+1);j++)
752  u_unzip[IDX(order,j,k)]=interp2D_out[IDX2D(j,k)];
753 
754  interpCpy[OCT_DIR_RIGHT_DOWN]=true;
755  interpCpy[OCT_DIR_RIGHT_UP]=true;
756  interpCpy[OCT_DIR_RIGHT_BACK]=true;
757  interpCpy[OCT_DIR_RIGHT_FRONT]=true;
758 
759  }
760 
761 
762  if(interpReq[OCT_DIR_DOWN])
763  {
764 
765  assert(binOp::getBit(child,1)==0);
766  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,0));
767 
768 
769  for(unsigned int k=0;k<(order+1);k++)
770  for(unsigned int i=0;i<(order+1);i++)
771  interp2D_in[IDX2D(i,k)]=parentVal[IDX(i,0,k)];
772 
773  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
774 
775  for(unsigned int k=0;k<(order+1);k++)
776  for(unsigned int i=0;i<(order+1);i++)
777  u_unzip[IDX(i,0,k)]=interp2D_out[IDX2D(i,k)];
778 
779  interpCpy[OCT_DIR_RIGHT_DOWN]=true;
780  interpCpy[OCT_DIR_LEFT_DOWN]=true;
781  interpCpy[OCT_DIR_DOWN_BACK]=true;
782  interpCpy[OCT_DIR_DOWN_FRONT]=true;
783 
784  }
785 
786 
787  if(interpReq[OCT_DIR_UP])
788  {
789 
790  assert(binOp::getBit(child,1)==1);
791  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,0));
792 
793 
794  for(unsigned int k=0;k<(order+1);k++)
795  for(unsigned int i=0;i<(order+1);i++)
796  interp2D_in[IDX2D(i,k)]=parentVal[IDX(i,order,k)];
797 
798  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
799 
800  for(unsigned int k=0;k<(order+1);k++)
801  for(unsigned int i=0;i<(order+1);i++)
802  u_unzip[IDX(i,order,k)]=interp2D_out[IDX2D(i,k)];
803 
804 
805  interpCpy[OCT_DIR_RIGHT_UP]=true;
806  interpCpy[OCT_DIR_LEFT_UP]=true;
807  interpCpy[OCT_DIR_UP_BACK]=true;
808  interpCpy[OCT_DIR_UP_FRONT]=true;
809 
810  }
811 
812 
813  if(interpReq[OCT_DIR_BACK])
814  {
815 
816  assert(binOp::getBit(child,2)==0);
817  cnum=(binOp::getBit(child,1)<<1)|(binOp::getBit(child,0));
818 
819  for(unsigned int j=0;j<(order+1);j++)
820  for(unsigned int i=0;i<(order+1);i++)
821  interp2D_in[IDX2D(i,j)]=parentVal[IDX(i,j,0)];
822 
823  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
824 
825 
826  for(unsigned int j=0;j<(order+1);j++)
827  for(unsigned int i=0;i<(order+1);i++)
828  u_unzip[IDX(i,j,0)]=interp2D_out[IDX2D(i,j)];
829 
830  interpCpy[OCT_DIR_RIGHT_BACK]=true;
831  interpCpy[OCT_DIR_LEFT_BACK]=true;
832  interpCpy[OCT_DIR_UP_BACK]=true;
833  interpCpy[OCT_DIR_DOWN_BACK]=true;
834 
835  }
836 
837 
838  if(interpReq[OCT_DIR_FRONT])
839  {
840 
841  assert(binOp::getBit(child,2)==1);
842  cnum=(binOp::getBit(child,1)<<1)|(binOp::getBit(child,0));
843 
844  for(unsigned int j=0;j<(order+1);j++)
845  for(unsigned int i=0;i<(order+1);i++)
846  interp2D_in[IDX2D(i,j)]=parentVal[IDX(i,j,order)];
847 
848  refEl->I2D_Parent2Child(interp2D_in,interp2D_out,cnum);
849 
850 
851  for(unsigned int j=0;j<(order+1);j++)
852  for(unsigned int i=0;i<(order+1);i++)
853  u_unzip[IDX(i,j,order)]=interp2D_out[IDX2D(i,j)];
854 
855 
856  interpCpy[OCT_DIR_RIGHT_FRONT]=true;
857  interpCpy[OCT_DIR_LEFT_FRONT]=true;
858  interpCpy[OCT_DIR_UP_FRONT]=true;
859  interpCpy[OCT_DIR_DOWN_FRONT]=true;
860 
861  }
862 
863  //edges copy.
864  if( interpReq[OCT_DIR_LEFT_DOWN] && (!interpCpy[OCT_DIR_LEFT_DOWN]))
865  {
866 
867  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,1)==0);
868  cnum=binOp::getBit(child,2);
869 
870  for(unsigned int k=0;k<(order+1);k++)
871  interp1D_in[k]=parentVal[IDX(0,0,k)];
872 
873  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
874 
875  for(unsigned int k=0;k<(order+1);k++)
876  u_unzip[IDX(0,0,k)]=interp1D_out[k];
877 
878  }
879 
880  if( interpReq[OCT_DIR_LEFT_UP] && (!interpCpy[OCT_DIR_LEFT_UP]))
881  {
882 
883  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,1)==1);
884  cnum=binOp::getBit(child,2);
885 
886  for(unsigned int k=0;k<(order+1);k++)
887  interp1D_in[k]=parentVal[IDX(0,order,k)];
888 
889  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
890 
891  for(unsigned int k=0;k<(order+1);k++)
892  u_unzip[IDX(0,order,k)]=interp1D_out[k];
893 
894  }
895 
896  if(interpReq[OCT_DIR_LEFT_BACK] && (!interpCpy[OCT_DIR_LEFT_BACK]))
897  {
898 
899  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,2)==0);
900  cnum=binOp::getBit(child,1);
901 
902  for(unsigned int j=0;j<(order+1);j++)
903  interp1D_in[j]=parentVal[IDX(0,j,0)];
904 
905  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
906 
907  for(unsigned int j=0;j<(order+1);j++)
908  u_unzip[IDX(0,j,0)]=interp1D_out[j];
909 
910 
911  }
912 
913  if(interpReq[OCT_DIR_LEFT_FRONT] && (!interpCpy[OCT_DIR_LEFT_FRONT]))
914  {
915 
916  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,2)==1);
917  cnum=binOp::getBit(child,1);
918 
919  for(unsigned int j=0;j<(order+1);j++)
920  interp1D_in[j]=parentVal[IDX(0,j,order)];
921 
922  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
923 
924 
925  for(unsigned int j=0;j<(order+1);j++)
926  u_unzip[IDX(0,j,order)]=interp1D_out[j];
927 
928 
929  }
930 
931 
932 
933  if(interpReq[OCT_DIR_RIGHT_DOWN] && (!interpCpy[OCT_DIR_RIGHT_DOWN]))
934  {
935 
936  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,1)==0);
937  cnum=binOp::getBit(child,2);
938 
939  for(unsigned int k=0;k<(order+1);k++)
940  interp1D_in[k]=parentVal[IDX(order,0,k)];
941 
942  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
943 
944  for(unsigned int k=0;k<(order+1);k++)
945  u_unzip[IDX(order,0,k)]=interp1D_out[k];
946 
947  }
948 
949  if(interpReq[OCT_DIR_RIGHT_UP] && (!interpCpy[OCT_DIR_RIGHT_UP]))
950  {
951 
952  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,1)==1);
953  cnum=binOp::getBit(child,2);
954 
955  for(unsigned int k=0;k<(order+1);k++)
956  interp1D_in[k]=parentVal[IDX(order,order,k)];
957 
958  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
959 
960  for(unsigned int k=0;k<(order+1);k++)
961  u_unzip[IDX(order,order,k)]=interp1D_out[k];
962 
963  }
964 
965  if(interpReq[OCT_DIR_RIGHT_BACK] && (!interpCpy[OCT_DIR_RIGHT_BACK]))
966  {
967  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,2)==0);
968  cnum=binOp::getBit(child,1);
969 
970  for(unsigned int j=0;j<(order+1);j++)
971  interp1D_in[j]=parentVal[IDX(order,j,0)];
972 
973  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
974 
975 
976  for(unsigned int j=0;j<(order+1);j++)
977  u_unzip[IDX(order,j,0)]=interp1D_out[j];
978 
979  }
980 
981  if(interpReq[OCT_DIR_RIGHT_FRONT] && (!interpCpy[OCT_DIR_RIGHT_FRONT]))
982  {
983 
984  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,2)==1);
985  cnum=binOp::getBit(child,1);
986 
987  for(unsigned int j=0;j<(order+1);j++)
988  interp1D_in[j]=parentVal[IDX(order,j,order)];
989 
990  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
991 
992  for(unsigned int j=0;j<(order+1);j++)
993  u_unzip[IDX(order,j,order)]=interp1D_out[j];
994 
995  }
996 
997 
998  if(interpReq[OCT_DIR_DOWN_BACK] && (!interpCpy[OCT_DIR_DOWN_BACK]))
999  {
1000 
1001  assert(binOp::getBit(child,1)==0 && binOp::getBit(child,2)==0);
1002  cnum=binOp::getBit(child,0);
1003 
1004  for(unsigned int i=0;i<(order+1);i++)
1005  interp1D_in[i]=parentVal[IDX(i,0,0)];
1006 
1007  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
1008 
1009  for(unsigned int i=0;i<(order+1);i++)
1010  u_unzip[IDX(i,0,0)]=interp1D_out[i];
1011 
1012 
1013  }
1014 
1015  if(interpReq[OCT_DIR_DOWN_FRONT] && (!interpCpy[OCT_DIR_DOWN_FRONT]))
1016  {
1017 
1018  assert(binOp::getBit(child,1)==0 && binOp::getBit(child,2)==1);
1019  cnum=binOp::getBit(child,0);
1020 
1021  for(unsigned int i=0;i<(order+1);i++)
1022  interp1D_in[i]=parentVal[IDX(i,0,order)];
1023 
1024  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
1025 
1026  for(unsigned int i=0;i<(order+1);i++)
1027  u_unzip[IDX(i,0,order)]=interp1D_out[i];
1028 
1029  }
1030 
1031  if(interpReq[OCT_DIR_UP_BACK] && (!interpCpy[OCT_DIR_UP_BACK]))
1032  {
1033 
1034  assert(binOp::getBit(child,1)==1 && binOp::getBit(child,2)==0);
1035  cnum=binOp::getBit(child,0);
1036 
1037  for(unsigned int i=0;i<(order+1);i++)
1038  interp1D_in[i]=parentVal[IDX(i,order,0)];
1039 
1040  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
1041 
1042  for(unsigned int i=0;i<(order+1);i++)
1043  u_unzip[IDX(i,order,0)]=interp1D_out[i];
1044 
1045  }
1046 
1047  if(interpReq[OCT_DIR_UP_FRONT] && (!interpCpy[OCT_DIR_UP_FRONT]))
1048  {
1049 
1050  assert(binOp::getBit(child,1)==1 && binOp::getBit(child,2)==1);
1051  cnum=binOp::getBit(child,0);
1052 
1053  for(unsigned int i=0;i<(order+1);i++)
1054  interp1D_in[i]=parentVal[IDX(i,order,order)];
1055 
1056  refEl->I1D_Parent2Child(interp1D_in,interp1D_out,cnum);
1057 
1058  for(unsigned int i=0;i<(order+1);i++)
1059  u_unzip[IDX(i,order,order)]=interp1D_out[i];
1060 
1061  }
1062 
1063 #ifdef MATVEC_PROFILE
1064  dendro::timer::sfcmatvec::t_p2cInterp.stop();
1065 #endif
1066 
1067  }
1068 
1069 
1070  if(mVecType==0) // mass matvec
1071  fem::operators::poisson::elementMassMatvec(u_unzip,childElem[child],refEl,v_unzip);
1072  else if(mVecType==1) // stifness matvec
1073  fem::operators::poisson::elementStiffnessMatvec(u_unzip,childElem[child],refEl,v_unzip);
1074 
1075 
1076 
1077  bool * accumStatus = new bool [nPe];
1078  bool * c2pInterp=new bool[nPe];
1079 
1080  for(unsigned int node=0;node<nPe;node++)
1081  {
1082  accumStatus[node]=false;
1083  c2pInterp[node]=false;
1084  }
1085 
1086 
1087  if(bCounts[child]!=nPe)
1088  {
1089 
1090 #ifdef MATVEC_PROFILE
1091  dendro::timer::sfcmatvec::t_c2pInterp.start();
1092 #endif
1093 
1094  unsigned int cnum;
1095 
1096  for(unsigned int node=0;node<nPe;node++)
1097  parentVal[node]=(da)0;
1098 
1099  bool vertexAcc[NUM_CHILDREN];
1100  for(unsigned int i=0;i<NUM_CHILDREN;i++)
1101  vertexAcc[i]=false;
1102 
1103  //face copy.
1104  if(interpReq[OCT_DIR_LEFT])
1105  {
1106 
1107  assert(binOp::getBit(child,0)==0);
1108  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,1));
1109  //std::cout<<"cnumL: "<<cnum<<" binOp::getBit(child,2) "<<binOp::getBit(child,2)<<" binOp::getBit(child,1)"<<binOp::getBit(child,1)<<std::endl;
1110  for(unsigned int k=0;k<(order+1);k++)
1111  for(unsigned int j=0;j<(order+1);j++)
1112  interp2D_in[IDX2D(j,k)]=v_unzip[IDX(0,j,k)];
1113 
1114  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1115 
1116 
1117  for(unsigned int k=1;k<(order);k++)
1118  for(unsigned int j=1;j<(order);j++)
1119  {
1120  parentVal[IDX(0,j,k)]=interp2D_out[IDX2D(j,k)];
1121  c2pInterp[IDX(0,j,k)]=true;
1122  }
1123 
1124 
1125  }
1126 
1127 
1128  if(interpReq[OCT_DIR_RIGHT])
1129  {
1130 
1131  assert(binOp::getBit(child,0)==1);
1132  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,1));
1133 
1134  for(unsigned int k=0;k<(order+1);k++)
1135  for(unsigned int j=0;j<(order+1);j++)
1136  interp2D_in[IDX2D(j,k)]=v_unzip[IDX(order,j,k)];
1137 
1138  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1139 
1140  for(unsigned int k=1;k<(order);k++)
1141  for(unsigned int j=1;j<(order);j++)
1142  {
1143  parentVal[IDX(order,j,k)]=interp2D_out[IDX2D(j,k)];
1144  c2pInterp[IDX(order,j,k)]=true;
1145  }
1146 
1147 
1148  }
1149 
1150 
1151  if(interpReq[OCT_DIR_DOWN])
1152  {
1153 
1154  assert(binOp::getBit(child,1)==0);
1155  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,0));
1156 
1157 
1158  for(unsigned int k=0;k<(order+1);k++)
1159  for(unsigned int i=0;i<(order+1);i++)
1160  interp2D_in[IDX2D(i,k)]=v_unzip[IDX(i,0,k)];
1161 
1162  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1163 
1164  for(unsigned int k=1;k<(order);k++)
1165  for(unsigned int i=1;i<(order);i++)
1166  {
1167  parentVal[IDX(i,0,k)]=interp2D_out[IDX2D(i,k)];
1168  c2pInterp[IDX(i,0,k)]=true;
1169  }
1170 
1171 
1172 
1173  }
1174 
1175 
1176  if(interpReq[OCT_DIR_UP])
1177  {
1178 
1179  assert(binOp::getBit(child,1)==1);
1180  cnum=(binOp::getBit(child,2)<<1)|(binOp::getBit(child,0));
1181 
1182 
1183  for(unsigned int k=0;k<(order+1);k++)
1184  for(unsigned int i=0;i<(order+1);i++)
1185  interp2D_in[IDX2D(i,k)]=v_unzip[IDX(i,order,k)];
1186 
1187  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1188 
1189  for(unsigned int k=1;k<(order);k++)
1190  for(unsigned int i=1;i<(order);i++)
1191  {
1192  parentVal[IDX(i,order,k)]=interp2D_out[IDX2D(i,k)];
1193  c2pInterp[IDX(i,order,k)]=true;
1194  }
1195 
1196 
1197 
1198  }
1199 
1200 
1201  if(interpReq[OCT_DIR_BACK])
1202  {
1203 
1204  assert(binOp::getBit(child,2)==0);
1205  cnum=(binOp::getBit(child,1)<<1)|(binOp::getBit(child,0));
1206 
1207  for(unsigned int j=0;j<(order+1);j++)
1208  for(unsigned int i=0;i<(order+1);i++)
1209  interp2D_in[IDX2D(i,j)]=v_unzip[IDX(i,j,0)];
1210 
1211  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1212 
1213 
1214  for(unsigned int j=1;j<(order);j++)
1215  for(unsigned int i=1;i<(order);i++)
1216  {
1217  parentVal[IDX(i,j,0)]=interp2D_out[IDX2D(i,j)];
1218  c2pInterp[IDX(i,j,0)]=true;
1219  }
1220 
1221 
1222 
1223 
1224  }
1225 
1226 
1227  if(interpReq[OCT_DIR_FRONT])
1228  {
1229 
1230  assert(binOp::getBit(child,2)==1);
1231  cnum=(binOp::getBit(child,1)<<1)|(binOp::getBit(child,0));
1232 
1233  for(unsigned int j=0;j<(order+1);j++)
1234  for(unsigned int i=0;i<(order+1);i++)
1235  interp2D_in[IDX2D(i,j)]=v_unzip[IDX(i,j,order)];
1236 
1237  refEl->I2D_Child2Parent(interp2D_in,interp2D_out,cnum);
1238 
1239 
1240  for(unsigned int j=1;j<(order);j++)
1241  for(unsigned int i=1;i<(order);i++)
1242  {
1243  parentVal[IDX(i,j,order)]=interp2D_out[IDX2D(i,j)];
1244  c2pInterp[IDX(i,j,order)]=true;
1245  }
1246 
1247 
1248 
1249 
1250  }
1251 
1252  //edges copy.
1253  if( interpReq[OCT_DIR_LEFT_DOWN])
1254  {
1255 
1256  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,1)==0);
1257  cnum=binOp::getBit(child,2);
1258 
1259  for(unsigned int k=0;k<(order+1);k++)
1260  interp1D_in[k]=v_unzip[IDX(0,0,k)];
1261 
1262  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1263 
1264  for(unsigned int k=1;k<(order);k++)
1265  {
1266  parentVal[IDX(0,0,k)]=interp1D_out[k];
1267  c2pInterp[IDX(0,0,k)]=true;
1268  }
1269 
1270 
1271  if(!vertexAcc[0])
1272  {
1273  parentVal[IDX(0,0,0)]=interp1D_out[0];
1274  c2pInterp[IDX(0,0,0)]=true;
1275  }
1276 
1277 
1278  if(!vertexAcc[4])
1279  {
1280  parentVal[IDX(0,0,order)]=interp1D_out[order];
1281  c2pInterp[IDX(0,0,order)]=true;
1282  }
1283 
1284 
1285  vertexAcc[0]=true;
1286  vertexAcc[4]=true;
1287 
1288  }
1289 
1290  if( interpReq[OCT_DIR_LEFT_UP])
1291  {
1292 
1293  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,1)==1);
1294  cnum=binOp::getBit(child,2);
1295 
1296  for(unsigned int k=0;k<(order+1);k++)
1297  interp1D_in[k]=v_unzip[IDX(0,order,k)];
1298 
1299  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1300 
1301  for(unsigned int k=1;k<(order);k++)
1302  {
1303  parentVal[IDX(0,order,k)]=interp1D_out[k];
1304  c2pInterp[IDX(0,order,k)]=true;
1305  }
1306 
1307 
1308  if(!vertexAcc[2])
1309  {
1310  parentVal[IDX(0,order,0)]=interp1D_out[0];
1311  c2pInterp[IDX(0,order,0)]=true;
1312  }
1313 
1314 
1315  if(!vertexAcc[6])
1316  {
1317  parentVal[IDX(0,order,order)]=interp1D_out[order];
1318  c2pInterp[IDX(0,order,order)]=true;
1319  }
1320 
1321 
1322  vertexAcc[2]=true;
1323  vertexAcc[6]=true;
1324 
1325  }
1326 
1327  if(interpReq[OCT_DIR_LEFT_BACK])
1328  {
1329 
1330  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,2)==0);
1331  cnum=binOp::getBit(child,1);
1332 
1333  for(unsigned int j=0;j<(order+1);j++)
1334  interp1D_in[j]=v_unzip[IDX(0,j,0)];
1335 
1336  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1337 
1338  for(unsigned int j=1;j<(order);j++)
1339  {
1340  parentVal[IDX(0,j,0)]=interp1D_out[j];
1341  c2pInterp[IDX(0,j,0)]=true;
1342  }
1343 
1344 
1345 
1346  if(!vertexAcc[0])
1347  {
1348  parentVal[IDX(0,0,0)]=interp1D_out[0];
1349  c2pInterp[IDX(0,0,0)]=true;
1350  }
1351 
1352 
1353  if(!vertexAcc[2])
1354  {
1355  parentVal[IDX(0,order,0)]=interp1D_out[order];
1356  c2pInterp[IDX(0,order,0)]=true;
1357  }
1358 
1359 
1360 
1361  vertexAcc[0]=true;
1362  vertexAcc[2]=true;
1363 
1364  }
1365 
1366  if(interpReq[OCT_DIR_LEFT_FRONT])
1367  {
1368 
1369  assert(binOp::getBit(child,0)==0 && binOp::getBit(child,2)==1);
1370  cnum=binOp::getBit(child,1);
1371 
1372  for(unsigned int j=0;j<(order+1);j++)
1373  interp1D_in[j]=v_unzip[IDX(0,j,order)];
1374 
1375  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1376 
1377 
1378  for(unsigned int j=1;j<(order);j++)
1379  {
1380  parentVal[IDX(0,j,order)]=interp1D_out[j];
1381  c2pInterp[IDX(0,j,order)]=true;
1382  }
1383 
1384 
1385  if(!vertexAcc[4])
1386  {
1387  parentVal[IDX(0,0,order)]=interp1D_out[0];
1388  c2pInterp[IDX(0,0,order)]=true;
1389  }
1390 
1391 
1392  if(!vertexAcc[6])
1393  {
1394  parentVal[IDX(0,order,order)]=interp1D_out[order];
1395  c2pInterp[IDX(0,order,order)]=true;
1396  }
1397 
1398 
1399  vertexAcc[4]=true;
1400  vertexAcc[6]=true;
1401 
1402  }
1403 
1404 
1405 
1406  if(interpReq[OCT_DIR_RIGHT_DOWN])
1407  {
1408 
1409  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,1)==0);
1410  cnum=binOp::getBit(child,2);
1411 
1412  for(unsigned int k=0;k<(order+1);k++)
1413  interp1D_in[k]=v_unzip[IDX(order,0,k)];
1414 
1415  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1416 
1417  for(unsigned int k=1;k<(order);k++)
1418  {
1419  parentVal[IDX(order,0,k)]=interp1D_out[k];
1420  c2pInterp[IDX(order,0,k)]=true;
1421  }
1422 
1423 
1424 
1425  if(!vertexAcc[1])
1426  {
1427  parentVal[IDX(order,0,0)]=interp1D_out[0];
1428  c2pInterp[IDX(order,0,0)]=true;
1429  }
1430 
1431 
1432  if(!vertexAcc[5])
1433  {
1434  parentVal[IDX(order,0,order)]=interp1D_out[order];
1435  c2pInterp[IDX(order,0,order)]=true;
1436  }
1437 
1438 
1439  vertexAcc[1]=true;
1440  vertexAcc[5]=true;
1441 
1442  }
1443 
1444  if(interpReq[OCT_DIR_RIGHT_UP])
1445  {
1446 
1447  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,1)==1);
1448  cnum=binOp::getBit(child,2);
1449 
1450  for(unsigned int k=0;k<(order+1);k++)
1451  interp1D_in[k]=v_unzip[IDX(order,order,k)];
1452 
1453  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1454 
1455  for(unsigned int k=1;k<order;k++)
1456  {
1457  parentVal[IDX(order,order,k)]=interp1D_out[k];
1458  c2pInterp[IDX(order,order,k)]=true;
1459  }
1460 
1461 
1462 
1463  if(!vertexAcc[3])
1464  {
1465  parentVal[IDX(order,order,0)]=interp1D_out[0];
1466  c2pInterp[IDX(order,order,0)]=true;
1467  }
1468 
1469 
1470  if(!vertexAcc[7])
1471  {
1472  parentVal[IDX(order,order,order)]=interp1D_out[order];
1473  c2pInterp[IDX(order,order,order)]=true;
1474  }
1475 
1476 
1477 
1478  vertexAcc[3]=true;
1479  vertexAcc[7]=true;
1480 
1481  }
1482 
1483  if(interpReq[OCT_DIR_RIGHT_BACK])
1484  {
1485  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,2)==0);
1486  cnum=binOp::getBit(child,1);
1487 
1488  for(unsigned int j=0;j<(order+1);j++)
1489  interp1D_in[j]=v_unzip[IDX(order,j,0)];
1490 
1491  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1492 
1493 
1494  for(unsigned int j=1;j<(order);j++)
1495  {
1496  parentVal[IDX(order,j,0)]=interp1D_out[j];
1497  c2pInterp[IDX(order,j,0)]=true;
1498  }
1499 
1500 
1501  if(!vertexAcc[1])
1502  {
1503  parentVal[IDX(order,0,0)]=interp1D_out[0];
1504  c2pInterp[IDX(order,0,0)]=true;
1505  }
1506 
1507 
1508  if(!vertexAcc[3])
1509  {
1510  parentVal[IDX(order,order,0)]=interp1D_out[order];
1511  c2pInterp[IDX(order,order,0)]=true;
1512  }
1513 
1514 
1515 
1516  vertexAcc[1]=true;
1517  vertexAcc[3]=true;
1518 
1519 
1520  }
1521 
1522  if(interpReq[OCT_DIR_RIGHT_FRONT])
1523  {
1524 
1525  assert(binOp::getBit(child,0)==1 && binOp::getBit(child,2)==1);
1526  cnum=binOp::getBit(child,1);
1527 
1528  for(unsigned int j=0;j<(order+1);j++)
1529  interp1D_in[j]=v_unzip[IDX(order,j,order)];
1530 
1531  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1532 
1533  for(unsigned int j=1;j<(order);j++)
1534  {
1535  parentVal[IDX(order,j,order)]=interp1D_out[j];
1536  c2pInterp[IDX(order,j,order)]=true;
1537  }
1538 
1539 
1540  if(!vertexAcc[5])
1541  {
1542  parentVal[IDX(order,0,order)]=interp1D_out[0];
1543  c2pInterp[IDX(order,0,order)]=true;
1544  }
1545 
1546 
1547  if(!vertexAcc[7])
1548  {
1549  parentVal[IDX(order,order,order)]=interp1D_out[order];
1550  c2pInterp[IDX(order,order,order)]=true;
1551  }
1552 
1553 
1554 
1555  vertexAcc[5]=true;
1556  vertexAcc[7]=true;
1557 
1558  }
1559 
1560 
1561  if(interpReq[OCT_DIR_DOWN_BACK])
1562  {
1563 
1564  assert(binOp::getBit(child,1)==0 && binOp::getBit(child,2)==0);
1565  cnum=binOp::getBit(child,0);
1566 
1567  for(unsigned int i=0;i<(order+1);i++)
1568  interp1D_in[i]=v_unzip[IDX(i,0,0)];
1569 
1570  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1571 
1572  for(unsigned int i=1;i<(order);i++)
1573  {
1574  parentVal[IDX(i,0,0)]=interp1D_out[i];
1575  c2pInterp[IDX(i,0,0)]=true;
1576  }
1577 
1578 
1579  if(!vertexAcc[0])
1580  {
1581  parentVal[IDX(0,0,0)]=interp1D_out[0];
1582  c2pInterp[IDX(0,0,0)]=true;
1583  }
1584 
1585 
1586  if(!vertexAcc[1])
1587  {
1588  parentVal[IDX(order,0,0)]=interp1D_out[order];
1589  c2pInterp[IDX(order,0,0)]=true;
1590  }
1591 
1592 
1593 
1594  vertexAcc[0]=true;
1595  vertexAcc[1]=true;
1596 
1597  }
1598 
1599  if(interpReq[OCT_DIR_DOWN_FRONT])
1600  {
1601 
1602  assert(binOp::getBit(child,1)==0 && binOp::getBit(child,2)==1);
1603  cnum=binOp::getBit(child,0);
1604 
1605  for(unsigned int i=0;i<(order+1);i++)
1606  interp1D_in[i]=v_unzip[IDX(i,0,order)];
1607 
1608  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1609 
1610  for(unsigned int i=1;i<(order);i++)
1611  {
1612  parentVal[IDX(i,0,order)]=interp1D_out[i];
1613  c2pInterp[IDX(i,0,order)]=true;
1614  }
1615 
1616 
1617  if(!vertexAcc[4])
1618  {
1619  parentVal[IDX(0,0,order)]=interp1D_out[0];
1620  c2pInterp[IDX(0,0,order)]=true;
1621  }
1622 
1623 
1624  if(!vertexAcc[5])
1625  {
1626  parentVal[IDX(order,0,order)]=interp1D_out[order];
1627  c2pInterp[IDX(order,0,order)]=true;
1628  }
1629 
1630 
1631 
1632  vertexAcc[4]=true;
1633  vertexAcc[5]=true;
1634 
1635 
1636  }
1637 
1638  if(interpReq[OCT_DIR_UP_BACK])
1639  {
1640 
1641  assert(binOp::getBit(child,1)==1 && binOp::getBit(child,2)==0);
1642  cnum=binOp::getBit(child,0);
1643 
1644  for(unsigned int i=0;i<(order+1);i++)
1645  interp1D_in[i]=v_unzip[IDX(i,order,0)];
1646 
1647  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1648 
1649  for(unsigned int i=1;i<(order);i++)
1650  {
1651  parentVal[IDX(i,order,0)]=interp1D_out[i];
1652  c2pInterp[IDX(i,order,0)]=true;
1653  }
1654 
1655 
1656 
1657  if(!vertexAcc[2])
1658  {
1659  parentVal[IDX(0,order,0)]=interp1D_out[0];
1660  c2pInterp[IDX(0,order,0)]=true;
1661  }
1662 
1663 
1664 
1665  if(!vertexAcc[3])
1666  {
1667  parentVal[IDX(order,order,0)]=interp1D_out[order];
1668  c2pInterp[IDX(order,order,0)]=true;
1669  }
1670 
1671 
1672 
1673  vertexAcc[2]=true;
1674  vertexAcc[3]=true;
1675 
1676 
1677  }
1678 
1679  if(interpReq[OCT_DIR_UP_FRONT])
1680  {
1681 
1682  assert(binOp::getBit(child,1)==1 && binOp::getBit(child,2)==1);
1683  cnum=binOp::getBit(child,0);
1684 
1685  for(unsigned int i=0;i<(order+1);i++)
1686  interp1D_in[i]=v_unzip[IDX(i,order,order)];
1687 
1688  refEl->I1D_Child2Parent(interp1D_in,interp1D_out,cnum);
1689 
1690 
1691  for(unsigned int i=1;i<(order);i++)
1692  {
1693  parentVal[IDX(i,order,order)]=interp1D_out[i];
1694  c2pInterp[IDX(i,order,order)]=true;
1695  }
1696 
1697 
1698 
1699 
1700  if(!vertexAcc[6])
1701  {
1702  parentVal[IDX(0,order,order)]=interp1D_out[0];
1703  c2pInterp[IDX(0,order,order)]=true;
1704  }
1705 
1706 
1707  if(!vertexAcc[7])
1708  {
1709  parentVal[IDX(order,order,order)]=interp1D_out[order];
1710  c2pInterp[IDX(order,order,order)]=true;
1711  }
1712 
1713 
1714 
1715  vertexAcc[6]=true;
1716  vertexAcc[7]=true;
1717 
1718 
1719  }
1720 
1721 
1722 
1723  /*for(unsigned int node=0;node<nPe;node++)
1724  out[parentIndex[node]]+=parentVal[node];*/
1725 
1726 
1727 #ifdef MATVEC_PROFILE
1728  dendro::timer::sfcmatvec::t_c2pInterp.stop();
1729 #endif
1730 
1731 
1732 #ifdef MATVEC_PROFILE
1733  dendro::timer::sfcmatvec::t_accum.start();
1734 #endif
1735  for(unsigned int p=pBegin;p<pEnd;p++) {
1736  if (fem::seq::computeIJK(pt_coords[p], pOctant, dx, logDx, ijk) == PT_INTERNAL_DG_NODE)
1737  {
1738  out[p]+=parentVal[IDX(ijk[0],ijk[1],ijk[2])];
1739  if(c2pInterp[IDX(ijk[0],ijk[1],ijk[2])])accumStatus[IDX(ijk[0],ijk[1],ijk[2])]=true;
1740  }
1741 
1742  }
1743 
1744 #ifdef MATVEC_PROFILE
1745  dendro::timer::sfcmatvec::t_accum.stop();
1746 #endif
1747 
1748 
1749 
1750  }
1751 
1752 
1753 #ifdef MATVEC_PROFILE
1754  dendro::timer::sfcmatvec::t_accum.start();
1755 #endif
1756 
1757  for(unsigned int p=0;p<bCounts[child];p++)
1758  {
1759  pt_status=fem::seq::computeIJK(pt_internal[child][p],childElem[child],dxb2,logDxb2,ijk);
1760  if(pt_status==PT_INTERNAL_DG_NODE && !accumStatus[IDX(ijk[0],ijk[1],ijk[2])])
1761  v_internal[child][p]=v_unzip[IDX(ijk[0],ijk[1],ijk[2])];
1762 
1763  }
1764 
1765 
1766 #ifdef MATVEC_PROFILE
1767  dendro::timer::sfcmatvec::t_accum.stop();
1768 #endif
1769 
1770 
1771 #ifdef MATVEC_PROFILE
1772  dendro::timer::sfcmatvec::t_malloc.start();
1773 #endif
1774 
1775  delete [] accumStatus;
1776  delete [] c2pInterp;
1777 
1778  delete [] unzipStatus;
1779  delete [] interp3D_out;
1780 
1781  delete [] interp2D_in;
1782  delete [] interp2D_out;
1783 
1784  delete [] interp1D_in;
1785  delete [] interp1D_out;
1786 
1787  delete [] u_unzip;
1788  delete [] v_unzip;
1789 
1790 #ifdef MATVEC_PROFILE
1791  dendro::timer::sfcmatvec::t_malloc.stop();
1792 #endif
1793 
1794 
1795 
1796 
1797  }
1798 
1799 /*
1800 #ifdef MATVEC_PROFILE
1801  dendro::timer::sfcmatvec::t_accum.start();
1802 #endif
1803 
1804  for(unsigned int p=pBegin;p<pEnd;p++)
1805  {
1806  pt_status=fem::seq::computeIJK(pt_coords[p],childElem[child],dxb2,logDxb2,ijk);
1807  if(pt_status!=PT_EXTERNAL)
1808  {
1809  out[p]+=v_internal[child][counts[child]];
1810  counts[child]++;
1811  }
1812 
1813  }
1814 
1815 
1816 #ifdef MATVEC_PROFILE
1817  dendro::timer::sfcmatvec::t_accum.stop();
1818 #endif
1819 
1820 
1821 
1822 
1823 
1824 
1825 
1826  delete [] u_internal[child];
1827  delete [] v_internal[child];
1828  delete [] pt_internal[child];
1829 */
1830 
1831 
1832 
1833  }
1834 
1835 
1836 #ifdef MATVEC_PROFILE
1837  dendro::timer::sfcmatvec::t_pts_p1_accum.start();
1838 #endif
1839  // pass 1
1840  for(unsigned int p=pBegin;p<pEnd;p++)
1841  {
1842 
1843  ijk[2]=0;
1844  ijk[1]=0;
1845  ijk[0]=0;
1846 
1847  if(pt_coords[p].zint()>pMid[2]) ijk[2]=2;
1848  if(pt_coords[p].yint()>pMid[1]) ijk[1]=2;
1849  if(pt_coords[p].xint()>pMid[0]) ijk[0]=2;
1850 
1851  if(pt_coords[p].zint()<pMid[2]) ijk[2]=0;
1852  if(pt_coords[p].yint()<pMid[1]) ijk[1]=0;
1853  if(pt_coords[p].xint()<pMid[0]) ijk[0]=0;
1854 
1855  if(pt_coords[p].zint()==pMid[2]) ijk[2]=1;
1856  if(pt_coords[p].yint()==pMid[1]) ijk[1]=1;
1857  if(pt_coords[p].xint()==pMid[0]) ijk[0]=1;
1858 
1859 
1860  if(ijk[0]!=1 && ijk[1]!=1 && ijk[2]!=1) {
1861 
1862  cnum = (((ijk[2] >> 1u) << 2u) | ((ijk[1] >> 1u) << 1u) | ((ijk[0] >> 1u)));
1863  out[p]+=v_internal[cnum][counts[cnum]];
1864  counts[cnum]++;
1865 
1866  }
1867 
1868  }
1869 #ifdef MATVEC_PROFILE
1870  dendro::timer::sfcmatvec::t_pts_p1_accum.stop();
1871 #endif
1872 
1873 #ifdef MATVEC_PROFILE
1874  dendro::timer::sfcmatvec::t_pts_p2_accum.start();
1875 #endif
1876  // pass 2;
1877  for(unsigned int index=0;index<duplicateCoordIndex.size();index++) {
1878 
1879  ijk[2] = 0;
1880  ijk[1] = 0;
1881  ijk[0] = 0;
1882 
1883  if (pt_coords[duplicateCoordIndex[index]].zint() > pMid[2]) ijk[2] = 2;
1884  if (pt_coords[duplicateCoordIndex[index]].yint() > pMid[1]) ijk[1] = 2;
1885  if (pt_coords[duplicateCoordIndex[index]].xint() > pMid[0]) ijk[0] = 2;
1886 
1887  if (pt_coords[duplicateCoordIndex[index]].zint() < pMid[2]) ijk[2] = 0;
1888  if (pt_coords[duplicateCoordIndex[index]].yint() < pMid[1]) ijk[1] = 0;
1889  if (pt_coords[duplicateCoordIndex[index]].xint() < pMid[0]) ijk[0] = 0;
1890 
1891  if (pt_coords[duplicateCoordIndex[index]].zint() == pMid[2]) ijk[2] = 1;
1892  if (pt_coords[duplicateCoordIndex[index]].yint() == pMid[1]) ijk[1] = 1;
1893  if (pt_coords[duplicateCoordIndex[index]].xint() == pMid[0]) ijk[0] = 1;
1894 
1895  for(unsigned int child=0;child<NUM_CHILDREN;child++)
1896  {
1897  cnum=SFC_MATVEC_BUCKET_TABLE[(ijk[2]*SFC_MVEC_TABLE_Rz+ijk[1]*SFC_MVEC_TABLE_Ry+ijk[0])][child];
1898  if(cnum==UCHAR_MAX)
1899  break;
1900 
1901  out[duplicateCoordIndex[index]]+=v_internal[cnum][counts[cnum]];
1902  counts[cnum]++;
1903 
1904 
1905  }
1906 
1907  }
1908 
1909 #ifdef MATVEC_PROFILE
1910  dendro::timer::sfcmatvec::t_pts_p2_accum.stop();
1911 #endif
1912 
1913 
1914 
1915 #ifdef MATVEC_PROFILE
1916  dendro::timer::sfcmatvec::t_malloc.start();
1917 #endif
1918 
1919 
1920  for(unsigned int child=0;child<NUM_CHILDREN;child++)
1921  {
1922  delete [] u_internal[child];
1923  delete [] v_internal[child];
1924  delete [] pt_internal[child];
1925  }
1926 
1927 
1928  delete [] parentVal;
1929  delete [] u_internal;
1930  delete [] v_internal;
1931  delete [] pt_internal;
1932 
1933 #ifdef MATVEC_PROFILE
1934  dendro::timer::sfcmatvec::t_malloc.stop();
1935 #endif
1936 
1937 
1938 
1939 }
1940 
1941 
1942 #endif //SFCSORTBENCH_FEMMATVEC_H
RefElement::I2D_Parent2Child
void I2D_Parent2Child(const double *in, double *out, unsigned int childNum) const
Definition: refel.h:460
ot::TreeNode::minX
unsigned int minX() const
returns true min corner(anchor) and the max corner coordinates of a octant.
Definition: TreeNode.h:164
binOp::getNextHighestPowerOfTwo
int getNextHighestPowerOfTwo(unsigned int n)
Definition: binUtils.cpp:61
ot::TreeNode
A class to manage octants.
Definition: TreeNode.h:35
RefElement::I1D_Parent2Child
void I1D_Parent2Child(const double *in, double *out, unsigned int childNUm) const
Definition: refel.h:555
fem
workspace variables allocations
Definition: matvec.h:44
binOp::fastLog2
unsigned int fastLog2(unsigned int num)
Definition: binUtils.cpp:15
RefElement::I1D_Child2Parent
void I1D_Child2Parent(const double *in, double *out, unsigned int childNUm) const
Definition: refel.h:596
RefElement::I2D_Child2Parent
void I2D_Child2Parent(const double *in, double *out, unsigned int childNum) const
Definition: refel.h:505
ot::TreeNode::getLevel
unsigned int getLevel() const
return the level of the of the octant
Definition: TreeNode.h:387
seq
Collection of Generic Sequential Functions.
binUtils.h
A set of efficient functions that use binary operations to perform some small computations.
binOp::getBit
unsigned int getBit(T val, unsigned int i)
gets the i^th bit on the value val
Definition: binUtils.h:70
RefElement
Definition: refel.h:69
Generated by   doxygen 1.8.13