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

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

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::FieldUtils::ProcessEquiSpacedOutput:
Inheritance graph
[legend]
Collaboration diagram for Nektar::FieldUtils::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::FieldUtils::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
- Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT Module (FieldSharedPtr p_f)
 
FIELD_UTILS_EXPORT void RegisterConfig (string key, string value)
 Register a configuration option with a module. More...
 
FIELD_UTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
FIELD_UTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
FIELD_UTILS_EXPORT bool GetRequireEquiSpaced (void)
 
FIELD_UTILS_EXPORT void SetRequireEquiSpaced (bool pVal)
 
FIELD_UTILS_EXPORT void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 

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::FieldUtils::Module
 Module ()
 

Additional Inherited Members

- Protected Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr m_f
 Field object. More...
 
map< string, ConfigOptionm_config
 List of configuration values. More...
 
bool m_requireEquiSpaced
 

Detailed Description

This processing module interpolates one field to another.

Definition at line 49 of file ProcessEquiSpacedOutput.h.

Constructor & Destructor Documentation

Nektar::FieldUtils::ProcessEquiSpacedOutput::ProcessEquiSpacedOutput ( FieldSharedPtr  f)

Definition at line 60 of file ProcessEquiSpacedOutput.cpp.

References Nektar::FieldUtils::Module::m_config.

61  : ProcessModule(f)
62 {
63  f->m_setUpEquiSpacedFields = true;
64 
65  m_config["tetonly"] =
66  ConfigOption(true, "NotSet", "Only process tetrahedral elements");
67 
68  m_config["modalenergy"] =
69  ConfigOption(true, "NotSet", "Write output as modal energy");
70 }
map< string, ConfigOption > m_config
List of configuration values.
Nektar::FieldUtils::ProcessEquiSpacedOutput::~ProcessEquiSpacedOutput ( )
virtual

Definition at line 72 of file ProcessEquiSpacedOutput.cpp.

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

Definition at line 71 of file ProcessEquiSpacedOutput.h.

71 {};

Member Function Documentation

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

Creates an instance of this class.

Definition at line 53 of file ProcessEquiSpacedOutput.h.

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

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

Definition at line 848 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::FieldUtils::Module::m_f.

Referenced by SetupEquiSpacedField().

852 {
854  e = m_f->m_exp[0]->GetExp(n);
855 
856  switch (e->DetShapeType())
857  {
859  {
860  int np0 = e->GetBasis(0)->GetNumPoints();
861  int np1 = e->GetBasis(1)->GetNumPoints();
862  int np = max(np0, np1);
863 
864  // to ensure points are correctly projected to np need
865  // to increase the order slightly of coordinates
866  LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
867  LibUtilities::PointsKey pb(np, e->GetPointsType(1));
868  Array<OneD, NekDouble> tophys(np * (np + 1));
869 
870  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
871  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, np, pb);
872  StdRegions::StdTriExp OrthoExp(Ba, Bb);
873 
874  // interpolate points to new phys points!
875  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
876  e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
877  tophys);
878 
879  OrthoExp.FwdTrans(tophys, coeffs);
880  break;
881  }
883  {
884  int np0 = e->GetBasis(0)->GetNumPoints();
885  int np1 = e->GetBasis(1)->GetNumPoints();
886  int np = max(np0, np1);
887 
888  LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
889  LibUtilities::PointsKey pb(np + 1, e->GetPointsType(1));
890  Array<OneD, NekDouble> tophys((np + 1) * (np + 1));
891 
892  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
893  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, np, pb);
894  StdRegions::StdQuadExp OrthoExp(Ba, Bb);
895 
896  // interpolate points to new phys points!
897  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
898  e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
899  tophys);
900 
901  OrthoExp.FwdTrans(phys, coeffs);
902  break;
903  }
904  default:
905  ASSERTL0(false, "Shape needs setting up");
906  break;
907  }
908 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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
FieldSharedPtr m_f
Field object.
virtual std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::GetModuleName ( )
inlinevirtual

Implements Nektar::FieldUtils::Module.

Definition at line 65 of file ProcessEquiSpacedOutput.h.

66  {
67  return "ProcessEquiSpacedOutput";
68  }
void Nektar::FieldUtils::ProcessEquiSpacedOutput::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::FieldUtils::Module.

Reimplemented in Nektar::FieldUtils::ProcessIsoContour.

Definition at line 76 of file ProcessEquiSpacedOutput.cpp.

References SetupEquiSpacedField().

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

Definition at line 461 of file ProcessEquiSpacedOutput.cpp.

References ASSERTL0, Nektar::StdRegions::eForwards, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::FieldUtils::Module::m_f, npts, CellMLToNektar.cellml_metadata::p, Vmath::Sadd(), and Vmath::Vcopy().

Referenced by SetupEquiSpacedField().

462 {
464  int nel = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
465  int nPlanes = m_f->m_exp[0]->GetZIDs().num_elements();
466  int npts = m_f->m_fieldPts->GetNpoints();
467  int nptsPerPlane = npts / nPlanes;
468  int coordim = 3;
469 
470  if (m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
472  m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
474  {
475  // Write points with extra plane
476  Array<OneD, Array<OneD, NekDouble> > pts;
477  m_f->m_fieldPts->GetPts(pts);
478  Array<OneD, Array<OneD, NekDouble> > newPts(pts.num_elements());
479  for (int i = 0; i < pts.num_elements(); i++)
480  {
481  newPts[i] = Array<OneD, NekDouble>(npts + nptsPerPlane);
482  // Copy old points
483  Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
484  // Get new plane
485  Array<OneD, NekDouble> extraPlane(nptsPerPlane, newPts[i] + npts);
486  if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
487  {
488  if ( i == (coordim-1))
489  {
490  // z-coordinate
491  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
492  pts[i], 1, extraPlane, 1);
493  }
494  else
495  {
496  Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
497  }
498  }
499  else
500  {
501  // Determine to and from rank for communication
502  int size = m_f->m_comm->GetColumnComm()->GetSize();
503  int rank = m_f->m_comm->GetColumnComm()->GetRank();
504  int fromRank = (rank + 1) % size;
505  int toRank = (rank == 0) ? size - 1 : rank - 1;
506  // Communicate using SendRecv
507  Array<OneD, NekDouble> send(nptsPerPlane, pts[i]);
508  m_f->m_comm->GetColumnComm()->SendRecv(toRank, send, fromRank,
509  extraPlane);
510  // Adjust z-coordinate of last process
511  if (i == (coordim-1) && (rank == size - 1) )
512  {
513  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
514  extraPlane, 1, extraPlane, 1);
515  }
516  }
517  }
518  m_f->m_fieldPts->SetPts(newPts);
519  // Now extend plane connectivity
520  vector<Array<OneD, int> > oldConn;
521  Array<OneD, int> conn;
522  m_f->m_fieldPts->GetConnectivity(oldConn);
523  vector<Array<OneD, int> > newConn = oldConn;
524  int connPerPlane = oldConn.size() / nPlanes;
525  for (int i = 0; i < connPerPlane; i++)
526  {
527  conn = Array<OneD, int>(oldConn[i].num_elements());
528  for (int j = 0; j < conn.num_elements(); j++)
529  {
530  conn[j] = oldConn[i][j] + npts;
531  }
532  newConn.push_back(conn);
533  }
534  m_f->m_fieldPts->SetConnectivity(newConn);
535 
536  nPlanes++;
537  npts += nptsPerPlane;
538  }
539 
540  vector<Array<OneD, int> > oldConn;
541  vector<Array<OneD, int> > newConn;
542  Array<OneD, int> conn;
543  m_f->m_fieldPts->GetConnectivity(oldConn);
544 
545  // 2DH1D case
546  if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock)
547  {
548  for(int i = 0; i < nel; ++i)
549  {
550  int nLines = oldConn[i].num_elements()/2;
551  // Create array for new connectivity
552  // (2 triangles between each plane for each line)
553  conn = Array<OneD, int> (2*3*nLines*(nPlanes-1));
554  int cnt = 0;
555  for(int n = 0; n<nLines; n++)
556  {
557  // Define new connectivity with triangles
558  for ( int p = 0; p<nPlanes-1; p++)
559  {
560  conn[cnt++] = oldConn[i+ p *nel][2*n+0];
561  conn[cnt++] = oldConn[i+ p *nel][2*n+1];
562  conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];
563 
564  conn[cnt++] = oldConn[i+ p *nel][2*n+1];
565  conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];
566  conn[cnt++] = oldConn[i+(p+1)*nel][2*n+1];
567  }
568  }
569  newConn.push_back(conn);
570  }
571  }
572  else if(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock)
573  {
574  // Get maximum number of points per direction in the mesh
575  int maxN = 0;
576  for(int i = 0; i < nel; ++i)
577  {
578  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
579 
580  int np0 = e->GetBasis(0)->GetNumPoints();
581  int np1 = e->GetBasis(1)->GetNumPoints();
582  int np = max(np0,np1);
583  maxN = max(np, maxN);
584  }
585 
586  // Create a global numbering for points in plane 0
587  Array<OneD, int> vId(nptsPerPlane);
588  int cnt1=0;
589  int cnt2=0;
590  for(int i = 0; i < nel; ++i)
591  {
592  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
593  int np0 = e->GetBasis(0)->GetNumPoints();
594  int np1 = e->GetBasis(1)->GetNumPoints();
595  int np = max(np0,np1);
596  switch(e->DetShapeType())
597  {
599  {
600  int newpoints = LibUtilities::StdTriData::
602 
603  // Vertex numbering
604  vId[cnt1] = e->GetGeom()->GetVid(0);
605  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
606  vId[cnt1+newpoints-1] = e->GetGeom()->GetVid(2);
607 
608  // Edge numbering
609  StdRegions::Orientation edgeOrient;
610  int edge1 = 0;
611  int edge2 = 0;
612  for (int n = 1; n < np-1; n++)
613  {
614  // Edge 0
615  edgeOrient = e->GetGeom()->GetEorient(0);
616  if (edgeOrient==StdRegions::eForwards)
617  {
618  vId[cnt1+n] = 4*nel +
619  maxN*e->GetGeom()->GetEid(0) + n;
620  }
621  else
622  {
623  vId[cnt1+np-1-n] = 4*nel +
624  maxN*e->GetGeom()->GetEid(0) + n;
625  }
626 
627  // Edge 1
628  edgeOrient = e->GetGeom()->GetEorient(1);
629  if (edgeOrient==StdRegions::eForwards)
630  {
631  edge1 += np-n;
632  vId[cnt1+np-1+edge1] = 4*nel +
633  maxN*e->GetGeom()->GetEid(1) + n;
634  }
635  else
636  {
637  edge1 += n;
638  vId[cnt1+newpoints-1-edge1] = 4*nel +
639  maxN*e->GetGeom()->GetEid(1) + n;
640  }
641 
642  // Edge 2
643  edgeOrient = e->GetGeom()->GetEorient(2);
644  if (edgeOrient==StdRegions::eForwards)
645  {
646  edge2 += n+1;
647  vId[cnt1+newpoints-1-edge2] = 4*nel +
648  maxN*e->GetGeom()->GetEid(2) + n;
649  }
650  else
651  {
652  edge2 += np+1-n;
653  vId[cnt1+edge2] = 4*nel +
654  maxN*e->GetGeom()->GetEid(2) + n;
655  }
656  }
657 
658  // Interior numbering
659  edge2 = 0;
660  for (int n = 1; n < np-1; n++)
661  {
662  edge2 += np+1-n;
663  for (int m = 1; m < np-n-1; m++)
664  {
665  vId[cnt1+edge2+m] = 4*nel + maxN*4*nel + cnt2;
666  cnt2++;
667  }
668  }
669  cnt1+= newpoints;
670  }
671  break;
673  {
674  int newpoints = LibUtilities::StdQuadData::
676  // Vertex numbering
677  vId[cnt1] = e->GetGeom()->GetVid(0);
678  vId[cnt1+np-1] = e->GetGeom()->GetVid(1);
679  vId[cnt1+np*np-1] = e->GetGeom()->GetVid(2);
680  vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);
681 
682  // Edge numbering
683  StdRegions::Orientation edgeOrient;
684  for (int n = 1; n < np-1; n++)
685  {
686  // Edge 0
687  edgeOrient = e->GetGeom()->GetEorient(0);
688  if (edgeOrient==StdRegions::eForwards)
689  {
690  vId[cnt1+n] = 4*nel +
691  maxN*e->GetGeom()->GetEid(0) + n;
692  }
693  else
694  {
695  vId[cnt1+np-1-n] = 4*nel +
696  maxN*e->GetGeom()->GetEid(0) + n;
697  }
698 
699  // Edge 1
700  edgeOrient = e->GetGeom()->GetEorient(1);
701  if (edgeOrient==StdRegions::eForwards)
702  {
703  vId[cnt1+np-1+n*np] = 4*nel +
704  maxN*e->GetGeom()->GetEid(1) + n;
705  }
706  else
707  {
708  vId[cnt1+np*np-1-n*np] = 4*nel +
709  maxN*e->GetGeom()->GetEid(1) + n;
710  }
711 
712  // Edge 2
713  edgeOrient = e->GetGeom()->GetEorient(2);
714  if (edgeOrient==StdRegions::eForwards)
715  {
716  vId[cnt1+np*np-1-n] = 4*nel +
717  maxN*e->GetGeom()->GetEid(2) + n;
718  }
719  else
720  {
721  vId[cnt1+np*(np-1)+n] = 4*nel +
722  maxN*e->GetGeom()->GetEid(2) + n;
723  }
724 
725  // Edge 3
726  edgeOrient = e->GetGeom()->GetEorient(3);
727  if (edgeOrient==StdRegions::eForwards)
728  {
729  vId[cnt1+np*(np-1)-n*np] = 4*nel +
730  maxN*e->GetGeom()->GetEid(3) + n;
731  }
732  else
733  {
734  vId[cnt1+n*np] = 4*nel +
735  maxN*e->GetGeom()->GetEid(3) + n;
736  }
737  }
738 
739  // Interior numbering
740  for (int n = 1; n < np-1; n++)
741  {
742  for (int m = 1; m < np-1; m++)
743  {
744  vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
745  cnt2++;
746  }
747  }
748  cnt1+= newpoints;
749  }
750  break;
751  default:
752  {
753  ASSERTL0(false,"Points not known");
754  }
755  }
756  }
757 
758  // Crete new connectivity using homogeneous information
759  for(int i = 0; i < nel; ++i)
760  {
761  int nTris = oldConn[i].num_elements()/3;
762  // Create array for new connectivity
763  // (3 tetrahedra between each plane for each triangle)
764  conn = Array<OneD, int> (4*3*nTris*(nPlanes-1));
765  cnt2=0;
766  for(int n = 0; n<nTris; n++)
767  {
768  // Get id of vertices in this triangle (on plane 0)
769  int vId0 = vId[oldConn[i][3*n+0]];
770  int vId1 = vId[oldConn[i][3*n+1]];
771  int vId2 = vId[oldConn[i][3*n+2]];
772 
773  // Determine ordering based on global Id of pts
774  Array<OneD, int> rot(3);
775  if ( (vId0<vId1) && (vId0<vId2))
776  {
777  rot[0] = 0;
778  if (vId1<vId2)
779  {
780  rot[1] = 1;
781  rot[2] = 2;
782  }
783  else
784  {
785  rot[1] = 2;
786  rot[2] = 1;
787  }
788  }
789  else if ( (vId1<vId0) && (vId1<vId2))
790  {
791  rot[0] = 1;
792  if (vId0<vId2)
793  {
794  rot[1] = 0;
795  rot[2] = 2;
796  }
797  else
798  {
799  rot[1] = 2;
800  rot[2] = 0;
801  }
802  }
803  else
804  {
805  rot[0] = 2;
806  if (vId0<vId1)
807  {
808  rot[1] = 0;
809  rot[2] = 1;
810  }
811  else
812  {
813  rot[1] = 1;
814  rot[2] = 0;
815  }
816  }
817 
818  // Define new connectivity with tetrahedra
819  for ( int p = 0; p<nPlanes-1; p++)
820  {
821  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
822  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
823  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[2]];
824  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
825 
826  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
827  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[1]];
828  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
829  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
830 
831  conn[cnt2++] = oldConn[i+ p *nel][3*n+rot[0]];
832  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
833  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
834  conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
835  }
836  }
837  newConn.push_back(conn);
838  }
839  }
840  else
841  {
842  ASSERTL0( false, "Points type not supported.");
843  }
844 
845  m_f->m_fieldPts->SetConnectivity(newConn);
846 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Fourier Expansion .
Definition: BasisType.h:52
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:66
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:315
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:1061
FieldSharedPtr m_f
Field object.
void Nektar::FieldUtils::ProcessEquiSpacedOutput::SetupEquiSpacedField ( void  )
protected

Definition at line 81 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::ePtsSegBlock, 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::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, npts, Nektar::NullNekDouble1DArray, Nektar::LibUtilities::NullPtsField, and SetHomogeneousConnectivity().

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

82 {
83  if (m_f->m_verbose)
84  {
85  if (m_f->m_comm->TreatAsRankZero())
86  {
87  cout << "Interpolating fields to equispaced..." << endl;
88  }
89  }
90  int nel = m_f->m_exp[0]->GetExpSize();
91  if (!nel)
92  {
93  m_f->m_fieldPts = LibUtilities::NullPtsField;
94  return;
95  }
96 
97  int coordim = m_f->m_exp[0]->GetCoordim(0);
98  int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
99  int npts = m_f->m_exp[0]->GetTotPoints();
100  Array<OneD, Array<OneD, NekDouble> > coords(3);
101 
102  // Check if we have a homogeneous expansion
103  bool homogeneous1D = false;
104  if (m_f->m_fielddef.size())
105  {
106  if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
107  {
108  coordim++;
109  shapedim++;
110  homogeneous1D = true;
111  }
112  else if (m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
113  {
114  ASSERTL0(false, "Homegeneous2D case not supported");
115  }
116  }
117  else
118  {
119  if (m_f->m_session->DefinesSolverInfo("HOMOGENEOUS"))
120  {
121  std::string HomoStr = m_f->m_session->GetSolverInfo("HOMOGENEOUS");
122 
123  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
124  (HomoStr == "1D") || (HomoStr == "Homo1D"))
125  {
126  coordim++;
127  shapedim++;
128  homogeneous1D = true;
129  }
130  if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
131  (HomoStr == "2D") || (HomoStr == "Homo2D"))
132  {
133  ASSERTL0(false, "Homegeneous2D case not supported");
134  }
135  }
136  }
137 
138  // set up the number of points in each element
139  int newpoints = 0;
140  int newtotpoints = 0;
141 
142  Array<OneD, int> conn;
143  int prevNcoeffs = 0;
144  int prevNpoints = 0;
145  int cnt = 0;
146 
147  // identify face 1 connectivity for prisms
148  map<int, StdRegions::Orientation> face0orient;
149  set<int> prismorient;
151 
152  // prepare PtsField
153  vector<std::string> fieldNames;
154  vector<int> ppe;
155  vector<Array<OneD, int> > ptsConn;
156  int nfields;
157 
158  for (int i = 0; i < nel; ++i)
159  {
160  e = m_f->m_exp[0]->GetExp(i);
161  if (e->DetShapeType() == LibUtilities::ePrism)
162  {
163  StdRegions::Orientation forient = e->GetForient(0);
164  int fid = e->GetGeom()->GetFid(0);
165  if (face0orient.count(fid))
166  { // face 1 meeting face 1 so reverse this id
167  prismorient.insert(i);
168  }
169  else
170  {
171  // just store if Dir 1 is fwd or bwd
172  if ((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
176  {
177  face0orient[fid] = StdRegions::eBwd;
178  }
179  else
180  {
181  face0orient[fid] = StdRegions::eFwd;
182  }
183  }
184  }
185  }
186 
187  for (int i = 0; i < nel; ++i)
188  {
189  e = m_f->m_exp[0]->GetExp(i);
190  if (e->DetShapeType() == LibUtilities::ePrism)
191  {
192  int fid = e->GetGeom()->GetFid(2);
193  // check to see if face 2 meets face 1
194  if (face0orient.count(fid))
195  {
196  // check to see how face 2 is orientated
197  StdRegions::Orientation forient2 = e->GetForient(2);
198  StdRegions::Orientation forient0 = face0orient[fid];
199 
200  // If dir 1 or forient2 is bwd then check agains
201  // face 1 value
202  if ((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
203  (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
204  (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
206  {
207  if (forient0 == StdRegions::eFwd)
208  {
209  prismorient.insert(i);
210  }
211  }
212  else
213  {
214  if (forient0 == StdRegions::eBwd)
215  {
216  prismorient.insert(i);
217  }
218  }
219  }
220  }
221  }
222 
223  for (int i = 0; i < nel; ++i)
224  {
225  e = m_f->m_exp[0]->GetExp(i);
226  if (m_config["tetonly"].m_beenSet)
227  {
228  if (m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
230  {
231  continue;
232  }
233  }
234 
235  switch (e->DetShapeType())
236  {
238  {
239  int npoints0 = e->GetBasis(0)->GetNumPoints();
240 
241  newpoints =
243  }
244  break;
246  {
247  int np0 = e->GetBasis(0)->GetNumPoints();
248  int np1 = e->GetBasis(1)->GetNumPoints();
249  int np = max(np0, np1);
250  newpoints =
252  }
253  break;
255  {
256  int np0 = e->GetBasis(0)->GetNumPoints();
257  int np1 = e->GetBasis(1)->GetNumPoints();
258  int np = max(np0, np1);
259 
260  newpoints =
262  }
263  break;
265  {
266  int np0 = e->GetBasis(0)->GetNumPoints();
267  int np1 = e->GetBasis(1)->GetNumPoints();
268  int np2 = e->GetBasis(2)->GetNumPoints();
269  int np = max(np0, max(np1, np2));
270 
272  np, np, np);
273  }
274  break;
276  {
277  int np0 = e->GetBasis(0)->GetNumPoints();
278  int np1 = e->GetBasis(1)->GetNumPoints();
279  int np2 = e->GetBasis(2)->GetNumPoints();
280  int np = max(np0, max(np1, np2));
281 
283  np, np, np);
284  }
285  break;
287  {
288  int np0 = e->GetBasis(0)->GetNumPoints();
289  int np1 = e->GetBasis(1)->GetNumPoints();
290  int np2 = e->GetBasis(2)->GetNumPoints();
291  int np = max(np0, max(np1, np2));
292 
294  np, np, np);
295  }
296  break;
298  {
299  int np0 = e->GetBasis(0)->GetNumPoints();
300  int np1 = e->GetBasis(1)->GetNumPoints();
301  int np2 = e->GetBasis(2)->GetNumPoints();
302  int np = max(np0, max(np1, np2));
303 
305  np, np, np);
306  }
307  break;
308  default:
309  {
310  ASSERTL0(false, "Points not known");
311  }
312  }
313 
314  ppe.push_back(newpoints);
315  newtotpoints += newpoints;
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  Array<OneD, NekDouble> tmp;
388 
389  for (int n = 0; n < coordim; ++n)
390  {
391  cnt = 0;
392  int cnt1 = 0;
393  for (int i = 0; i < nel; ++i)
394  {
395  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
396  coords[n] + cnt, 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, tmp = pts[coordim + n] + cnt1);
429  cnt1 += ppe[i];
430  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
431  }
432  }
433 
434  // Set up Variable string.
435  fieldNames.push_back(m_f->m_fielddef[0]->m_fields[n]);
436  }
437  }
438 
440  coordim, fieldNames, pts);
441  if (shapedim == 1)
442  {
443  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsSegBlock);
444  }
445 
446  if (shapedim == 2)
447  {
448  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
449  }
450  else if (shapedim == 3)
451  {
452  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
453  }
454  m_f->m_fieldPts->SetConnectivity(ptsConn);
455  if (homogeneous1D)
456  {
458  }
459 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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.
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:186
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 GenOrthoModes(int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
static PtsFieldSharedPtr NullPtsField
Definition: PtsField.h:179
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:151
FieldSharedPtr m_f
Field object.

Member Data Documentation

ModuleKey Nektar::FieldUtils::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.