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

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::FieldUtils::ProcessEquiSpacedOutput:

Collaboration diagram for Nektar::FieldUtils::ProcessEquiSpacedOutput:

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< Module >	create (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, ConfigOption >	m_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.

     : ProcessModule(f)
 {
     f->m_setUpEquiSpacedFields = true;
 
     m_config["tetonly"] =
         ConfigOption(true, "NotSet", "Only process tetrahedral elements");
 
     m_config["modalenergy"] =
         ConfigOption(true, "NotSet", "Write output as modal energy");
 }

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

     {
         return MemoryManager<ProcessEquiSpacedOutput>::AllocateSharedPtr(f);
     }

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

 {
     LocalRegions::ExpansionSharedPtr e;
     e = m_f->m_exp[0]->GetExp(n);
 
     switch (e->DetShapeType())
     {
         case LibUtilities::eTriangle:
         {
             int np0 = e->GetBasis(0)->GetNumPoints();
             int np1 = e->GetBasis(1)->GetNumPoints();
             int np  = max(np0, np1);
 
             // to ensure points are correctly projected to np need
             // to increase the order slightly of coordinates
             LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
             LibUtilities::PointsKey pb(np, e->GetPointsType(1));
             Array<OneD, NekDouble> tophys(np * (np + 1));
 
             LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
             LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, np, pb);
             StdRegions::StdTriExp OrthoExp(Ba, Bb);
 
             // interpolate points to new phys points!
             LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
                                    e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
                                    tophys);
 
             OrthoExp.FwdTrans(tophys, coeffs);
             break;
         }
         case LibUtilities::eQuadrilateral:
         {
             int np0 = e->GetBasis(0)->GetNumPoints();
             int np1 = e->GetBasis(1)->GetNumPoints();
             int np  = max(np0, np1);
 
             LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
             LibUtilities::PointsKey pb(np + 1, e->GetPointsType(1));
             Array<OneD, NekDouble> tophys((np + 1) * (np + 1));
 
             LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
             LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, np, pb);
             StdRegions::StdQuadExp OrthoExp(Ba, Bb);
 
             // interpolate points to new phys points!
             LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
                                    e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
                                    tophys);
 
             OrthoExp.FwdTrans(phys, coeffs);
             break;
         }
         default:
             ASSERTL0(false, "Shape needs setting up");
             break;
     }
 }

virtual std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::GetModuleName ( )

inlinevirtual

Implements Nektar::FieldUtils::Module.

Definition at line 65 of file ProcessEquiSpacedOutput.h.

     {
         return "ProcessEquiSpacedOutput";
     }

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

 {
     SetupEquiSpacedField();
 }

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

protected

Definition at line 461 of file ProcessEquiSpacedOutput.cpp.

Referenced by SetupEquiSpacedField().

 {
     LocalRegions::ExpansionSharedPtr e;
     int nel          = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
     int nPlanes      = m_f->m_exp[0]->GetZIDs().num_elements();
     int npts         = m_f->m_fieldPts->GetNpoints();
     int nptsPerPlane = npts / nPlanes;
     int coordim      = 3;
 
     if (m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
             LibUtilities::eFourier &&
         m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
             LibUtilities::eFourierEvenlySpaced)
     {
         // Write points with extra plane
         Array<OneD, Array<OneD, NekDouble> > pts;
         m_f->m_fieldPts->GetPts(pts);
         Array<OneD, Array<OneD, NekDouble> > newPts(pts.num_elements());
         for (int i = 0; i < pts.num_elements(); i++)
         {
             newPts[i] = Array<OneD, NekDouble>(npts + nptsPerPlane);
             // Copy old points
             Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
             // Get new plane
             Array<OneD, NekDouble> extraPlane(nptsPerPlane, newPts[i] + npts);
             if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
             {
                 if ( i == (coordim-1))
                 {
                     // z-coordinate
                     Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
                                 pts[i], 1, extraPlane, 1);
                 }
                 else
                 {
                     Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
                 }
             }
             else
             {
                 // Determine to and from rank for communication
                 int size     = m_f->m_comm->GetColumnComm()->GetSize();
                 int rank     = m_f->m_comm->GetColumnComm()->GetRank();
                 int fromRank = (rank + 1) % size;
                 int toRank   = (rank == 0) ? size - 1 : rank - 1;
                 // Communicate using SendRecv
                 Array<OneD, NekDouble> send(nptsPerPlane, pts[i]);
                 m_f->m_comm->GetColumnComm()->SendRecv(toRank, send, fromRank,
                                                        extraPlane);
                 // Adjust z-coordinate of last process
                 if (i == (coordim-1) && (rank == size - 1) )
                 {
                     Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
                                 extraPlane, 1, extraPlane, 1);
                 }
             }
         }
         m_f->m_fieldPts->SetPts(newPts);
         // Now extend plane connectivity
         vector<Array<OneD, int> > oldConn;
         Array<OneD, int> conn;
         m_f->m_fieldPts->GetConnectivity(oldConn);
         vector<Array<OneD, int> > newConn = oldConn;
         int connPerPlane = oldConn.size() / nPlanes;
         for (int i = 0; i < connPerPlane; i++)
         {
             conn = Array<OneD, int>(oldConn[i].num_elements());
             for (int j = 0; j < conn.num_elements(); j++)
             {
                 conn[j] = oldConn[i][j] + npts;
             }
             newConn.push_back(conn);
         }
         m_f->m_fieldPts->SetConnectivity(newConn);
 
         nPlanes++;
         npts += nptsPerPlane;
     }
 
     vector<Array<OneD, int> > oldConn;
     vector<Array<OneD, int> > newConn;
     Array<OneD, int>          conn;
     m_f->m_fieldPts->GetConnectivity(oldConn);
 
     // 2DH1D case
     if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock)
     {
         for(int i = 0; i < nel; ++i)
         {
             int nLines = oldConn[i].num_elements()/2;
             // Create array for new connectivity
             // (2 triangles between each plane for each line)
             conn = Array<OneD, int> (2*3*nLines*(nPlanes-1));
             int cnt = 0;
             for(int n = 0; n<nLines; n++)
             {
                 // Define new connectivity with triangles
                 for ( int p = 0; p<nPlanes-1; p++)
                 {
                     conn[cnt++] = oldConn[i+  p  *nel][2*n+0];
                     conn[cnt++] = oldConn[i+  p  *nel][2*n+1];
                     conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];
 
                     conn[cnt++] = oldConn[i+  p  *nel][2*n+1];
                     conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];
                     conn[cnt++] = oldConn[i+(p+1)*nel][2*n+1];
                 }
             }
             newConn.push_back(conn);
         }
     }
     else if(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock)
     {
         // Get maximum number of points per direction in the mesh
         int maxN = 0;
         for(int i = 0; i < nel; ++i)
         {
             e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
 
             int np0 = e->GetBasis(0)->GetNumPoints();
             int np1 = e->GetBasis(1)->GetNumPoints();
             int np = max(np0,np1);
             maxN = max(np, maxN);
         }
 
         // Create a global numbering for points in plane 0
         Array<OneD, int> vId(nptsPerPlane);
         int cnt1=0;
         int cnt2=0;
         for(int i = 0; i < nel; ++i)
         {
             e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
             int np0 = e->GetBasis(0)->GetNumPoints();
             int np1 = e->GetBasis(1)->GetNumPoints();
             int np = max(np0,np1);
             switch(e->DetShapeType())
             {
             case LibUtilities::eTriangle:
                 {
                     int newpoints = LibUtilities::StdTriData::
                                         getNumberOfCoefficients(np,np);
 
                     // Vertex numbering
                     vId[cnt1]                 = e->GetGeom()->GetVid(0);
                     vId[cnt1+np-1]            = e->GetGeom()->GetVid(1);
                     vId[cnt1+newpoints-1]     = e->GetGeom()->GetVid(2);
 
                     // Edge numbering
                     StdRegions::Orientation             edgeOrient;
                     int edge1 = 0;
                     int edge2 = 0;
                     for (int n = 1; n < np-1; n++)
                     {
                         // Edge 0
                         edgeOrient          = e->GetGeom()->GetEorient(0);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             vId[cnt1+n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(0) + n;
                         }
                         else
                         {
                             vId[cnt1+np-1-n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(0) + n;
                         }
 
                         // Edge 1
                         edgeOrient          = e->GetGeom()->GetEorient(1);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             edge1 += np-n;
                             vId[cnt1+np-1+edge1] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(1) + n;
                         }
                         else
                         {
                             edge1 += n;
                             vId[cnt1+newpoints-1-edge1] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(1) + n;
                         }
 
                         // Edge 2
                         edgeOrient          = e->GetGeom()->GetEorient(2);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             edge2 += n+1;
                             vId[cnt1+newpoints-1-edge2] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(2) + n;
                         }
                         else
                         {
                             edge2 += np+1-n;
                             vId[cnt1+edge2] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(2) + n;
                         }
                     }
 
                     // Interior numbering
                     int mode = np+1;
                     for (int n = 1; n < np-1; n++)
                     {
                         for (int m = 1; m < np-n-1; m++)
                         {
                             vId[cnt1+mode] = 4*nel + maxN*4*nel + cnt2;
                             cnt2++;
                             mode++;
                         }
                     }
                     cnt1+= newpoints;
                 }
                 break;
             case LibUtilities::eQuadrilateral:
                 {
                     int newpoints = LibUtilities::StdQuadData::
                                         getNumberOfCoefficients(np,np);
                     // Vertex numbering
                     vId[cnt1]           = e->GetGeom()->GetVid(0);
                     vId[cnt1+np-1]      = e->GetGeom()->GetVid(1);
                     vId[cnt1+np*np-1]   = e->GetGeom()->GetVid(2);
                     vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);
 
                     // Edge numbering
                     StdRegions::Orientation             edgeOrient;
                     for (int n = 1; n < np-1; n++)
                     {
                         // Edge 0
                         edgeOrient          = e->GetGeom()->GetEorient(0);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             vId[cnt1+n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(0) + n;
                         }
                         else
                         {
                             vId[cnt1+np-1-n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(0) + n;
                         }
 
                         // Edge 1
                         edgeOrient          = e->GetGeom()->GetEorient(1);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             vId[cnt1+np-1+n*np] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(1) + n;
                         }
                         else
                         {
                             vId[cnt1+np*np-1-n*np] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(1) + n;
                         }
 
                         // Edge 2
                         edgeOrient          = e->GetGeom()->GetEorient(2);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             vId[cnt1+np*np-1-n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(2) + n;
                         }
                         else
                         {
                             vId[cnt1+np*(np-1)+n] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(2) + n;
                         }
 
                         // Edge 3
                         edgeOrient          = e->GetGeom()->GetEorient(3);
                         if (edgeOrient==StdRegions::eForwards)
                         {
                             vId[cnt1+np*(np-1)-n*np] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(3) + n;
                         }
                         else
                         {
                             vId[cnt1+n*np] = 4*nel +
                                             maxN*e->GetGeom()->GetEid(3) + n;
                         }
                     }
 
                     // Interior numbering
                     for (int n = 1; n < np-1; n++)
                     {
                         for (int m = 1; m < np-1; m++)
                         {
                             vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
                             cnt2++;
                         }
                     }
                     cnt1+= newpoints;
                 }
                 break;
             default:
                 {
                     ASSERTL0(false,"Points not known");
                 }
             }
         }
 
         // Crete new connectivity using homogeneous information
         for(int i = 0; i < nel; ++i)
         {
             int nTris = oldConn[i].num_elements()/3;
             // Create array for new connectivity
             // (3 tetrahedra between each plane for each triangle)
             conn = Array<OneD, int> (4*3*nTris*(nPlanes-1));
             cnt2=0;
             for(int n = 0; n<nTris; n++)
             {
                 // Get id of vertices in this triangle (on plane 0)
                 int vId0 = vId[oldConn[i][3*n+0]];
                 int vId1 = vId[oldConn[i][3*n+1]];
                 int vId2 = vId[oldConn[i][3*n+2]];
 
                 // Determine ordering based on global Id of pts
                 Array<OneD, int> rot(3);
                 if ( (vId0<vId1) && (vId0<vId2))
                 {
                     rot[0] = 0;
                     if (vId1<vId2)
                     {
                         rot[1] = 1;
                         rot[2] = 2;
                     }
                     else
                     {
                         rot[1] = 2;
                         rot[2] = 1;
                     }
                 }
                 else if ( (vId1<vId0) && (vId1<vId2))
                 {
                     rot[0] = 1;
                     if (vId0<vId2)
                     {
                         rot[1] = 0;
                         rot[2] = 2;
                     }
                     else
                     {
                         rot[1] = 2;
                         rot[2] = 0;
                     }
                 }
                 else
                 {
                     rot[0] = 2;
                     if (vId0<vId1)
                     {
                         rot[1] = 0;
                         rot[2] = 1;
                     }
                     else
                     {
                         rot[1] = 1;
                         rot[2] = 0;
                     }
                 }
 
                 // Define new connectivity with tetrahedra
                 for ( int p = 0; p<nPlanes-1; p++)
                 {
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[1]];
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[2]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
 
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[1]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
 
                     conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
                     conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
                 }
             }
             newConn.push_back(conn);
         }
     }
     else
     {
         ASSERTL0( false, "Points type not supported.");
     }
 
     m_f->m_fieldPts->SetConnectivity(newConn);
 }

void Nektar::FieldUtils::ProcessEquiSpacedOutput::SetupEquiSpacedField ( void )

protected

Definition at line 81 of file ProcessEquiSpacedOutput.cpp.

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

 {
     if (m_f->m_verbose)
     {
         if (m_f->m_comm->TreatAsRankZero())
         {
             cout << "Interpolating fields to equispaced..." << endl;
         }
     }
     int nel = m_f->m_exp[0]->GetExpSize();
     if (!nel)
     {
         m_f->m_fieldPts = LibUtilities::NullPtsField;
         return;
     }
 
     int coordim  = m_f->m_exp[0]->GetCoordim(0);
     int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
     int npts     = m_f->m_exp[0]->GetTotPoints();
     Array<OneD, Array<OneD, NekDouble> > coords(3);
 
     // Check if we have a homogeneous expansion
     bool homogeneous1D = false;
     if (m_f->m_fielddef.size())
     {
         if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
         {
             coordim++;
             shapedim++;
             homogeneous1D = true;
         }
         else if (m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
         {
             ASSERTL0(false, "Homegeneous2D case not supported");
         }
     }
     else
     {
         if (m_f->m_session->DefinesSolverInfo("HOMOGENEOUS"))
         {
             std::string HomoStr = m_f->m_session->GetSolverInfo("HOMOGENEOUS");
 
             if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
                 (HomoStr == "1D") || (HomoStr == "Homo1D"))
             {
                 coordim++;
                 shapedim++;
                 homogeneous1D = true;
             }
             if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
                 (HomoStr == "2D") || (HomoStr == "Homo2D"))
             {
                 ASSERTL0(false, "Homegeneous2D case not supported");
             }
         }
     }
 
     // set up the number of points in each element
     int newpoints    = 0;
     int newtotpoints = 0;
 
     Array<OneD, int> conn;
     int prevNcoeffs = 0;
     int prevNpoints = 0;
     int cnt         = 0;
 
     // identify face 1 connectivity for prisms
     map<int, StdRegions::Orientation> face0orient;
     set<int> prismorient;
     LocalRegions::ExpansionSharedPtr e;
 
     // prepare PtsField
     vector<std::string> fieldNames;
     vector<int> ppe;
     vector<Array<OneD, int> > ptsConn;
     int nfields;
 
     for (int i = 0; i < nel; ++i)
     {
         e = m_f->m_exp[0]->GetExp(i);
         if (e->DetShapeType() == LibUtilities::ePrism)
         {
             StdRegions::Orientation forient = e->GetForient(0);
             int fid                         = e->GetGeom()->GetFid(0);
             if (face0orient.count(fid))
             { // face 1 meeting face 1 so reverse this id
                 prismorient.insert(i);
             }
             else
             {
                 // just store if Dir 1 is fwd or bwd
                 if ((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
                     (forient == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
                     (forient == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
                     (forient == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
                 {
                     face0orient[fid] = StdRegions::eBwd;
                 }
                 else
                 {
                     face0orient[fid] = StdRegions::eFwd;
                 }
             }
         }
     }
 
     for (int i = 0; i < nel; ++i)
     {
         e = m_f->m_exp[0]->GetExp(i);
         if (e->DetShapeType() == LibUtilities::ePrism)
         {
             int fid = e->GetGeom()->GetFid(2);
             // check to see if face 2 meets face 1
             if (face0orient.count(fid))
             {
                 // check to see how face 2 is orientated
                 StdRegions::Orientation forient2 = e->GetForient(2);
                 StdRegions::Orientation forient0 = face0orient[fid];
 
                 // If dir 1 or forient2 is bwd then check agains
                 // face 1 value
                 if ((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
                     (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
                     (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
                     (forient2 == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
                 {
                     if (forient0 == StdRegions::eFwd)
                     {
                         prismorient.insert(i);
                     }
                 }
                 else
                 {
                     if (forient0 == StdRegions::eBwd)
                     {
                         prismorient.insert(i);
                     }
                 }
             }
         }
     }
 
     for (int i = 0; i < nel; ++i)
     {
         e = m_f->m_exp[0]->GetExp(i);
         if (m_config["tetonly"].m_beenSet)
         {
             if (m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
                 LibUtilities::eTetrahedron)
             {
                 continue;
             }
         }
 
         switch (e->DetShapeType())
         {
             case LibUtilities::eSegment:
             {
                 int npoints0 = e->GetBasis(0)->GetNumPoints();
 
                 newpoints =
                     LibUtilities::StdSegData::getNumberOfCoefficients(npoints0);
             }
             break;
             case LibUtilities::eTriangle:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np  = max(np0, np1);
                 newpoints =
                     LibUtilities::StdTriData::getNumberOfCoefficients(np, np);
             }
             break;
             case LibUtilities::eQuadrilateral:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np  = max(np0, np1);
 
                 newpoints =
                     LibUtilities::StdQuadData::getNumberOfCoefficients(np, np);
             }
             break;
             case LibUtilities::eTetrahedron:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np2 = e->GetBasis(2)->GetNumPoints();
                 int np  = max(np0, max(np1, np2));
 
                 newpoints = LibUtilities::StdTetData::getNumberOfCoefficients(
                     np, np, np);
             }
             break;
             case LibUtilities::ePrism:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np2 = e->GetBasis(2)->GetNumPoints();
                 int np  = max(np0, max(np1, np2));
 
                 newpoints = LibUtilities::StdPrismData::getNumberOfCoefficients(
                     np, np, np);
             }
             break;
             case LibUtilities::ePyramid:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np2 = e->GetBasis(2)->GetNumPoints();
                 int np  = max(np0, max(np1, np2));
 
                 newpoints = LibUtilities::StdPyrData::getNumberOfCoefficients(
                     np, np, np);
             }
             break;
             case LibUtilities::eHexahedron:
             {
                 int np0 = e->GetBasis(0)->GetNumPoints();
                 int np1 = e->GetBasis(1)->GetNumPoints();
                 int np2 = e->GetBasis(2)->GetNumPoints();
                 int np  = max(np0, max(np1, np2));
 
                 newpoints = LibUtilities::StdHexData::getNumberOfCoefficients(
                     np, np, np);
             }
             break;
             default:
             {
                 ASSERTL0(false, "Points not known");
             }
         }
 
         ppe.push_back(newpoints);
         newtotpoints += newpoints;
 
         if (e->DetShapeType() == LibUtilities::ePrism)
         {
             bool standard = true;
 
             if (prismorient.count(i))
             {
                 standard = false; // reverse direction
             }
 
             e->GetSimplexEquiSpacedConnectivity(conn, standard);
         }
         else
         {
 
             if ((prevNcoeffs != e->GetNcoeffs()) ||
                 (prevNpoints != e->GetTotPoints()))
             {
                 prevNcoeffs = e->GetNcoeffs();
                 prevNpoints = e->GetTotPoints();
 
                 e->GetSimplexEquiSpacedConnectivity(conn);
             }
         }
         Array<OneD, int> newconn(conn.num_elements());
         for (int j = 0; j < conn.num_elements(); ++j)
         {
             newconn[j] = conn[j] + cnt;
         }
 
         ptsConn.push_back(newconn);
         cnt += newpoints;
     }
 
     if (m_f->m_fielddef.size())
     {
         nfields = m_f->m_exp.size();
     }
     else // just the mesh points
     {
         nfields = 0;
     }
 
     Array<OneD, Array<OneD, NekDouble> > pts(nfields + coordim);
 
     for (int i = 0; i < nfields + coordim; ++i)
     {
         pts[i] = Array<OneD, NekDouble>(newtotpoints);
     }
 
     // Interpolate coordinates
     for (int i = 0; i < coordim; ++i)
     {
         coords[i] = Array<OneD, NekDouble>(npts);
     }
 
     for (int i = coordim; i < 3; ++i)
     {
         coords[i] = NullNekDouble1DArray;
     }
 
     m_f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
 
     int nq1 = m_f->m_exp[0]->GetTotPoints();
 
     Array<OneD, NekDouble> x1(nq1);
     Array<OneD, NekDouble> y1(nq1);
     Array<OneD, NekDouble> z1(nq1);
 
     m_f->m_exp[0]->GetCoords(x1, y1, z1);
 
     Array<OneD, NekDouble> tmp;
 
     for (int n = 0; n < coordim; ++n)
     {
         cnt      = 0;
         int cnt1 = 0;
         for (int i = 0; i < nel; ++i)
         {
             m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
                 coords[n] + cnt, tmp = pts[n] + cnt1);
             cnt1 += ppe[i];
             cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
         }
     }
 
     if (m_f->m_fielddef.size())
     {
         ASSERTL0(m_f->m_fielddef[0]->m_fields.size() == m_f->m_exp.size(),
                  "More expansion defined than fields");
 
         for (int n = 0; n < m_f->m_exp.size(); ++n)
         {
             cnt      = 0;
             int cnt1 = 0;
 
             if (m_config["modalenergy"].m_beenSet)
             {
                 Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
                 for (int i = 0; i < nel; ++i)
                 {
                     GenOrthoModes(i, phys + cnt, tmp = pts[coordim + n] + cnt1);
                     cnt1 += ppe[i];
                     cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
                 }
             }
             else
             {
                 Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
                 for (int i = 0; i < nel; ++i)
                 {
                     m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
                         phys + cnt, tmp = pts[coordim + n] + cnt1);
                     cnt1 += ppe[i];
                     cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
                 }
             }
 
             // Set up Variable string.
             fieldNames.push_back(m_f->m_fielddef[0]->m_fields[n]);
         }
     }
 
     m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
         coordim, fieldNames, pts);
     if (shapedim == 1)
     {
         m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsSegBlock);
     }
 
     if (shapedim == 2)
     {
         m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
     }
     else if (shapedim == 3)
     {
         m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
     }
     m_f->m_fieldPts->SetConnectivity(ptsConn);
     if (homogeneous1D)
     {
         SetHomogeneousConnectivity();
     }
 }

Member Data Documentation

ModuleKey Nektar::FieldUtils::ProcessEquiSpacedOutput::className

static

Initial value:

=
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "equispacedoutput"),
        ProcessEquiSpacedOutput::create,
        "Write data as equi-spaced output using simplices to represent the "
        "data for connecting points")

Definition at line 57 of file ProcessEquiSpacedOutput.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation