Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | List of all members
Nektar::Utilities::ProcessEquiSpacedOutput Class Reference

This processing module interpolates one field to another. More...

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::Utilities::ProcessEquiSpacedOutput:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessEquiSpacedOutput:
Collaboration graph
[legend]

Public Member Functions

 ProcessEquiSpacedOutput (FieldSharedPtr f)
 
virtual ~ProcessEquiSpacedOutput ()
 
virtual void Process (po::variables_map &vm)
 Write mesh to output file. More...
 
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
 ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
 
void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
void PrintConfig ()
 Print out all configuration options for a module. More...
 
void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
bool GetRequireEquiSpaced (void)
 
void SetRequireEquiSpaced (bool pVal)
 
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 
 Module (MeshSharedPtr p_m)
 
virtual void Process ()=0
 
void RegisterConfig (std::string key, std::string value)
 
void PrintConfig ()
 
void SetDefaults ()
 
MeshSharedPtr GetMesh ()
 
virtual void ProcessVertices ()
 Extract element vertices. More...
 
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual void ProcessElements ()
 Generate element IDs. More...
 
virtual void ProcessComposites ()
 Generate composites. More...
 
virtual void ClearElementLinks ()
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Protected Member Functions

 ProcessEquiSpacedOutput ()
 
void SetupEquiSpacedField (void)
 
void SetHomogeneousConnectivity (void)
 
void GenOrthoModes (int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
 
- Protected Member Functions inherited from Nektar::Utilities::Module
 Module ()
 
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 

Additional Inherited Members

- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::map< std::string,
ConfigOption
m_config
 List of configuration values. More...
 

Detailed Description

This processing module interpolates one field to another.

Definition at line 49 of file ProcessEquiSpacedOutput.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessEquiSpacedOutput::ProcessEquiSpacedOutput ( FieldSharedPtr  f)

Definition at line 60 of file ProcessEquiSpacedOutput.cpp.

References Nektar::Utilities::Module::m_config.

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 }
map< string, ConfigOption > m_config
List of configuration values.
Nektar::Utilities::ProcessEquiSpacedOutput::~ProcessEquiSpacedOutput ( )
virtual

Definition at line 72 of file ProcessEquiSpacedOutput.cpp.

73 {
74 }
Nektar::Utilities::ProcessEquiSpacedOutput::ProcessEquiSpacedOutput ( )
inlineprotected

Definition at line 65 of file ProcessEquiSpacedOutput.h.

65 {};

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessEquiSpacedOutput::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 53 of file ProcessEquiSpacedOutput.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

53  {
56  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Nektar::Utilities::ProcessEquiSpacedOutput::GenOrthoModes ( int  n,
const Array< OneD, const NekDouble > &  phys,
Array< OneD, NekDouble > &  coeffs 
)
protected

Definition at line 736 of file ProcessEquiSpacedOutput.cpp.

References ASSERTL0, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::LibUtilities::Interp2D(), and Nektar::Utilities::Module::m_f.

Referenced by SetupEquiSpacedField().

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 
758  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,np,pa);
759  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,np,pb);
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 
780  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,np,pa);
781  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,np,pb);
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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
Principle Orthogonal Functions .
Definition: BasisType.h:46
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void Nektar::Utilities::ProcessEquiSpacedOutput::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Reimplemented in Nektar::Utilities::ProcessIsoContour.

Definition at line 76 of file ProcessEquiSpacedOutput.cpp.

References SetupEquiSpacedField().

void Nektar::Utilities::ProcessEquiSpacedOutput::SetHomogeneousConnectivity ( void  )
protected

Definition at line 456 of file ProcessEquiSpacedOutput.cpp.

References ASSERTL0, Nektar::StdRegions::eForwards, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::Utilities::Module::m_f, and npts.

Referenced by SetupEquiSpacedField().

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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
FieldSharedPtr m_f
Field object.
static std::string npts
Definition: InputFld.cpp:43
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
void Nektar::Utilities::ProcessEquiSpacedOutput::SetupEquiSpacedField ( void  )
protected

Definition at line 82 of file ProcessEquiSpacedOutput.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::StdRegions::eBwd, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eFwd, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, GenOrthoModes(), Nektar::LibUtilities::StdSegData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::LibUtilities::StdHexData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTetData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, npts, Nektar::NullNekDouble1DArray, and SetHomogeneousConnectivity().

Referenced by Process(), and Nektar::Utilities::ProcessIsoContour::Process().

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();
93  Array<OneD, Array<OneD, NekDouble> > coords(3);
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 
386  Array<OneD, NekDouble> tmp;
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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:286
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
map< string, ConfigOption > m_config
List of configuration values.
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.
static std::string npts
Definition: InputFld.cpp:43
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
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

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessEquiSpacedOutput::className
static
Initial value:
=
ModuleKey(eProcessModule, "equispacedoutput"),
"Write data as equi-spaced output using simplices to represent the data for connecting points")

Definition at line 57 of file ProcessEquiSpacedOutput.h.