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

int main(int argc, char *argv[])
{
    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);

        // Execute driver
        drv->Execute();

        // Finalise session
        session->Finalise();
    }
    catch (const std::runtime_error& e)
    {
        return 1;
    }
    catch (const std::string& eStr)
    {
        cout << "Error: " << eStr << endl;
    }

    return 0;
}

//.........这里部分代码省略.........
    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());
    //-------------------------------------------

    //--------------------------------------------
    // Calculate L_inf error
    cout << "L infinity error: " << Exp->Linf(Sol->GetPhys()) << endl;
    cout << "L 2 error:        " << Exp->L2  (Sol->GetPhys()) << endl;
    //--------------------------------------------

    vSession->Finalise();

    return 0;
}

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

    MultiRegions::ExpList1DSharedPtr Exp,Sol;
    int i,j;
    int nq;
    int coordim;
    Array<OneD, NekDouble> sol;
    Array<OneD, NekDouble> xc0,xc1,xc2;

    // read in mesh
    SpatialDomains::MeshGraphSharedPtr graph1D = MemoryManager<SpatialDomains::MeshGraph1D>::AllocateSharedPtr(vSession);

    // Define Expansion
    const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
    LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
    int nmodes = bkey0.GetNumModes();

    Exp = MemoryManager<MultiRegions::ExpList1D>::AllocateSharedPtr(vSession,bkey0,graph1D);

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

    // 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 < nmodes; ++j)
        {
            sol[i] += pow(xc0[i],j);
            sol[i] += pow(xc1[i],j);
            sol[i] += pow(xc2[i],j);
        }
    }

    //---------------------------------------------
    // Set up ExpList1D containing the solution
    Sol = MemoryManager<MultiRegions::ExpList1D>::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());
    //-------------------------------------------

    //--------------------------------------------
    // Calculate L_inf error
    if (vSession->GetComm()->GetRank() == 0)
    {
        cout << "L infinity error: " << Exp->Linf(Sol->GetPhys()) << endl;
        cout << "L 2 error:        " << Exp->L2  (Sol->GetPhys()) << endl;
    }
    //--------------------------------------------

    vSession->Finalise();

    return 0;
}

//.........这里部分代码省略.........
                            = 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);
    //----------------------------------------------
    Timing("Linear Advection Solve ..");
    
    //----------------------------------------------
    // Backward Transform Solution to get solved values 
    Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
    //----------------------------------------------
    
    //----------------------------------------------
    // See if there is an exact solution, if so 
    // evaluate and plot errors
    LibUtilities::EquationSharedPtr ex_sol = vSession->GetFunction("ExactSolution",0);

    if(ex_sol)
    {
        //----------------------------------------------
        // evaluate exact solution 
        ex_sol->Evaluate(xc0,xc1,xc2,fce);
        //----------------------------------------------

        //--------------------------------------------
        // Calculate L_inf error 
        Fce->SetPhys(fce);
        Fce->SetPhysState(true);


        cout << "L infinity error: " << Exp->Linf(Fce->GetPhys()) << endl;
        cout << "L 2 error:        " << Exp->L2  (Fce->GetPhys()) << endl;
        //--------------------------------------------        
    }
    //----------------------------------------------        

    vSession->Finalise();

    return 0;
}

//.........这里部分代码省略.........
    //----------------------------------------------
    // 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);
    //-------------------------------------------

    //----------------------------------------------
    // Write solution
    //ofstream outfile2("ProjectContFile3D.dat");
    //Exp->WriteToFile(outfile2,eGnuplot);
    //outfile2.close();
    //----------------------------------------------

    //--------------------------------------------
    // Calculate L_inf error
    cout << "L infinity error: " << Exp->Linf(Fce->GetPhys()) << endl;
    cout << "L 2 error:        " << Exp->L2  (Fce->GetPhys()) << endl;
    //--------------------------------------------

    vSession->Finalise();

    return 0;
}

//.........这里部分代码省略.........
	const LibUtilities::PointsKey PkeyT3(num_points,LibUtilities::eGaussRadauMAlpha2Beta0);//need to doublecheck this one
	LibUtilities::BasisKeyVector  BkeyT;
	BkeyT.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyT1));
	BkeyT.push_back(LibUtilities::BasisKey(LibUtilities::eModified_B, num_modes+1, PkeyT2));
	BkeyT.push_back(LibUtilities::BasisKey(LibUtilities::eModified_C, num_modes+1, PkeyT3));
	//Prism
	const LibUtilities::PointsKey PkeyP1(num_points+1,LibUtilities::eGaussLobattoLegendre);
	const LibUtilities::PointsKey PkeyP2(num_points+1,LibUtilities::eGaussLobattoLegendre);
	const LibUtilities::PointsKey PkeyP3(num_points,LibUtilities::eGaussRadauMAlpha1Beta0);//need to doublecheck this one
	LibUtilities::BasisKeyVector  BkeyP;
	BkeyP.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyP1));
	BkeyP.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyP2));
	BkeyP.push_back(LibUtilities::BasisKey(LibUtilities::eModified_B, num_modes+1, PkeyP3));
	//Hexahedron
	const LibUtilities::PointsKey PkeyH(num_points+1,LibUtilities::eGaussLobattoLegendre);
	LibUtilities::BasisKeyVector  BkeyH;
	BkeyH.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyH));
	BkeyH.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyH));
	BkeyH.push_back(LibUtilities::BasisKey(LibUtilities::eModified_A, num_modes+1, PkeyH));


	graph3D->SetBasisKey(LibUtilities::eTetrahedron, BkeyT);
	graph3D->SetBasisKey(LibUtilities::ePrism, BkeyP);
	graph3D->SetBasisKey(LibUtilities::eHexahedron, BkeyH);

	MultiRegions::DisContField3DSharedPtr PostProc = 
		MemoryManager<MultiRegions::DisContField3D>::AllocateSharedPtr(vSession,graph3D,vSession->GetVariable(0));

	int ErrorCoordim = PostProc->GetCoordim(0);
	int ErrorNq      = PostProc->GetTotPoints();

	Array<OneD,NekDouble> ErrorXc0(ErrorNq,0.0);
	Array<OneD,NekDouble> ErrorXc1(ErrorNq,0.0);
	Array<OneD,NekDouble> ErrorXc2(ErrorNq,0.0);

	switch(ErrorCoordim)
	{
		case 1:
			PostProc->GetCoords(ErrorXc0);
			break;
		case 2:
			PostProc->GetCoords(ErrorXc0,ErrorXc1);
			break;
		case 3:
			PostProc->GetCoords(ErrorXc0,ErrorXc1,ErrorXc2);
			break;
	}
        
        
	// evaluate exact solution 
	Array<OneD,NekDouble> ppSol(ErrorNq);
	ex_sol->Evaluate(ErrorXc0,ErrorXc1,ErrorXc2,ppSol);

	// calcualte spectral/hp approximation on the quad points of this new
	// expansion basis
	std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef 
		= Exp->GetFieldDefinitions();
	std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
	std::string fieldstr = "u";

	for(i = 0; i < FieldDef.size(); ++i)
	{
		FieldDef[i]->m_fields.push_back(fieldstr);
		Exp->AppendFieldData(FieldDef[i], FieldData[i]);
		PostProc->ExtractDataToCoeffs(FieldDef[i],FieldData[i],fieldstr,PostProc->UpdateCoeffs());
	}

	// Interpolation of trace 
	std::vector<LibUtilities::FieldDefinitionsSharedPtr> TraceDef 
		= Exp->GetTrace()->GetFieldDefinitions();
	std::vector<std::vector<NekDouble> > TraceData(TraceDef.size());
	for(i = 0; i < TraceDef.size(); ++i)
	{
		TraceDef[i]->m_fields.push_back(fieldstr);
		Exp->GetTrace()->AppendFieldData(TraceDef[i], TraceData[i]);
		PostProc->GetTrace()->ExtractDataToCoeffs(TraceDef[i],TraceData[i],fieldstr,PostProc->GetTrace()->UpdateCoeffs());
	}
        
	PostProc->BwdTrans_IterPerExp(PostProc->GetCoeffs(),PostProc->UpdatePhys());

	PostProc->EvaluateHDGPostProcessing(PostProc->UpdateCoeffs());
	PostProc->BwdTrans_IterPerExp(PostProc->GetCoeffs(),PostProc->UpdatePhys());
	
	NekDouble vLinfError = Exp->Linf(Exp->GetPhys(), fce);
	NekDouble vL2Error   = Exp->L2  (Exp->GetPhys(), fce);
	NekDouble L2ErrorPostProc = PostProc->L2(PostProc->GetPhys(), ppSol);
	NekDouble LinfErrorPostProc = PostProc->Linf(PostProc->GetPhys(), ppSol); 

	if (vSession->GetComm()->GetRank() == 0)
	{
		cout << "L infinity error : " << vLinfError << endl;
		cout << "L 2 error        : " << vL2Error   << endl;
		cout << "Postprocessed L infinity error : " << LinfErrorPostProc << endl;
		cout << "Postprocessed L 2 error        : " << L2ErrorPostProc   << endl;
	}

	vSession->Finalise();
    
    return 0;
}

//.........这里部分代码省略.........
    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 ..");

#if 0
    for(i = 0; i < 100; ++i)
    {
        Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
    }
    
    Timing("100 Helmholtz Solves:... ");
#endif 

    //----------------------------------------------
    // Backward Transform Solution to get solved values at 
    Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys());
    //----------------------------------------------
    Timing("Backward Transform ..");
    
    //-----------------------------------------------
    // Write solution to file
    string out = vSession->GetSessionName() + ".fld";
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
                                                = Exp->GetFieldDefinitions();
    std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());

    for(i = 0; i < FieldDef.size(); ++i)
    {
        FieldDef[i]->m_fields.push_back("u");
        Exp->AppendFieldData(FieldDef[i], FieldData[i]);
    }
    fld->Write(out, FieldDef, FieldData);
    //-----------------------------------------------
    
    //----------------------------------------------
    // See if there is an exact solution, if so 
    // evaluate and plot errors
    LibUtilities::EquationSharedPtr ex_sol =
        vSession->GetFunction("ExactSolution", 0);
    
    if(ex_sol)
    {
        //----------------------------------------------
        // evaluate exact solution 
        ex_sol->Evaluate(xc0, xc1, xc2,  fce);

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

        //--------------------------------------------
        // Calculate L_inf error 
        Fce->SetPhys(fce);
        Fce->SetPhysState(true);

        NekDouble vLinfError = Exp->Linf(Exp->GetPhys(), Fce->GetPhys());
        NekDouble vL2Error   = Exp->L2  (Exp->GetPhys(), Fce->GetPhys());
        NekDouble vH1Error   = Exp->H1  (Exp->GetPhys(), Fce->GetPhys());

        if (vComm->GetRank() == 0)
        {
            cout << "L infinity error: " << vLinfError << endl;
            cout << "L 2 error       : " << vL2Error   << endl;
            cout << "H 1 error       : " << vH1Error   << endl;
        }
        //--------------------------------------------        
    }
    
    Timing("Output ..");

    //----------------------------------------------        
    
    vSession->Finalise();
    
    return 0;
}