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

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::Utilities::ProcessEquiSpacedOutput:

Collaboration diagram for Nektar::Utilities::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::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< 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::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, ConfigOption >	m_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.

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

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

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

     {
         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::Utilities::ProcessEquiSpacedOutput::GetModuleName ( )

inlinevirtual

Implements Nektar::Utilities::Module.

Reimplemented in Nektar::Utilities::ProcessIsoContour.

Definition at line 65 of file ProcessEquiSpacedOutput.h.

         {
             return "ProcessEquiSpacedOutput";
         }

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

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

 {
     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;
 
     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 == 2)
                 {
                     // 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 == 2 && (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;
     }
 
     // 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
     vector<Array<OneD, int> > oldConn;
     vector<Array<OneD, int> > newConn;
     Array<OneD, int>          conn;
     m_f->m_fieldPts->GetConnectivity(oldConn);
 
     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);
     }
     m_f->m_fieldPts->SetConnectivity(newConn);
 }

void Nektar::Utilities::ProcessEquiSpacedOutput::SetupEquiSpacedField ( void )

protected

Definition at line 82 of file ProcessEquiSpacedOutput.cpp.

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

 {
     if(m_f->m_verbose)
     {
         if(m_f->m_comm->TreatAsRankZero())
         {
             cout << "Interpolating fields to equispaced..." << endl;
         }
     }
 
     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);
 
     int nel = m_f->m_exp[0]->GetExpSize();
 
     // 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)
         {
             ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
             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"))
             {
                 ASSERTL0(shapedim==2, "Homogeneous only implemented for 3DH1D");
                 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 == 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::Utilities::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