Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessEquiSpacedOutput.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessEquiSpacedOutput.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Set up fields as interpolation to equispaced output
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 #include <string>
36 #include <iostream>
37 using namespace std;
38 
40 
44 #include <boost/math/special_functions/fpclassify.hpp>
45 #include <StdRegions/StdQuadExp.h>
46 #include <StdRegions/StdTriExp.h>
47 
48 namespace Nektar
49 {
50 namespace Utilities
51 {
52 
53 ModuleKey ProcessEquiSpacedOutput::className =
55  ModuleKey(eProcessModule, "equispacedoutput"),
56  ProcessEquiSpacedOutput::create,
57  "Write data as equi-spaced output using simplices to represent the data for connecting points");
58 
59 
60 ProcessEquiSpacedOutput::ProcessEquiSpacedOutput(FieldSharedPtr f)
61  : ProcessModule(f)
62 {
63  f->m_setUpEquiSpacedFields = true;
64 
65  m_config["tetonly"] = ConfigOption(true, "NotSet",
66  "Only process tetrahedral elements");
67 
68  m_config["modalenergy"] = ConfigOption(true,"NotSet","Write output as modal energy");
69 
70 }
71 
73 {
74 }
75 
76 void ProcessEquiSpacedOutput::Process(po::variables_map &vm)
77 {
79 }
80 
81 
83 {
84  if(m_f->m_verbose)
85  {
86  if(m_f->m_comm->GetRank() == 0)
87  {
88  cout << "Interpolating fields to equispaced..." << endl;
89  }
90  }
91 
92  int coordim = m_f->m_exp[0]->GetCoordim(0);
93  int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
94  int npts = m_f->m_exp[0]->GetTotPoints();
96 
97  int nel = m_f->m_exp[0]->GetExpSize();
98 
99  // Check if we have a homogeneous expansion
100  bool homogeneous1D = false;
101  if (m_f->m_fielddef.size())
102  {
103  if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
104  {
105  ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
106  coordim++;
107  shapedim++;
108  homogeneous1D = true;
109  }
110  else if(m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
111  {
112  ASSERTL0(false, "Homegeneous2D case not supported");
113  }
114  }
115  else
116  {
117  if(m_f->m_session->DefinesSolverInfo("HOMOGENEOUS"))
118  {
119  std::string HomoStr = m_f->m_session->GetSolverInfo("HOMOGENEOUS");
120 
121  if((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D")
122  || (HomoStr == "1D") || (HomoStr == "Homo1D"))
123  {
124  ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
125  coordim++;
126  shapedim++;
127  homogeneous1D = true;
128  }
129  if((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D")
130  || (HomoStr == "2D") || (HomoStr == "Homo2D"))
131  {
132  ASSERTL0(false, "Homegeneous2D case not supported");
133  }
134  }
135  }
136 
137  // set up the number of points in each element
138  int newpoints = 0;
139  int newtotpoints = 0;
140 
141  Array<OneD,int> conn;
142  int prevNcoeffs = 0;
143  int prevNpoints = 0;
144  int cnt = 0;
145 
146  // identify face 1 connectivity for prisms
147  map<int,StdRegions::Orientation > face0orient;
148  set<int> prismorient;
150 
151  // prepare PtsField
152  vector<std::string> fieldNames;
153  vector<int> ppe;
154  vector<Array<OneD, int> > ptsConn;
155  int nfields;
156 
157  for(int i = 0; i < nel; ++i)
158  {
159  e = m_f->m_exp[0]->GetExp(i);
160  if(e->DetShapeType() == LibUtilities::ePrism)
161  {
162  StdRegions::Orientation forient = e->GetForient(0);
163  int fid = e->GetGeom()->GetFid(0);
164  if(face0orient.count(fid))
165  { // face 1 meeting face 1 so reverse this id
166  prismorient.insert(i);
167  }
168  else
169  {
170  // just store if Dir 1 is fwd or bwd
171  if((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
175  {
176  face0orient[fid] = StdRegions::eBwd;
177  }
178  else
179  {
180  face0orient[fid] = StdRegions::eFwd;
181  }
182  }
183  }
184  }
185 
186  for(int i = 0; i < nel; ++i)
187  {
188  e = m_f->m_exp[0]->GetExp(i);
189  if(e->DetShapeType() == LibUtilities::ePrism)
190  {
191  int fid = e->GetGeom()->GetFid(2);
192  // check to see if face 2 meets face 1
193  if(face0orient.count(fid))
194  {
195  // check to see how face 2 is orientated
196  StdRegions::Orientation forient2 = e->GetForient(2);
197  StdRegions::Orientation forient0 = face0orient[fid];
198 
199  // If dir 1 or forient2 is bwd then check agains
200  // face 1 value
201  if((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
202  (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
203  (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
205  {
206  if(forient0 == StdRegions::eFwd)
207  {
208  prismorient.insert(i);
209  }
210  }
211  else
212  {
213  if(forient0 == StdRegions::eBwd)
214  {
215  prismorient.insert(i);
216  }
217  }
218  }
219  }
220  }
221 
222  for(int i = 0; i < nel; ++i)
223  {
224  e = m_f->m_exp[0]->GetExp(i);
225  if(m_config["tetonly"].m_beenSet)
226  {
227  if(m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
229  {
230  continue;
231  }
232  }
233 
234  switch(e->DetShapeType())
235  {
237  {
238  int npoints0 = e->GetBasis(0)->GetNumPoints();
239 
240  newpoints = LibUtilities::StdSegData::
241  getNumberOfCoefficients(npoints0);
242  }
243  break;
245  {
246  int np0 = e->GetBasis(0)->GetNumPoints();
247  int np1 = e->GetBasis(1)->GetNumPoints();
248  int np = max(np0,np1);
249  newpoints = LibUtilities::StdTriData::
251  }
252  break;
254  {
255  int np0 = e->GetBasis(0)->GetNumPoints();
256  int np1 = e->GetBasis(1)->GetNumPoints();
257  int np = max(np0,np1);
258 
259  newpoints = LibUtilities::StdQuadData::
261  }
262  break;
264  {
265  int np0 = e->GetBasis(0)->GetNumPoints();
266  int np1 = e->GetBasis(1)->GetNumPoints();
267  int np2 = e->GetBasis(2)->GetNumPoints();
268  int np = max(np0,max(np1,np2));
269 
270  newpoints = LibUtilities::StdTetData::
271  getNumberOfCoefficients(np,np,np);
272  }
273  break;
275  {
276  int np0 = e->GetBasis(0)->GetNumPoints();
277  int np1 = e->GetBasis(1)->GetNumPoints();
278  int np2 = e->GetBasis(2)->GetNumPoints();
279  int np = max(np0,max(np1,np2));
280 
281  newpoints = LibUtilities::StdPrismData::
282  getNumberOfCoefficients(np,np,np);
283  }
284  break;
286  {
287  int np0 = e->GetBasis(0)->GetNumPoints();
288  int np1 = e->GetBasis(1)->GetNumPoints();
289  int np2 = e->GetBasis(2)->GetNumPoints();
290  int np = max(np0,max(np1,np2));
291 
292  newpoints = LibUtilities::StdPyrData::
293  getNumberOfCoefficients(np,np,np);
294  }
295  break;
297  {
298  int np0 = e->GetBasis(0)->GetNumPoints();
299  int np1 = e->GetBasis(1)->GetNumPoints();
300  int np2 = e->GetBasis(2)->GetNumPoints();
301  int np = max(np0,max(np1,np2));
302 
303  newpoints = LibUtilities::StdHexData::
304  getNumberOfCoefficients(np,np,np);
305  }
306  break;
307  default:
308  {
309  ASSERTL0(false,"Points not known");
310  }
311  }
312 
313  ppe.push_back(newpoints);
314  newtotpoints += newpoints;
315 
316 
317  if(e->DetShapeType() == LibUtilities::ePrism)
318  {
319  bool standard = true;
320 
321  if(prismorient.count(i))
322  {
323  standard = false; // reverse direction
324  }
325 
326  e->GetSimplexEquiSpacedConnectivity(conn,standard);
327  }
328  else
329  {
330 
331  if((prevNcoeffs != e->GetNcoeffs()) ||
332  (prevNpoints != e->GetTotPoints()))
333  {
334  prevNcoeffs = e->GetNcoeffs();
335  prevNpoints = e->GetTotPoints();
336 
337  e->GetSimplexEquiSpacedConnectivity(conn);
338  }
339  }
340  Array<OneD, int> newconn(conn.num_elements());
341  for(int j = 0; j < conn.num_elements(); ++j)
342  {
343  newconn[j] = conn[j] + cnt;
344  }
345 
346  ptsConn.push_back(newconn);
347  cnt += newpoints;
348  }
349 
350  if(m_f->m_fielddef.size())
351  {
352  nfields = m_f->m_exp.size();
353  }
354  else // just the mesh points
355  {
356  nfields = 0;
357  }
358 
359  Array<OneD, Array<OneD, NekDouble> > pts(nfields + coordim);
360 
361  for(int i = 0; i < nfields + coordim; ++i)
362  {
363  pts[i] = Array<OneD, NekDouble>(newtotpoints);
364  }
365 
366  // Interpolate coordinates
367  for(int i = 0; i < coordim; ++i)
368  {
369  coords[i] = Array<OneD, NekDouble>(npts);
370  }
371 
372  for(int i = coordim; i < 3; ++i)
373  {
374  coords[i] = NullNekDouble1DArray;
375  }
376 
377  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
378 
379  int nq1 = m_f->m_exp[0]->GetTotPoints();
380 
381  Array<OneD, NekDouble> x1(nq1);
382  Array<OneD, NekDouble> y1(nq1);
383  Array<OneD, NekDouble> z1(nq1);
384 
385  m_f->m_exp[0]->GetCoords(x1, y1, z1);
386 
387 
389 
390  for(int n = 0; n < coordim; ++n)
391  {
392  cnt = 0;
393  int cnt1 = 0;
394  for(int i = 0; i < nel; ++i)
395  {
396  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
397  coords[n] + cnt,
398  tmp = pts[n] + cnt1);
399  cnt1 += ppe[i];
400  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
401  }
402  }
403 
404  if(m_f->m_fielddef.size())
405  {
406  ASSERTL0(m_f->m_fielddef[0]->m_fields.size() == m_f->m_exp.size(),
407  "More expansion defined than fields");
408 
409  for(int n = 0; n < m_f->m_exp.size(); ++n)
410  {
411  cnt = 0;
412  int cnt1 = 0;
413 
414  if(m_config["modalenergy"].m_beenSet)
415  {
416  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
417  for(int i = 0; i < nel; ++i)
418  {
419  GenOrthoModes(i,phys+cnt,tmp = pts[coordim + n] + cnt1);
420  cnt1 += ppe[i];
421  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
422  }
423  }
424  else
425  {
426  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
427  for(int i = 0; i < nel; ++i)
428  {
429  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
430  phys + cnt,
431  tmp = pts[coordim + n] + cnt1);
432  cnt1 += ppe[i];
433  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
434  }
435  }
436 
437  // Set up Variable string.
438  fieldNames.push_back(m_f->m_fielddef[0]->m_fields[n]);
439  }
440  }
441 
442  m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(coordim, fieldNames, pts);
443  if (shapedim == 2)
444  {
445  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
446  }
447  else if (shapedim == 3)
448  {
449  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
450  }
451  m_f->m_fieldPts->SetConnectivity(ptsConn);
452  if (homogeneous1D)
453  {
455  }
456 }
457 
459 {
461  int nel = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
462  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
463  int npts = m_f->m_fieldPts->GetNpoints();
464  int nptsPerPlane = npts/nPlanes;
465 
466  if ( m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
468  && m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
470  {
471  // Write points with extra plane
473  m_f->m_fieldPts->GetPts(pts);
474  Array<OneD, Array<OneD, NekDouble> > newPts(pts.num_elements());
475  for (int i = 0; i < pts.num_elements(); i++)
476  {
477  newPts[i] = Array<OneD, NekDouble> (npts + nptsPerPlane);
478  // Copy old points
479  Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
480  // Get new plane
481  Array<OneD, NekDouble> extraPlane (nptsPerPlane, newPts[i]+npts);
482  if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
483  {
484  if ( i == 2)
485  {
486  // z-coordinate
487  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
488  pts[i], 1, extraPlane, 1);
489  }
490  else
491  {
492  Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
493  }
494  }
495  else
496  {
497  // Determine to and from rank for communication
498  int size = m_f->m_comm->GetColumnComm()->GetSize();
499  int rank = m_f->m_comm->GetColumnComm()->GetRank();
500  int fromRank = (rank+1) % size;
501  int toRank = (rank == 0) ? size-1 : rank-1;
502  // Communicate using SendRecv
503  Array<OneD, NekDouble> send (nptsPerPlane, pts[i]);
504  m_f->m_comm->GetColumnComm()->SendRecv(toRank, send,
505  fromRank, extraPlane);
506  // Adjust z-coordinate of last process
507  if (i == 2 && (rank == size - 1) )
508  {
509  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
510  extraPlane, 1, extraPlane, 1);
511  }
512  }
513  }
514  m_f->m_fieldPts->SetPts(newPts);
515  // Now extend plane connectivity
516  vector<Array<OneD, int> > oldConn;
517  Array<OneD, int> conn;
518  m_f->m_fieldPts->GetConnectivity(oldConn);
519  vector<Array<OneD, int> > newConn = oldConn;
520  int connPerPlane = oldConn.size()/nPlanes;
521  for (int i = 0; i < connPerPlane; i++)
522  {
523  conn = Array<OneD, int> (oldConn[i].num_elements());
524  for (int j = 0; j < conn.num_elements(); j++)
525  {
526  conn[j] = oldConn[i][j] + npts;
527  }
528  newConn.push_back(conn);
529  }
530  m_f->m_fieldPts->SetConnectivity(newConn);
531 
532  nPlanes++;
533  npts += nptsPerPlane;
534  }
535 
536  // Get maximum number of points per direction in the mesh
537  int maxN = 0;
538  for(int i = 0; i < nel; ++i)
539  {
540  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
541 
542  int np0 = e->GetBasis(0)->GetNumPoints();
543  int np1 = e->GetBasis(1)->GetNumPoints();
544  int np = max(np0,np1);
545  maxN = max(np, maxN);
546  }
547 
548  // Create a global numbering for points in plane 0
549  Array<OneD, int> vId(nptsPerPlane);
550  int cnt1=0;
551  int cnt2=0;
552  for(int i = 0; i < nel; ++i)
553  {
554  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
555  int np0 = e->GetBasis(0)->GetNumPoints();
556  int np1 = e->GetBasis(1)->GetNumPoints();
557  int np = max(np0,np1);
558  switch(e->DetShapeType())
559  {
561  {
562  int newpoints = LibUtilities::StdTriData::
564 
565  // Vertex numbering
566  vId[cnt1] = e->GetGeom()->GetVid(0);
567  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
568  vId[cnt1+newpoints-1] = e->GetGeom()->GetVid(2);
569 
570  // Edge numbering
571  StdRegions::Orientation edgeOrient;
572  int edge1 = 0;
573  int edge2 = 0;
574  for (int n = 1; n < np-1; n++)
575  {
576  // Edge 0
577  edgeOrient = e->GetGeom()->GetEorient(0);
578  if (edgeOrient==StdRegions::eForwards)
579  {
580  vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
581  }
582  else
583  {
584  vId[cnt1+np-1-n] = 4*nel +
585  maxN*e->GetGeom()->GetEid(0) + n;
586  }
587 
588  // Edge 1
589  edgeOrient = e->GetGeom()->GetEorient(1);
590  if (edgeOrient==StdRegions::eForwards)
591  {
592  edge1 += np-n;
593  vId[cnt1+np-1+edge1] = 4*nel +
594  maxN*e->GetGeom()->GetEid(1) + n;
595  }
596  else
597  {
598  edge1 += n;
599  vId[cnt1+newpoints-1-edge1] = 4*nel +
600  maxN*e->GetGeom()->GetEid(1) + n;
601  }
602 
603  // Edge 2
604  edgeOrient = e->GetGeom()->GetEorient(2);
605  if (edgeOrient==StdRegions::eForwards)
606  {
607  edge2 += n+1;
608  vId[cnt1+newpoints-1-edge2] = 4*nel +
609  maxN*e->GetGeom()->GetEid(2) + n;
610  }
611  else
612  {
613  edge2 += np+1-n;
614  vId[cnt1+edge2] = 4*nel +
615  maxN*e->GetGeom()->GetEid(2) + n;
616  }
617  }
618 
619  // Interior numbering
620  int mode = np+1;
621  for (int n = 1; n < np-1; n++)
622  {
623  for (int m = 1; m < np-n-1; m++)
624  {
625  vId[cnt1+mode] = 4*nel + maxN*4*nel + cnt2;
626  cnt2++;
627  mode++;
628  }
629  }
630  cnt1+= newpoints;
631  }
632  break;
634  {
635  int newpoints = LibUtilities::StdQuadData::
637  // Vertex numbering
638  vId[cnt1] = e->GetGeom()->GetVid(0);
639  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
640  vId[cnt1+np*np-1] = e->GetGeom()->GetVid(2);
641  vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);
642 
643  // Edge numbering
644  StdRegions::Orientation edgeOrient;
645  for (int n = 1; n < np-1; n++)
646  {
647  // Edge 0
648  edgeOrient = e->GetGeom()->GetEorient(0);
649  if (edgeOrient==StdRegions::eForwards)
650  {
651  vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
652  }
653  else
654  {
655  vId[cnt1+np-1-n] = 4*nel +
656  maxN*e->GetGeom()->GetEid(0) + n;
657  }
658 
659  // Edge 1
660  edgeOrient = e->GetGeom()->GetEorient(1);
661  if (edgeOrient==StdRegions::eForwards)
662  {
663  vId[cnt1+np-1+n*np] = 4*nel +
664  maxN*e->GetGeom()->GetEid(1) + n;
665  }
666  else
667  {
668  vId[cnt1+np*np-1-n*np] = 4*nel +
669  maxN*e->GetGeom()->GetEid(1) + n;
670  }
671 
672  // Edge 2
673  edgeOrient = e->GetGeom()->GetEorient(2);
674  if (edgeOrient==StdRegions::eForwards)
675  {
676  vId[cnt1+np*np-1-n] = 4*nel +
677  maxN*e->GetGeom()->GetEid(2) + n;
678  }
679  else
680  {
681  vId[cnt1+np*(np-1)+n] = 4*nel +
682  maxN*e->GetGeom()->GetEid(2) + n;
683  }
684 
685  // Edge 3
686  edgeOrient = e->GetGeom()->GetEorient(3);
687  if (edgeOrient==StdRegions::eForwards)
688  {
689  vId[cnt1+np*(np-1)-n*np] = 4*nel +
690  maxN*e->GetGeom()->GetEid(3) + n;
691  }
692  else
693  {
694  vId[cnt1+n*np] = 4*nel +
695  maxN*e->GetGeom()->GetEid(3) + n;
696  }
697  }
698 
699  // Interior numbering
700  for (int n = 1; n < np-1; n++)
701  {
702  for (int m = 1; m < np-1; m++)
703  {
704  vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
705  cnt2++;
706  }
707  }
708  cnt1+= newpoints;
709  }
710  break;
711  default:
712  {
713  ASSERTL0(false,"Points not known");
714  }
715  }
716  }
717 
718  // Crete new connectivity using homogeneous information
719  vector<Array<OneD, int> > oldConn;
720  vector<Array<OneD, int> > newConn;
721  Array<OneD, int> conn;
722  m_f->m_fieldPts->GetConnectivity(oldConn);
723 
724  for(int i = 0; i < nel; ++i)
725  {
726  int nTris = oldConn[i].num_elements()/3;
727  // Create array for new connectivity
728  // (3 tetrahedra between each plane for each triangle)
729  conn = Array<OneD, int> (4*3*nTris*(nPlanes-1));
730  cnt2=0;
731  for(int n = 0; n<nTris; n++)
732  {
733  // Get id of vertices in this triangle (on plane 0)
734  int vId0 = vId[oldConn[i][3*n+0]];
735  int vId1 = vId[oldConn[i][3*n+1]];
736  int vId2 = vId[oldConn[i][3*n+2]];
737 
738  // Determine ordering based on global Id of pts
739  Array<OneD, int> rot(3);
740  if ( (vId0<vId1) && (vId0<vId2))
741  {
742  rot[0] = 0;
743  if (vId1<vId2)
744  {
745  rot[1] = 1;
746  rot[2] = 2;
747  }
748  else
749  {
750  rot[1] = 2;
751  rot[2] = 1;
752  }
753  }
754  else if ( (vId1<vId0) && (vId1<vId2))
755  {
756  rot[0] = 1;
757  if (vId0<vId2)
758  {
759  rot[1] = 0;
760  rot[2] = 2;
761  }
762  else
763  {
764  rot[1] = 2;
765  rot[2] = 0;
766  }
767  }
768  else
769  {
770  rot[0] = 2;
771  if (vId0<vId1)
772  {
773  rot[1] = 0;
774  rot[2] = 1;
775  }
776  else
777  {
778  rot[1] = 1;
779  rot[2] = 0;
780  }
781  }
782 
783  // Define new connectivity with tetrahedra
784  for ( int p = 0; p<nPlanes-1; p++)
785  {
786  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
787  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
788  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[2]];
789  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
790 
791  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
792  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
793  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
794  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
795 
796  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
797  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
798  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
799  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
800  }
801  }
802  newConn.push_back(conn);
803  }
804  m_f->m_fieldPts->SetConnectivity(newConn);
805 }
806 
807 
809  int n,
810  const Array<OneD,const NekDouble> &phys,
811  Array<OneD, NekDouble> &coeffs)
812  {
814  e = m_f->m_exp[0]->GetExp(n);
815 
816  switch(e->DetShapeType())
817  {
819  {
820  int np0 = e->GetBasis(0)->GetNumPoints();
821  int np1 = e->GetBasis(1)->GetNumPoints();
822  int np = max(np0,np1);
823 
824  // to ensure points are correctly projected to np need
825  // to increase the order slightly of coordinates
826  LibUtilities::PointsKey pa(np+1,e->GetPointsType(0));
827  LibUtilities::PointsKey pb(np,e->GetPointsType(1));
828  Array<OneD, NekDouble> tophys(np*(np+1));
829 
832  StdRegions::StdTriExp OrthoExp(Ba,Bb);
833 
834  // interpolate points to new phys points!
835  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
836  e->GetBasis(1)->GetBasisKey(),
837  phys,Ba,Bb,tophys);
838 
839  OrthoExp.FwdTrans(tophys,coeffs);
840  break;
841  }
843  {
844  int np0 = e->GetBasis(0)->GetNumPoints();
845  int np1 = e->GetBasis(1)->GetNumPoints();
846  int np = max(np0,np1);
847 
848  LibUtilities::PointsKey pa(np+1,e->GetPointsType(0));
849  LibUtilities::PointsKey pb(np+1,e->GetPointsType(1));
850  Array<OneD, NekDouble> tophys((np+1)*(np+1));
851 
854  StdRegions::StdQuadExp OrthoExp(Ba,Bb);
855 
856  // interpolate points to new phys points!
857  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
858  e->GetBasis(1)->GetBasisKey(),
859  phys,Ba,Bb,tophys);
860 
861  OrthoExp.FwdTrans(phys,coeffs);
862  break;
863  }
864  default:
865  ASSERTL0(false,"Shape needs setting up");
866  break;
867  }
868  }
869 
870 }
871 }
872 
873 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:286
pair< ModuleType, string > ModuleKey
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
void GenOrthoModes(int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
Fourier Expansion .
Definition: BasisType.h:52
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:186
FieldSharedPtr m_f
Field object.
Principle Orthogonal Functions .
Definition: BasisType.h:47
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
static std::string npts
Definition: InputFld.cpp:43
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Represents a command-line configuration option.
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:151
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Describes the specification for a Basis.
Definition: Basis.h:50
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...