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...
 
virtual std::string GetModuleName ()
 
- 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 71 of file ProcessEquiSpacedOutput.h.

71 {};

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 808 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().

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 
830  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,np,pa);
831  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,np,pb);
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 
852  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,np,pa);
853  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,np,pb);
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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
virtual std::string Nektar::Utilities::ProcessEquiSpacedOutput::GetModuleName ( )
inlinevirtual

Implements Nektar::Utilities::Module.

Reimplemented in Nektar::Utilities::ProcessIsoContour.

Definition at line 65 of file ProcessEquiSpacedOutput.h.

66  {
67  return "ProcessEquiSpacedOutput";
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 458 of file ProcessEquiSpacedOutput.cpp.

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

Referenced by SetupEquiSpacedField().

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
472  Array<OneD, Array<OneD, NekDouble> > pts;
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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Fourier Expansion .
Definition: BasisType.h:52
FieldSharedPtr m_f
Field object.
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
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
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
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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  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();
95  Array<OneD, Array<OneD, NekDouble> > coords(3);
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 
388  Array<OneD, NekDouble> tmp;
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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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.