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 
85  if(m_f->m_verbose)
86  {
87  cout << "Interpolating fields to equispaced" << endl;
88  }
89 
90  int coordim = m_f->m_exp[0]->GetCoordim(0);
91  int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
92  int npts = m_f->m_exp[0]->GetTotPoints();
94 
95  int nel = m_f->m_exp[0]->GetExpSize();
96 
97  // Check if we have a homogeneous expansion
98  bool homogeneous1D = false;
99  if (m_f->m_fielddef.size())
100  {
101  if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
102  {
103  ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
104  coordim++;
105  shapedim++;
106  homogeneous1D = true;
107  }
108  else if(m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
109  {
110  ASSERTL0(false, "Homegeneous2D case not supported");
111  }
112  }
113  else
114  {
115  if(m_f->m_session->DefinesSolverInfo("HOMOGENEOUS"))
116  {
117  std::string HomoStr = m_f->m_session->GetSolverInfo("HOMOGENEOUS");
118 
119  if((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D")
120  || (HomoStr == "1D") || (HomoStr == "Homo1D"))
121  {
122  ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
123  coordim++;
124  shapedim++;
125  homogeneous1D = true;
126  }
127  if((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D")
128  || (HomoStr == "2D") || (HomoStr == "Homo2D"))
129  {
130  ASSERTL0(false, "Homegeneous2D case not supported");
131  }
132  }
133  }
134 
135  // set up the number of points in each element
136  int newpoints = 0;
137  int newtotpoints = 0;
138 
139  Array<OneD,int> conn;
140  int prevNcoeffs = 0;
141  int prevNpoints = 0;
142  int cnt = 0;
143 
144  // identify face 1 connectivity for prisms
145  map<int,StdRegions::Orientation > face0orient;
146  set<int> prismorient;
148 
149  // prepare PtsField
150  vector<std::string> fieldNames;
151  vector<int> ppe;
152  vector<Array<OneD, int> > ptsConn;
153  int nfields;
154 
155  for(int i = 0; i < nel; ++i)
156  {
157  e = m_f->m_exp[0]->GetExp(i);
158  if(e->DetShapeType() == LibUtilities::ePrism)
159  {
160  StdRegions::Orientation forient = e->GetForient(0);
161  int fid = e->GetGeom()->GetFid(0);
162  if(face0orient.count(fid))
163  { // face 1 meeting face 1 so reverse this id
164  prismorient.insert(i);
165  }
166  else
167  {
168  // just store if Dir 1 is fwd or bwd
169  if((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
173  {
174  face0orient[fid] = StdRegions::eBwd;
175  }
176  else
177  {
178  face0orient[fid] = StdRegions::eFwd;
179  }
180  }
181  }
182  }
183 
184  for(int i = 0; i < nel; ++i)
185  {
186  e = m_f->m_exp[0]->GetExp(i);
187  if(e->DetShapeType() == LibUtilities::ePrism)
188  {
189  int fid = e->GetGeom()->GetFid(2);
190  // check to see if face 2 meets face 1
191  if(face0orient.count(fid))
192  {
193  // check to see how face 2 is orientated
194  StdRegions::Orientation forient2 = e->GetForient(2);
195  StdRegions::Orientation forient0 = face0orient[fid];
196 
197  // If dir 1 or forient2 is bwd then check agains
198  // face 1 value
199  if((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
200  (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
201  (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
203  {
204  if(forient0 == StdRegions::eFwd)
205  {
206  prismorient.insert(i);
207  }
208  }
209  else
210  {
211  if(forient0 == StdRegions::eBwd)
212  {
213  prismorient.insert(i);
214  }
215  }
216  }
217  }
218  }
219 
220  for(int i = 0; i < nel; ++i)
221  {
222  e = m_f->m_exp[0]->GetExp(i);
223  if(m_config["tetonly"].m_beenSet)
224  {
225  if(m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
227  {
228  continue;
229  }
230  }
231 
232  switch(e->DetShapeType())
233  {
235  {
236  int npoints0 = e->GetBasis(0)->GetNumPoints();
237 
238  newpoints = LibUtilities::StdSegData::
239  getNumberOfCoefficients(npoints0);
240  }
241  break;
243  {
244  int np0 = e->GetBasis(0)->GetNumPoints();
245  int np1 = e->GetBasis(1)->GetNumPoints();
246  int np = max(np0,np1);
247  newpoints = LibUtilities::StdTriData::
249  }
250  break;
252  {
253  int np0 = e->GetBasis(0)->GetNumPoints();
254  int np1 = e->GetBasis(1)->GetNumPoints();
255  int np = max(np0,np1);
256 
257  newpoints = LibUtilities::StdQuadData::
259  }
260  break;
262  {
263  int np0 = e->GetBasis(0)->GetNumPoints();
264  int np1 = e->GetBasis(1)->GetNumPoints();
265  int np2 = e->GetBasis(2)->GetNumPoints();
266  int np = max(np0,max(np1,np2));
267 
268  newpoints = LibUtilities::StdTetData::
269  getNumberOfCoefficients(np,np,np);
270  }
271  break;
273  {
274  int np0 = e->GetBasis(0)->GetNumPoints();
275  int np1 = e->GetBasis(1)->GetNumPoints();
276  int np2 = e->GetBasis(2)->GetNumPoints();
277  int np = max(np0,max(np1,np2));
278 
279  newpoints = LibUtilities::StdPrismData::
280  getNumberOfCoefficients(np,np,np);
281  }
282  break;
284  {
285  int np0 = e->GetBasis(0)->GetNumPoints();
286  int np1 = e->GetBasis(1)->GetNumPoints();
287  int np2 = e->GetBasis(2)->GetNumPoints();
288  int np = max(np0,max(np1,np2));
289 
290  newpoints = LibUtilities::StdPyrData::
291  getNumberOfCoefficients(np,np,np);
292  }
293  break;
295  {
296  int np0 = e->GetBasis(0)->GetNumPoints();
297  int np1 = e->GetBasis(1)->GetNumPoints();
298  int np2 = e->GetBasis(2)->GetNumPoints();
299  int np = max(np0,max(np1,np2));
300 
301  newpoints = LibUtilities::StdHexData::
302  getNumberOfCoefficients(np,np,np);
303  }
304  break;
305  default:
306  {
307  ASSERTL0(false,"Points not known");
308  }
309  }
310 
311  ppe.push_back(newpoints);
312  newtotpoints += newpoints;
313 
314 
315  if(e->DetShapeType() == LibUtilities::ePrism)
316  {
317  bool standard = true;
318 
319  if(prismorient.count(i))
320  {
321  standard = false; // reverse direction
322  }
323 
324  e->GetSimplexEquiSpacedConnectivity(conn,standard);
325  }
326  else
327  {
328 
329  if((prevNcoeffs != e->GetNcoeffs()) ||
330  (prevNpoints != e->GetTotPoints()))
331  {
332  prevNcoeffs = e->GetNcoeffs();
333  prevNpoints = e->GetTotPoints();
334 
335  e->GetSimplexEquiSpacedConnectivity(conn);
336  }
337  }
338  Array<OneD, int> newconn(conn.num_elements());
339  for(int j = 0; j < conn.num_elements(); ++j)
340  {
341  newconn[j] = conn[j] + cnt;
342  }
343 
344  ptsConn.push_back(newconn);
345  cnt += newpoints;
346  }
347 
348  if(m_f->m_fielddef.size())
349  {
350  nfields = m_f->m_exp.size();
351  }
352  else // just the mesh points
353  {
354  nfields = 0;
355  }
356 
357  Array<OneD, Array<OneD, NekDouble> > pts(nfields + coordim);
358 
359  for(int i = 0; i < nfields + coordim; ++i)
360  {
361  pts[i] = Array<OneD, NekDouble>(newtotpoints);
362  }
363 
364  // Interpolate coordinates
365  for(int i = 0; i < coordim; ++i)
366  {
367  coords[i] = Array<OneD, NekDouble>(npts);
368  }
369 
370  for(int i = coordim; i < 3; ++i)
371  {
372  coords[i] = NullNekDouble1DArray;
373  }
374 
375  m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
376 
377  int nq1 = m_f->m_exp[0]->GetTotPoints();
378 
379  Array<OneD, NekDouble> x1(nq1);
380  Array<OneD, NekDouble> y1(nq1);
381  Array<OneD, NekDouble> z1(nq1);
382 
383  m_f->m_exp[0]->GetCoords(x1, y1, z1);
384 
385 
387 
388  for(int n = 0; n < coordim; ++n)
389  {
390  cnt = 0;
391  int cnt1 = 0;
392  for(int i = 0; i < nel; ++i)
393  {
394  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
395  coords[n] + cnt,
396  tmp = pts[n] + cnt1);
397  cnt1 += ppe[i];
398  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
399  }
400  }
401 
402  if(m_f->m_fielddef.size())
403  {
404  ASSERTL0(m_f->m_fielddef[0]->m_fields.size() == m_f->m_exp.size(),
405  "More expansion defined than fields");
406 
407  for(int n = 0; n < m_f->m_exp.size(); ++n)
408  {
409  cnt = 0;
410  int cnt1 = 0;
411 
412  if(m_config["modalenergy"].m_beenSet)
413  {
414  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
415  for(int i = 0; i < nel; ++i)
416  {
417  GenOrthoModes(i,phys+cnt,tmp = pts[coordim + n] + cnt1);
418  cnt1 += ppe[i];
419  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
420  }
421  }
422  else
423  {
424  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
425  for(int i = 0; i < nel; ++i)
426  {
427  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
428  phys + cnt,
429  tmp = pts[coordim + n] + cnt1);
430  cnt1 += ppe[i];
431  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
432  }
433  }
434 
435  // Set up Variable string.
436  fieldNames.push_back(m_f->m_fielddef[0]->m_fields[n]);
437  }
438  }
439 
440  m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(coordim, fieldNames, pts);
441  if (shapedim == 2)
442  {
443  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
444  }
445  else if (shapedim == 3)
446  {
447  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
448  }
449  m_f->m_fieldPts->SetConnectivity(ptsConn);
450  if (homogeneous1D)
451  {
453  }
454 }
455 
457 {
459  int nel = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
460  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
461  int npts = m_f->m_fieldPts->GetNpoints();
462  int nptsPerPlane = npts/nPlanes;
463 
464  // Get maximum number of points per direction in the mesh
465  int maxN = 0;
466  for(int i = 0; i < nel; ++i)
467  {
468  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
469 
470  int np0 = e->GetBasis(0)->GetNumPoints();
471  int np1 = e->GetBasis(1)->GetNumPoints();
472  int np = max(np0,np1);
473  maxN = max(np, maxN);
474  }
475 
476  // Create a global numbering for points in plane 0
477  Array<OneD, int> vId(nptsPerPlane);
478  int cnt1=0;
479  int cnt2=0;
480  for(int i = 0; i < nel; ++i)
481  {
482  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
483  int np0 = e->GetBasis(0)->GetNumPoints();
484  int np1 = e->GetBasis(1)->GetNumPoints();
485  int np = max(np0,np1);
486  switch(e->DetShapeType())
487  {
489  {
490  int newpoints = LibUtilities::StdTriData::
492 
493  // Vertex numbering
494  vId[cnt1] = e->GetGeom()->GetVid(0);
495  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
496  vId[cnt1+newpoints-1] = e->GetGeom()->GetVid(2);
497 
498  // Edge numbering
499  StdRegions::Orientation edgeOrient;
500  int edge1 = 0;
501  int edge2 = 0;
502  for (int n = 1; n < np-1; n++)
503  {
504  // Edge 0
505  edgeOrient = e->GetGeom()->GetEorient(0);
506  if (edgeOrient==StdRegions::eForwards)
507  {
508  vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
509  }
510  else
511  {
512  vId[cnt1+np-1-n] = 4*nel +
513  maxN*e->GetGeom()->GetEid(0) + n;
514  }
515 
516  // Edge 1
517  edgeOrient = e->GetGeom()->GetEorient(1);
518  if (edgeOrient==StdRegions::eForwards)
519  {
520  edge1 += np-n;
521  vId[cnt1+np-1+edge1] = 4*nel +
522  maxN*e->GetGeom()->GetEid(1) + n;
523  }
524  else
525  {
526  edge1 += n;
527  vId[cnt1+newpoints-1-edge1] = 4*nel +
528  maxN*e->GetGeom()->GetEid(1) + n;
529  }
530 
531  // Edge 2
532  edgeOrient = e->GetGeom()->GetEorient(2);
533  if (edgeOrient==StdRegions::eForwards)
534  {
535  edge2 += n+1;
536  vId[cnt1+newpoints-1-edge2] = 4*nel +
537  maxN*e->GetGeom()->GetEid(2) + n;
538  }
539  else
540  {
541  edge2 += np+1-n;
542  vId[cnt1+edge2] = 4*nel +
543  maxN*e->GetGeom()->GetEid(2) + n;
544  }
545  }
546 
547  // Interior numbering
548  int mode = np+1;
549  for (int n = 1; n < np-1; n++)
550  {
551  for (int m = 1; m < np-n-1; m++)
552  {
553  vId[cnt1+mode] = 4*nel + maxN*4*nel + cnt2;
554  cnt2++;
555  mode++;
556  }
557  }
558  cnt1+= newpoints;
559  }
560  break;
562  {
563  int newpoints = LibUtilities::StdQuadData::
565  // Vertex numbering
566  vId[cnt1] = e->GetGeom()->GetVid(0);
567  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
568  vId[cnt1+np*np-1] = e->GetGeom()->GetVid(2);
569  vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);
570 
571  // Edge numbering
572  StdRegions::Orientation edgeOrient;
573  for (int n = 1; n < np-1; n++)
574  {
575  // Edge 0
576  edgeOrient = e->GetGeom()->GetEorient(0);
577  if (edgeOrient==StdRegions::eForwards)
578  {
579  vId[cnt1+n] = 4*nel + maxN*e->GetGeom()->GetEid(0) + n;
580  }
581  else
582  {
583  vId[cnt1+np-1-n] = 4*nel +
584  maxN*e->GetGeom()->GetEid(0) + n;
585  }
586 
587  // Edge 1
588  edgeOrient = e->GetGeom()->GetEorient(1);
589  if (edgeOrient==StdRegions::eForwards)
590  {
591  vId[cnt1+np-1+n*np] = 4*nel +
592  maxN*e->GetGeom()->GetEid(1) + n;
593  }
594  else
595  {
596  vId[cnt1+np*np-1-n*np] = 4*nel +
597  maxN*e->GetGeom()->GetEid(1) + n;
598  }
599 
600  // Edge 2
601  edgeOrient = e->GetGeom()->GetEorient(2);
602  if (edgeOrient==StdRegions::eForwards)
603  {
604  vId[cnt1+np*np-1-n] = 4*nel +
605  maxN*e->GetGeom()->GetEid(2) + n;
606  }
607  else
608  {
609  vId[cnt1+np*(np-1)+n] = 4*nel +
610  maxN*e->GetGeom()->GetEid(2) + n;
611  }
612 
613  // Edge 3
614  edgeOrient = e->GetGeom()->GetEorient(3);
615  if (edgeOrient==StdRegions::eForwards)
616  {
617  vId[cnt1+np*(np-1)-n*np] = 4*nel +
618  maxN*e->GetGeom()->GetEid(3) + n;
619  }
620  else
621  {
622  vId[cnt1+n*np] = 4*nel +
623  maxN*e->GetGeom()->GetEid(3) + n;
624  }
625  }
626 
627  // Interior numbering
628  for (int n = 1; n < np-1; n++)
629  {
630  for (int m = 1; m < np-1; m++)
631  {
632  vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
633  cnt2++;
634  }
635  }
636  cnt1+= newpoints;
637  }
638  break;
639  default:
640  {
641  ASSERTL0(false,"Points not known");
642  }
643  }
644  }
645 
646  // Crete new connectivity using homogeneous information
647  vector<Array<OneD, int> > oldConn;
648  vector<Array<OneD, int> > newConn;
649  Array<OneD, int> conn;
650  m_f->m_fieldPts->GetConnectivity(oldConn);
651 
652  for(int i = 0; i < nel; ++i)
653  {
654  int nTris = oldConn[i].num_elements()/3;
655  // Create array for new connectivity
656  // (3 tetrahedra between each plane for each triangle)
657  conn = Array<OneD, int> (4*3*nTris*(nPlanes-1));
658  cnt2=0;
659  for(int n = 0; n<nTris; n++)
660  {
661  // Get id of vertices in this triangle (on plane 0)
662  int vId0 = vId[oldConn[i][3*n+0]];
663  int vId1 = vId[oldConn[i][3*n+1]];
664  int vId2 = vId[oldConn[i][3*n+2]];
665 
666  // Determine ordering based on global Id of pts
667  Array<OneD, int> rot(3);
668  if ( (vId0<vId1) && (vId0<vId2))
669  {
670  rot[0] = 0;
671  if (vId1<vId2)
672  {
673  rot[1] = 1;
674  rot[2] = 2;
675  }
676  else
677  {
678  rot[1] = 2;
679  rot[2] = 1;
680  }
681  }
682  else if ( (vId1<vId0) && (vId1<vId2))
683  {
684  rot[0] = 1;
685  if (vId0<vId2)
686  {
687  rot[1] = 0;
688  rot[2] = 2;
689  }
690  else
691  {
692  rot[1] = 2;
693  rot[2] = 0;
694  }
695  }
696  else
697  {
698  rot[0] = 2;
699  if (vId0<vId1)
700  {
701  rot[1] = 0;
702  rot[2] = 1;
703  }
704  else
705  {
706  rot[1] = 1;
707  rot[2] = 0;
708  }
709  }
710 
711  // Define new connectivity with tetrahedra
712  for ( int p = 0; p<nPlanes-1; p++)
713  {
714  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
715  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
716  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[2]];
717  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
718 
719  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
720  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
721  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
722  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
723 
724  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
725  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
726  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
727  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
728  }
729  }
730  newConn.push_back(conn);
731  }
732  m_f->m_fieldPts->SetConnectivity(newConn);
733 }
734 
735 
737  int n,
738  const Array<OneD,const NekDouble> &phys,
739  Array<OneD, NekDouble> &coeffs)
740  {
742  e = m_f->m_exp[0]->GetExp(n);
743 
744  switch(e->DetShapeType())
745  {
747  {
748  int np0 = e->GetBasis(0)->GetNumPoints();
749  int np1 = e->GetBasis(1)->GetNumPoints();
750  int np = max(np0,np1);
751 
752  // to ensure points are correctly projected to np need
753  // to increase the order slightly of coordinates
754  LibUtilities::PointsKey pa(np+1,e->GetPointsType(0));
755  LibUtilities::PointsKey pb(np,e->GetPointsType(1));
756  Array<OneD, NekDouble> tophys(np*(np+1));
757 
760  StdRegions::StdTriExp OrthoExp(Ba,Bb);
761 
762  // interpolate points to new phys points!
763  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
764  e->GetBasis(1)->GetBasisKey(),
765  phys,Ba,Bb,tophys);
766 
767  OrthoExp.FwdTrans(tophys,coeffs);
768  break;
769  }
771  {
772  int np0 = e->GetBasis(0)->GetNumPoints();
773  int np1 = e->GetBasis(1)->GetNumPoints();
774  int np = max(np0,np1);
775 
776  LibUtilities::PointsKey pa(np+1,e->GetPointsType(0));
777  LibUtilities::PointsKey pb(np+1,e->GetPointsType(1));
778  Array<OneD, NekDouble> tophys((np+1)*(np+1));
779 
782  StdRegions::StdQuadExp OrthoExp(Ba,Bb);
783 
784  // interpolate points to new phys points!
785  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
786  e->GetBasis(1)->GetBasisKey(),
787  phys,Ba,Bb,tophys);
788 
789  OrthoExp.FwdTrans(phys,coeffs);
790  break;
791  }
792  default:
793  ASSERTL0(false,"Shape needs setting up");
794  break;
795  }
796  }
797 
798 }
799 }
800 
801 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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)
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
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
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
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...