C++ SessionReaderSharedPtr::GetVariable方法代码示例

//.........这里部分代码省略.........
    vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
    vector<vector<NekDouble> >                      fieldData;

    LibUtilities::Import(fieldFile, fieldDef, fieldData);
    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    // Set up Expansion information
    vector< vector<LibUtilities::PointsType> > pointsType;
    for (i = 0; i < fieldDef.size(); ++i)
    {
        vector<LibUtilities::PointsType> ptype;
        for (j = 0; j < 3; ++j)
        {
            ptype.push_back(LibUtilities::ePolyEvenlySpaced);
        }
        pointsType.push_back(ptype);
    }
    graphShPt->SetExpansions(fieldDef, pointsType);

    //--------------------------------------------------------------------------


    //--------------------------------------------------------------------------
    // Define Expansion
    int nfields = fieldDef[0]->m_fields.size();
    Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
    Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);

    for(i = 0; i < pFields.num_elements(); i++)
    {
        pFields[i] = MemoryManager<MultiRegions
                                   ::DisContField3D>::AllocateSharedPtr(vSession, graphShPt,
                                                                        vSession->GetVariable(i));
    }

    MultiRegions::ExpList3DSharedPtr Exp3D;
    Exp3D = MemoryManager<MultiRegions::ExpList3D>
        ::AllocateSharedPtr(vSession, graphShPt);

    Exp[0] = Exp3D;

    for (i = 1; i < nfields; ++i)
    {
        Exp[i] = MemoryManager<MultiRegions::ExpList3D>
            ::AllocateSharedPtr(*Exp3D);
    }

    // Count of the point on the surface
    int nSurfacePts = 0;
    if (pFields[0]->GetBndCondExpansions().num_elements())
    {
        nSurfacePts = 0;
        cnt = 0;
        nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
        for (b = 0; b < nBndRegions; ++b)
        {
            nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
            for (e = 0; e < nBndEdges; ++e)
            {
                nBndEdgePts = pFields[0]->
                    GetBndCondExpansions()[b]->GetExp(e)->GetTotPoints();

                if (pFields[0]->GetBndConditions()[b]->
                        GetUserDefined() == "WallViscous" ||
                    pFields[0]->GetBndConditions()[b]->

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    MultiRegions::ContField1DSharedPtr Exp,Sol;

    int     i,j;
    int     order, nq;
    int     coordim;
    Array<OneD,NekDouble> sol;
    Array<OneD,NekDouble>  xc0,xc1,xc2;

    if(argc != 2)
    {
        fprintf(stderr,"Usage: ProjectCont1D mesh \n");
        exit(1);
    }

    //----------------------------------------------
    // read the problem parameters from input file
    SpatialDomains::MeshGraphSharedPtr graph1D = MemoryManager<SpatialDomains::MeshGraph1D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
    LibUtilities::BasisKey bkey0
                            = expansions.begin()->second->m_basisKeyVector[0];
    int nmodes = bkey0.GetNumModes(); 
    cout << "Solving 1D Continuous Projection"  << endl; 
    cout << "    Expansion  : (" << LibUtilities::BasisTypeMap[bkey0.GetBasisType()] <<")" << endl;
    cout << "    No. modes  : " << nmodes << endl;
    cout << endl;
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    Exp = MemoryManager<MultiRegions::ContField1D>
                ::AllocateSharedPtr(vSession,graph1D,vSession->GetVariable(0));
    //----------------------------------------------

    //----------------------------------------------
    // Define solution to be projected
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();
    order   = Exp->GetExp(0)->GetNcoeffs();

    // define coordinates and solution
    sol = Array<OneD,NekDouble>(nq);

    xc0 = Array<OneD,NekDouble>(nq);
    xc1 = Array<OneD,NekDouble>(nq);
    xc2 = Array<OneD,NekDouble>(nq);

    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        Vmath::Zero(nq,&xc1[0],1);
        Vmath::Zero(nq,&xc2[0],1);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        Vmath::Zero(nq,&xc2[0],1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }

    for(i = 0; i < nq; ++i)
    {
        sol[i] = 0.0;
        for(j = 0; j < order; ++j)
        {
            sol[i] += pow(xc0[i],j);
            sol[i] += pow(xc1[i],j);
            sol[i] += pow(xc2[i],j);
        }
    }
    //----------------------------------------------

    //----------------------------------------------
    // Setup Temporary expansion and put in solution
    Sol = MemoryManager<MultiRegions::ContField1D>
                                ::AllocateSharedPtr(*Exp);
    Sol->SetPhys(sol);
    //----------------------------------------------

    //---------------------------------------------
    // Project onto Expansion
    Exp->FwdTrans(Sol->GetPhys(), Exp->UpdateCoeffs());
    //---------------------------------------------

    //-------------------------------------------
    // Backward Transform Solution to get projected values
    Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys());
    //-------------------------------------------

//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    if(argc != 2)
    {
        fprintf(stderr,"Usage: ./Aliasing file.xml \n");
        fprintf(stderr,"\t Method will read intiial conditions section of .xml file for input \n");
        exit(1);
    }

    LibUtilities::SessionReaderSharedPtr session;
    string vDriverModule;
    DriverSharedPtr drv;  
    try
    {
        // Create session reader.
        session = LibUtilities::SessionReader::CreateInstance(argc, argv);
        
        // Create driver
        session->LoadSolverInfo("Driver", vDriverModule, "Standard");
        drv = GetDriverFactory().CreateInstance(vDriverModule, session);


        EquationSystemSharedPtr EqSys = drv->GetEqu()[0];
        IncNavierStokesSharedPtr IncNav = EqSys->as<IncNavierStokes>();
        
        IncNav->SetInitialConditions(0.0,false);
        Array<OneD, MultiRegions::ExpListSharedPtr> fields = IncNav->UpdateFields();
        
        int i;
        int nConvectiveFields = IncNav->GetNConvectiveFields();
        int nphys = fields[0]->GetTotPoints();
        Array<OneD, Array<OneD, NekDouble> > VelFields(nConvectiveFields);
        Array<OneD, Array<OneD, NekDouble> > NonLinear(nConvectiveFields);
        Array<OneD, Array<OneD, NekDouble> > NonLinearDealiased(nConvectiveFields);
        
        for(i = 0; i < nConvectiveFields; ++i)
        {
            VelFields[i] = fields[i]->UpdatePhys();
            NonLinear[i] = Array<OneD, NekDouble> (nphys);
            NonLinearDealiased[i] = Array<OneD, NekDouble> (nphys);
        }

        boost::shared_ptr<NavierStokesAdvection> A
            = boost::dynamic_pointer_cast<NavierStokesAdvection>(IncNav->GetAdvObject());

        if (!A)
        {
            cout << "Must use non-linear Navier-Stokes advection" << endl;
            exit(-1);
        }

        // calculate non-linear terms without dealiasing
        A->SetSpecHPDealiasing(false);
        A->Advect(nConvectiveFields, fields,
                                            VelFields, VelFields, 
                                            NonLinear, 0.0);


        // calculate non-linear terms with dealiasing
        A->SetSpecHPDealiasing(true);
        A->Advect(nConvectiveFields, fields,
                                            VelFields, VelFields, 
                                            NonLinearDealiased, 0.0);

        // Evaulate Difference and put into fields;
        for(i = 0; i < nConvectiveFields; ++i)
        {
            Vmath::Vsub(nphys,NonLinearDealiased[i],1,NonLinear[i],1,NonLinear[i],1);
            fields[i]->FwdTrans_IterPerExp(NonLinear[i],fields[i]->UpdateCoeffs());
            // Need to reset varibale name for output
            string name = "NL_Aliasing_"+session->GetVariable(i);
            session->SetVariable(i,name.c_str());
        }


        // Reset session name for output file
        std::string outname = IncNav->GetSessionName();
        
        outname += "_NonLinear_Aliasing";
        IncNav->ResetSessionName(outname);
        IncNav->Output();

    }
    catch (const std::runtime_error&)
    {
        return 1;
    }
    catch (const std::string& eStr)
    {
        cout << "Error: " << eStr << endl;
    }


    return 0;
}

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession = LibUtilities::SessionReader::CreateInstance(argc, argv);

    LibUtilities::CommSharedPtr vComm = vSession->GetComm();
    
    MultiRegions::ContField3DHomogeneous1DSharedPtr Exp_u, Exp_v, Exp_w;
    
    StdRegions::ConstFactorMap factors;
    FlagList flags;

    if( (argc != 2) && (argc != 3))
    {
        fprintf(stderr,"Usage: Deriv3DHomo2D meshfile [SysSolnType]   \n");
        exit(1);
    }

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    int nzpoints;
    NekDouble lz;
    int FFT;

    vSession->LoadParameter("HomModesZ", nzpoints);
	vSession->LoadParameter("LZ",        lz);
	vSession->LoadParameter("USEFFT",    FFT);
	
	bool useFFT = false;
	bool deal = false;
	if(FFT==1){useFFT = true;}
			
	
	const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierSingleModeSpaced);
	const LibUtilities::BasisKey  BkeyZ(LibUtilities::eFourierHalfModeRe,nzpoints,PkeyZ);
		
	Exp_u = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(0));
	Exp_v = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(1));
	Exp_w = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(2));
		

    //----------------------------------------------
    // Print summary of solution details
        flags.set(eUseGlobal, false);
		
    const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
	
    LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
    
	cout << "Calculating Derivatives (Homogeneous in z-plane):"  << endl;
	cout << "         Lz              : " << lz << endl;
    cout << "         N.modes         : " << bkey0.GetNumModes() << endl;
	cout << "         N.Z homo modes  : " << BkeyZ.GetNumModes() << endl;
    cout << endl;
    //----------------------------------------------

    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
	int nq  = Exp_u->GetTotPoints();
	
	Array<OneD,NekDouble>  xc0,xc1,xc2;
    
    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);

    Exp_u->GetCoords(xc0,xc1,xc2);
    //----------------------------------------------
    Array<OneD,NekDouble>  dudx,dvdy,dwdz;
	Array<OneD,NekDouble>  dump;
	dump = Array<OneD,NekDouble>(nq,0.0);
	dudx = Array<OneD,NekDouble>(nq,0.0);
	dvdy = Array<OneD,NekDouble>(nq,0.0);
	dwdz = Array<OneD,NekDouble>(nq,0.0);
    //----------------------------------------------
    // Define initial fields
    LibUtilities::EquationSharedPtr ffunc_u = vSession->GetFunction("InitialCondition", 0);
	LibUtilities::EquationSharedPtr ffunc_v = vSession->GetFunction("InitialCondition", 1);
	LibUtilities::EquationSharedPtr ffunc_w = vSession->GetFunction("InitialCondition", 2);
	
	LibUtilities::EquationSharedPtr exac_u = vSession->GetFunction("ExactSolution", 0);
	LibUtilities::EquationSharedPtr exac_v = vSession->GetFunction("ExactSolution", 1);
	LibUtilities::EquationSharedPtr exac_w = vSession->GetFunction("ExactSolution", 2);
    	

    ffunc_u->Evaluate(xc0,xc1,xc2,Exp_u->UpdatePhys());
    ffunc_v->Evaluate(xc0,xc1,xc2,Exp_v->UpdatePhys());
    ffunc_w->Evaluate(xc0,xc1,xc2,Exp_w->UpdatePhys());

    exac_u->Evaluate(xc0,xc1,xc2,dudx);
    exac_v->Evaluate(xc0,xc1,xc2,dvdy);
    exac_w->Evaluate(xc0,xc1,xc2,dwdz);

    //----------------------------------------------
	
    //Taking derivative and printing the error
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    MultiRegions::ContField2DSharedPtr Exp,Fce;
    int     nq,  coordim;
    Array<OneD,NekDouble>  fce; 
    Array<OneD,NekDouble>  xc0,xc1,xc2; 
    NekDouble  lambda;
    NekDouble  ax,ay;

    if((argc != 2)&&(argc != 3))
    {
        fprintf(stderr,"Usage: SteadyLinearAdvectionReaction2D  meshfile [SysSolnType]\n");
        exit(1);
    }

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Get Advection Velocity
    ax = vSession->GetParameter("Advection_x");
    ay = vSession->GetParameter("Advection_y");
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    lambda = vSession->GetParameter("Lambda");
    cout << "            Lambda         : " << lambda << endl;
    const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
    LibUtilities::BasisKey bkey0
                            = expansions.begin()->second->m_basisKeyVector[0];
    LibUtilities::BasisKey bkey1
                            = expansions.begin()->second->m_basisKeyVector[1];
    cout << "Solving Steady 2D LinearAdvection :"  << endl; 
    cout << "            Advection_x    : " << ax << endl; 
    cout << "            Advection_y    : " << ay << endl; 
    cout << "            Expansion      : (" << LibUtilities::BasisTypeMap[bkey0.GetBasisType()] <<","<< LibUtilities::BasisTypeMap[bkey1.GetBasisType()]  << ")" << endl;
    cout << "            No. modes      : " << bkey0.GetNumModes() << endl;
    cout << endl;
    //----------------------------------------------
   
    //----------------------------------------------
    // Define Expansion 
    Exp = MemoryManager<MultiRegions::ContField2D>::
        AllocateSharedPtr(vSession,graph2D,vSession->GetVariable(0));
    //----------------------------------------------

    Timing("Read files and define exp ..");
    
    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();
    
    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);
    
    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }

    Array<OneD, Array< OneD, NekDouble> > Vel(2);
    Vel[0] = Array<OneD, NekDouble> (nq,ax);
    Vel[1] = Array<OneD, NekDouble> (nq,ay);
    //----------------------------------------------
    
    //----------------------------------------------
    // Define forcing function for first variable defined in file 
    fce = Array<OneD,NekDouble>(nq);
    LibUtilities::EquationSharedPtr ffunc = vSession->GetFunction("Forcing",0);

    ffunc->Evaluate(xc0,xc1,xc2,fce);
    //----------------------------------------------

    //----------------------------------------------
    // Setup expansion containing the  forcing function
    Fce = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(*Exp);
    Fce->SetPhys(fce);
    //----------------------------------------------
    Timing("Define forcing ..");
  
    //----------------------------------------------
    // Helmholtz solution taking physical forcing 
    Exp->LinearAdvectionReactionSolve(Vel, Fce->GetPhys(), Exp->UpdateCoeffs(), lambda, MultiRegions::eGlobal);
    //----------------------------------------------
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    int i,j;
    int surfID;

    if(argc != 5)
    {
        fprintf(stderr,"Usage: FldAddScalGrad meshfile infld outfld BoundaryID\n");
        exit(1);
    }

    surfID = boost::lexical_cast<int>(argv[argc - 1]);

    argv[argc -1] = argv[argc - 2];

    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    //----------------------------------------------
    // Read in mesh from input file
    string meshfile(argv[argc-4]);
    SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Import field file.
    string fieldfile(argv[argc-3]);
    vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
    vector<vector<NekDouble> > fielddata;
    LibUtilities::Import(fieldfile,fielddef,fielddata);
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    int expdim  = graphShPt->GetMeshDimension();
    int nfields = 1;
    int addfields = 7;
    Array<OneD, MultiRegions::ExpListSharedPtr> exp(nfields + addfields);
    MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;

    switch(expdim)
    {
        case 1:
        {
            ASSERTL0(false,"Expansion dimension not recognised");
        }
        break;
        case 2:
        {
            ASSERTL0(false,"Expansion dimension not recognised");
        }
        break;
        case 3:
        {
            MultiRegions::ContField3DSharedPtr originalfield =
                MemoryManager<MultiRegions::ContField3D>
                ::AllocateSharedPtr(vSession, graphShPt, 
                                    vSession->GetVariable(0));

            m_locToGlobalMap = originalfield->GetLocalToGlobalMap();
            
            exp[0] = originalfield;
            for (i=0; i<addfields; i++)
            {
                exp[i+1] = MemoryManager<MultiRegions::ContField3D>
                    ::AllocateSharedPtr(*originalfield, graphShPt, 
                                        vSession->GetVariable(0));
            }
        }
        break;
        default:
            ASSERTL0(false,"Expansion dimension not recognised");
            break;
    }
    //----------------------------------------------

    //----------------------------------------------
    // Copy data from field file
    for(j = 0; j < nfields+addfields; ++j)
    {
        for(int i = 0; i < fielddata.size(); ++i)
        {
            exp[j]->ExtractDataToCoeffs(fielddef [i],
                                       fielddata[i],
                                        fielddef [i]->m_fields[0],
                                        exp[j]->UpdateCoeffs());
        }
        
        exp[j]->BwdTrans(exp[j]->GetCoeffs(),exp[j]->UpdatePhys());
    }


    //----------------------------------------------

    //----------------------------------------------
    int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
    nt = exp[0]->GetNpoints();
    Array<OneD, Array<OneD, NekDouble> > grad(expdim);
    Array<OneD, Array<OneD, NekDouble> > fgrad(expdim);
    Array<OneD, Array<OneD, NekDouble> > values(addfields);
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
        = LibUtilities::SessionReader::CreateInstance(argc, argv);

    LibUtilities::CommSharedPtr vComm = vSession->GetComm();
    string meshfile(argv[1]);

    MultiRegions::DisContField3DHomogeneous1DSharedPtr Exp,Fce;
    MultiRegions::ExpListSharedPtr DerExp1,DerExp2,DerExp3;
    int i, nq;
    Array<OneD,NekDouble>  fce;
    Array<OneD,NekDouble>  xc0,xc1,xc2;
    StdRegions::ConstFactorMap factors;
    NekDouble lz;

    if(argc != 2)
    {
        fprintf(stderr,"Usage: Helmholtz2D  meshfile\n");
        exit(1);
    }

    LibUtilities::FieldIOSharedPtr fld = MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(vComm);

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    int nplanes      = vSession->GetParameter("HomModesZ");
    lz     = vSession->GetParameter("LZ");
    bool useFFT = false;
    bool deal = false;
    const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::eFourierEvenlySpaced);
    const LibUtilities::BasisKey Bkey(LibUtilities::eFourier,nplanes,Pkey);
    Exp = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>::
          AllocateSharedPtr(vSession,Bkey,lz,useFFT,deal,graph2D,vSession->GetVariable(0));
    //----------------------------------------------
    Timing("Read files and define exp ..");



    //----------------------------------------------
    // Print summary of solution details
    factors[StdRegions::eFactorLambda]  = vSession->GetParameter("Lambda");
    factors[StdRegions::eFactorTau] = 1.0;

    const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
    LibUtilities::BasisKey bkey0
        = expansions.begin()->second->m_basisKeyVector[0];
    cout << "Solving 3D Helmholtz (Homogeneous in z-direction):"  << endl;
    cout << "         Lambda         : " << factors[StdRegions::eFactorLambda] << endl;
    cout << "         Lz             : " << lz << endl;
    cout << "         No. modes      : " << bkey0.GetNumModes() << endl;
    cout << "         No. hom. modes : " << Bkey.GetNumModes() << endl;
    cout << endl;
    //----------------------------------------------

    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    nq  = Exp->GetTotPoints();
    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);

    Exp->GetCoords(xc0,xc1,xc2);
    //----------------------------------------------

    //----------------------------------------------
    // Define forcing function for first variable defined in file
    fce = Array<OneD,NekDouble>(nq);
    LibUtilities::EquationSharedPtr ffunc
        = vSession->GetFunction("Forcing", 0);

    ffunc->Evaluate(xc0, xc1, xc2, fce);

    //----------------------------------------------


    //----------------------------------------------
    // Setup expansion containing the  forcing function
    Fce = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>::AllocateSharedPtr(*Exp);
    Fce->SetPhys(fce);
    //----------------------------------------------
    Timing("Define forcing ..");

    //----------------------------------------------
    // Helmholtz solution taking physical forcing
    Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
    //----------------------------------------------

    Timing("Helmholtz Solve ..");

#ifdef TIMING
    for(i = 0; i < 100; ++i)
    {
        Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
    }
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    int cnt;
    int id1, id2;
    int i, j, e, b;
    
    int nBndEdgePts, nBndEdges, nBndRegions;
        
    if (argc < 3)
    {
        fprintf(stderr,
                "Usage: ExtractSurface2DCFS meshfile fieldFile\n");
        fprintf(stderr,
                "Extracts a surface from a 2D fld file" 
                "(only for CompressibleFlowSolver and purely 2D .fld files)\n");
        exit(1);
    }

    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(3, argv);

    //--------------------------------------------------------------------------
    // Read in mesh from input file
    string meshfile(argv[1]);
    SpatialDomains::MeshGraphSharedPtr graphShPt = 
        SpatialDomains::MeshGraph::Read(vSession);
    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    // Import field file
    string                                          fieldFile(argv[2]);
    vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
    vector<vector<NekDouble> >                      fieldData;
    
    LibUtilities::Import(fieldFile, fieldDef, fieldData);
    //--------------------------------------------------------------------------

    //--------------------------------------------------------------------------
    // Set up Expansion information
    vector< vector<LibUtilities::PointsType> > pointsType;
    for (i = 0; i < fieldDef.size(); ++i)
    {
        vector<LibUtilities::PointsType> ptype;
        for (j = 0; j < 2; ++j)
        {
            ptype.push_back(LibUtilities::ePolyEvenlySpaced);
        }
        pointsType.push_back(ptype);
    }
    graphShPt->SetExpansions(fieldDef, pointsType);

    //--------------------------------------------------------------------------


    //--------------------------------------------------------------------------
    // Define Expansion
    int nfields = fieldDef[0]->m_fields.size();
    Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
    Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);
        
    for(i = 0; i < pFields.num_elements(); i++)
    {
        pFields[i] = MemoryManager<MultiRegions
        ::DisContField2D>::AllocateSharedPtr(vSession, graphShPt, 
                                             vSession->GetVariable(i));
    }
    
    MultiRegions::ExpList2DSharedPtr Exp2D;
    Exp2D = MemoryManager<MultiRegions::ExpList2D>
        ::AllocateSharedPtr(vSession, graphShPt);
    
    Exp[0] = Exp2D;

    for (i = 1; i < nfields; ++i)
    {
        Exp[i] = MemoryManager<MultiRegions::ExpList2D>
            ::AllocateSharedPtr(*Exp2D);
    }
    
    int nSolutionPts = pFields[0]->GetNpoints();
    int nTracePts    = pFields[0]->GetTrace()->GetTotPoints();

    Array<OneD, NekDouble> x(nSolutionPts);
    Array<OneD, NekDouble> y(nSolutionPts); 
    Array<OneD, NekDouble> z(nSolutionPts);
    
    Array<OneD, NekDouble> traceX(nTracePts);
    Array<OneD, NekDouble> traceY(nTracePts); 
    Array<OneD, NekDouble> traceZ(nTracePts);
    
    Array<OneD, NekDouble> surfaceX(nTracePts);
    Array<OneD, NekDouble> surfaceY(nTracePts); 
    Array<OneD, NekDouble> surfaceZ(nTracePts);

    pFields[0]->GetCoords(x, y, z);
    
    pFields[0]->ExtractTracePhys(x, traceX);
    pFields[0]->ExtractTracePhys(y, traceY);
    pFields[0]->ExtractTracePhys(z, traceZ);
    //--------------------------------------------------------------------------
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    MultiRegions::ContField3DSharedPtr Exp,Fce;
    int     i, j, nq,  coordim;
    Array<OneD,NekDouble>  fce;
    Array<OneD,NekDouble>  xc0,xc1,xc2;

    if(argc != 2)
    {
        fprintf(stderr,"Usage: ProjectCont3D  meshfile \n");
        exit(1);
    }

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph3D = MemoryManager<SpatialDomains::MeshGraph3D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
    LibUtilities::BasisKey bkey
                            = expansions.begin()->second->m_basisKeyVector[0];
    int nmodes =  bkey.GetNumModes();
    cout << "Solving 3D C0 continuous Projection"  << endl;
    cout << "    No. modes  : " << nmodes << endl;
    cout << endl;
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    Exp = MemoryManager<MultiRegions::ContField3D>
                                ::AllocateSharedPtr(vSession,graph3D,vSession->GetVariable(0));
    //----------------------------------------------

    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();

    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);

    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }
    //----------------------------------------------

    //----------------------------------------------
    // Define forcing function
    fce = Array<OneD,NekDouble>(nq);

    for(i = 0; i < nq; ++i)
    {
        fce[i] = 0.0;
        for(j = 0; j < nmodes; ++j)
        {
            fce[i] += pow(xc0[i],j);
            fce[i] += pow(xc1[i],j);
            fce[i] += pow(xc2[i],j);
        }
    }

    //---------------------------------------------
    // Set up ExpList1D containing the solution
    Fce = MemoryManager<MultiRegions::ContField3D>::AllocateSharedPtr(*Exp);
    Fce->SetPhys(fce);
    //---------------------------------------------

    //----------------------------------------------
    // Write solution
    //ofstream outfile("ProjectContFileOrig3D.dat");
    //Fce->WriteToFile(outfile,eGnuplot);
    //outfile.close();
    //----------------------------------------------

    //---------------------------------------------
    // Project onto Expansion
    Exp->FwdTrans(Fce->GetPhys(), Exp->UpdateCoeffs(), MultiRegions::eGlobal);    
    //---------------------------------------------

    //-------------------------------------------
    // Backward Transform Solution to get projected values
    Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
    //-------------------------------------------

    //----------------------------------------------
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    LibUtilities::CommSharedPtr vComm = vSession->GetComm();

    MultiRegions::DisContField3DSharedPtr Exp,Fce;
    int     i, nq, coordim;
    Array<OneD,NekDouble>  fce;
    Array<OneD,NekDouble>  xc0,xc1,xc2;
    StdRegions::ConstFactorMap factors;

    if(argc < 2)
    {
        fprintf(stderr,"Usage: PostProcHDG3D  meshfile [solntype]\n");
        exit(1);
    }

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph3D = MemoryManager<SpatialDomains::MeshGraph3D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
    factors[StdRegions::eFactorTau] = 1.0;
    const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
    LibUtilities::BasisKey bkey0
                            = expansions.begin()->second->m_basisKeyVector[0];

	//MAY NEED ADJUSTMENT FOR VARIOUS ELEMENT TYPES
	int num_modes = bkey0.GetNumModes();
	int num_points = bkey0.GetNumPoints();

    if (vComm->GetRank() == 0)
    {
        cout << "Solving 3D Helmholtz:"  << endl;
        cout << "         Lambda     : " << factors[StdRegions::eFactorLambda] << endl;
        cout << "         No. modes  : " << num_modes << endl;
        cout << "         No. points : " << num_points << endl;
        cout << endl;
    }

    //----------------------------------------------
    // Define Expansion
    //----------------------------------------------
    Exp = MemoryManager<MultiRegions::DisContField3D>::
        AllocateSharedPtr(vSession,graph3D,vSession->GetVariable(0));
    //----------------------------------------------
    Timing("Read files and define exp ..");

    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();

    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);

    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }
    //----------------------------------------------

    //----------------------------------------------
    // Define forcing function for first variable defined in file
    fce = Array<OneD,NekDouble>(nq);
    LibUtilities::EquationSharedPtr ffunc
                                    = vSession->GetFunction("Forcing", 0);

    ffunc->Evaluate(xc0, xc1, xc2, fce);

    //----------------------------------------------


    //----------------------------------------------
    // Setup expansion containing the  forcing function
    Fce = MemoryManager<MultiRegions::DisContField3D>::AllocateSharedPtr(*Exp);
    Fce->SetPhys(fce);
    //----------------------------------------------
    Timing("Define forcing ..");

    //----------------------------------------------
    // Helmholtz solution taking physical forcing
    Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
    //----------------------------------------------

    Timing("Helmholtz Solve ..");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
        case 0:
            type = StdRegions::eBwdTrans;
            break;
        case 1:
            type = StdRegions::eIProductWRTBase;
            break;
        case 2:
            type = StdRegions::eMass;
            break;
        case 3:
            type = StdRegions::eHelmholtz;
            break;
        default:
            cout << "Operator " << opToTest << " not defined." << endl;
    }

    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv, vFilenames);

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph3D = MemoryManager<SpatialDomains::MeshGraph3D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    lambda = vSession->GetParameter("Lambda");
    //----------------------------------------------

    //----------------------------------------------
    // Define Expansion
    Exp = MemoryManager<MultiRegions::ContField3D>
                    ::AllocateSharedPtr(vSession, graph3D, vSession->GetVariable(0));
    //----------------------------------------------
    int NumElements = Exp->GetExpSize();

    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();

    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);

    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }
    //----------------------------------------------

    //----------------------------------------------
    // Define forcing function for first variable defined in file
    fce = Array<OneD,NekDouble>(nq);
    LibUtilities::EquationSharedPtr ffunc = vSession->GetFunction("Forcing",0);
    for(i = 0; i < nq; ++i)
    {

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    LibUtilities::CommSharedPtr vComm = vSession->GetComm();

    MultiRegions::DisContField3DSharedPtr Exp, Fce;
    int     i, nq,  coordim;
    Array<OneD,NekDouble>  fce; 
    Array<OneD,NekDouble>  xc0,xc1,xc2; 
    StdRegions::ConstFactorMap factors;

    if(argc < 2)
    {
        fprintf(stderr,"Usage: HDGHelmholtz3D  meshfile [solntype]\n");
        exit(1);
    }

    LibUtilities::FieldIOSharedPtr fld = MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(vComm);

    //----------------------------------------------
    // Read in mesh from input file
    SpatialDomains::MeshGraphSharedPtr graph3D = 
        MemoryManager<SpatialDomains::MeshGraph3D>::AllocateSharedPtr(vSession);
    //----------------------------------------------

    //----------------------------------------------
    // Print summary of solution details
    factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
    factors[StdRegions::eFactorTau] = 1.0;
    const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
    LibUtilities::BasisKey bkey0
                            = expansions.begin()->second->m_basisKeyVector[0];

    if (vComm->GetRank() == 0)
    {
            cout << "Solving 3D Helmholtz:"  << endl;
            cout << "  - Communication: " 
                 << vSession->GetComm()->GetType() << " (" 
                 << vSession->GetComm()->GetSize() 
                 << " processes)" << endl;
            cout << "  - Solver type  : " 
                 << vSession->GetSolverInfo("GlobalSysSoln") << endl;
            cout << "  - Lambda       : " 
                 << factors[StdRegions::eFactorLambda] << endl;
            cout << "  - No. modes    : " 
                 << bkey0.GetNumModes() << endl;
            cout << endl;
    }
    //----------------------------------------------
   
    //----------------------------------------------
    // Define Expansion 
    Exp = MemoryManager<MultiRegions::DisContField3D>::
        AllocateSharedPtr(vSession,graph3D,vSession->GetVariable(0));
    //----------------------------------------------
    Timing("Read files and define exp ..");
    
    //----------------------------------------------
    // Set up coordinates of mesh for Forcing function evaluation
    coordim = Exp->GetCoordim(0);
    nq      = Exp->GetTotPoints();
    
    xc0 = Array<OneD,NekDouble>(nq,0.0);
    xc1 = Array<OneD,NekDouble>(nq,0.0);
    xc2 = Array<OneD,NekDouble>(nq,0.0);
    
    switch(coordim)
    {
    case 1:
        Exp->GetCoords(xc0);
        break;
    case 2:
        Exp->GetCoords(xc0,xc1);
        break;
    case 3:
        Exp->GetCoords(xc0,xc1,xc2);
        break;
    }
    //----------------------------------------------
    
    //----------------------------------------------
    // Define forcing function for first variable defined in file 
    fce = Array<OneD,NekDouble>(nq);
    LibUtilities::EquationSharedPtr ffunc = vSession->GetFunction("Forcing", 0);

    ffunc->Evaluate(xc0, xc1, xc2, fce);

    //----------------------------------------------


    //----------------------------------------------
    // Setup expansion containing the  forcing function
    Fce = MemoryManager<MultiRegions::DisContField3D>::AllocateSharedPtr(*Exp);
    Fce->SetPhys(fce);
    //----------------------------------------------
    Timing("Define forcing ..");
  
    //----------------------------------------------
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
    LibUtilities::SessionReaderSharedPtr vSession
            = LibUtilities::SessionReader::CreateInstance(argc, argv);

    LibUtilities::CommSharedPtr vComm = vSession->GetComm();
    MultiRegions::ContField1DSharedPtr Exp,Fce;
    int     i, nq,  coordim;
    Array<OneD,NekDouble>  fce;
    Array<OneD,NekDouble>  xc0,xc1,xc2;
    StdRegions::ConstFactorMap factors;

    if( (argc != 2) && (argc != 3) && (argc != 4))
    {
        fprintf(stderr,"Usage: Helmholtz1D  meshfile \n");
        exit(1);
    }

    try
    {
        LibUtilities::FieldIOSharedPtr fld =
            MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(vComm);

        //----------------------------------------------
        // Read in mesh from input file
        SpatialDomains::MeshGraphSharedPtr graph1D =
            SpatialDomains::MeshGraph::Read(vSession);
        //----------------------------------------------

        //----------------------------------------------
        // Print summary of solution details
        factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
        const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
        LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];

        if (vComm->GetRank() ==0)
        {
            cout << "Solving 1D Helmholtz: "  << endl;
            cout << "       Communication: " << vComm->GetType() << endl;
            cout << "       Solver type  : " << vSession->GetSolverInfo("GlobalSysSoln") << endl;
            cout << "       Lambda       : " << factors[StdRegions::eFactorLambda] << endl;
            cout << "       No. modes    : " << bkey0.GetNumModes() << endl;
        }
        //----------------------------------------------

        //----------------------------------------------
        // Define Expansion
        Exp = MemoryManager<MultiRegions::ContField1D>::
            AllocateSharedPtr(vSession,graph1D,vSession->GetVariable(0));
        //----------------------------------------------

        //----------------------------------------------
        // Set up coordinates of mesh for Forcing function evaluation
        coordim = Exp->GetCoordim(0);
        nq      = Exp->GetTotPoints();

        xc0 = Array<OneD,NekDouble>(nq);
        xc1 = Array<OneD,NekDouble>(nq);
        xc2 = Array<OneD,NekDouble>(nq);

        switch(coordim)
        {
        case 1:
            Exp->GetCoords(xc0);
            Vmath::Zero(nq,&xc1[0],1);
            Vmath::Zero(nq,&xc2[0],1);
            break;
        case 2:
            Exp->GetCoords(xc0,xc1);
            Vmath::Zero(nq,&xc2[0],1);
            break;
        case 3:
            Exp->GetCoords(xc0,xc1,xc2);
            break;
        }
        //----------------------------------------------

        //----------------------------------------------
        // Define forcing function for first variable defined in file
        fce = Array<OneD,NekDouble>(nq);
        LibUtilities::EquationSharedPtr ffunc
                                        = vSession->GetFunction("Forcing", 0);

        ffunc->Evaluate(xc0,xc1,xc2, fce);

        //----------------------------------------------

        //----------------------------------------------
        // Setup expansion containing the  forcing function
        Fce = MemoryManager<MultiRegions::ContField1D>::AllocateSharedPtr(*Exp);
        Fce->SetPhys(fce);
        //----------------------------------------------

        //----------------------------------------------
        //Helmholtz solution taking physical forcing after setting
        //initial condition to zero
        Vmath::Zero(Exp->GetNcoeffs(),Exp->UpdateCoeffs(),1);
        Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
        //----------------------------------------------

//.........这里部分代码省略.........