C++ SolutionPtr类代码示例

// TODO: Merge with ConcolicRuntime::printSolution and put somewhere accessible by both.
void ConcolicTargetGenerator::printSolution(const SolutionPtr solution) const
{
    QStringList symbols = solution->symbols();

    foreach (QString var, symbols) {
        Symbolvalue value = solution->findSymbol(var);
        assert(value.found);

        // TODO: Only object type should be seen here?
        switch (value.kind) {
        case Symbolic::INT:
            Log::debug(QString("    %1 = %2").arg(var).arg(value.u.integer).toStdString());
            break;
        case Symbolic::BOOL:
            Log::debug(QString("    %1 = %2").arg(var).arg(value.u.boolean ? "true" : "false").toStdString());
            break;
        case Symbolic::STRING:
            if (value.string.empty()) {
                Log::debug(QString("    %1 = \"\"").arg(var).toStdString());
            } else {
                Log::debug(QString("    %1 = \"%2\"").arg(var).arg(value.string.c_str()).toStdString());
            }
            break;
        case Symbolic::OBJECT:
            Log::debug(QString("    %1 -> %2").arg(var).arg(value.string.c_str()).toStdString());
            break;
        default:
            Log::fatal(QString("Unimplemented value type encountered for variable %1 (%2)").arg(var).arg(value.kind).toStdString());
            exit(1);
        }
    }

/**
 * Implementation of the Clarke and Wright savings algorithm
 * This function modifies an existing solution by merging tours from the same depot if this results in a net
 * cost reduction. The tours are merged in order of net reduction from highest to lowest. This is primarily
 * used to initiate the solver with a sensible first guess.
 */
SolutionPtr cw_savings(SolutionPtr prevsol)
{
    set<TourPtr> newsol;

    for (auto depot : prevsol->get_problem()->get_depots())
    {
        vector<TourPtr> tours_vec = prevsol->tours_from_depot(depot);
        set<TourPtr> tours(tours_vec.begin(), tours_vec.end());

        while (true)
        {
            double reduction = -1;
            TourPtr left, right;

            for (auto t1 : tours)
            for (auto t2 : tours)
            {
                if (t1 == t2)
                    continue;

                double test_reduction = triangle(*t1->last(), *depot, *t2->first());

                // Only accept this merge if it doesn't break the daily cap condition
                if (   test_reduction > reduction
                    && t1->duration() + t2->duration() - test_reduction <= 
                               prevsol->get_problem()->get_daily_cap())
                {
                    reduction = test_reduction;
                    left = t1;
                    right = t2;
                }
            }

            if (reduction >= 0)
            {
                tours.erase(left);
                tours.erase(right);
                tours.insert(TourPtr(new Tour(*left, *right)));
            }
            else
                break;
        }

        newsol.insert(tours.begin(), tours.end());
    }

    return SolutionPtr(new Solution(prevsol->get_problem(), newsol));
}

void writePatchValues(double xMin, double xMax, double yMin, double yMax,
                      SolutionPtr solution, VarPtr u1, string filename,
                      int numPoints=100)
{
  FieldContainer<double> points = pointGrid(xMin,xMax,yMin,yMax,numPoints);
  FieldContainer<double> values(numPoints*numPoints);
  solution->solutionValues(values, u1->ID(), points);
  ofstream fout(filename.c_str());
  fout << setprecision(15);

  fout << "X = zeros(" << numPoints << ",1);\n";
  //    fout << "Y = zeros(numPoints);\n";
  fout << "U = zeros(" << numPoints << "," << numPoints << ");\n";
  for (int i=0; i<numPoints; i++)
  {
    fout << "X(" << i+1 << ")=" << points(i,0) << ";\n";
  }
  for (int i=0; i<numPoints; i++)
  {
    fout << "Y(" << i+1 << ")=" << points(i,1) << ";\n";
  }

  for (int i=0; i<numPoints; i++)
  {
    for (int j=0; j<numPoints; j++)
    {
      int pointIndex = i*numPoints + j;
      fout << "U("<<i+1<<","<<j+1<<")=" << values(pointIndex) << ";" << endl;
    }
  }
  fout.close();
}

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL); // initialize MPI
  
  Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options
  
  int numElements = 3;
  double xLeft = 0.0, xRight = 1.0;
  int polyOrder = 1, delta_k = 1;
  
  cmdp.setOption("numElements", &numElements );
  cmdp.setOption("polyOrder", &polyOrder );
  cmdp.setOption("delta_k", &delta_k );
  cmdp.setOption("xLeft", &xLeft );
  cmdp.setOption("xRight", &xRight );
  
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
  {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }
  int spaceDim = 1;
  bool conformingTraces = true; // conformingTraces argument has no effect in 1D
  PoissonFormulation poissonForm(spaceDim, conformingTraces);
  
  MeshPtr mesh = MeshFactory::intervalMesh(poissonForm.bf(), xLeft, xRight, numElements, polyOrder + 1, delta_k); // 1D equispaced
  
  RHSPtr rhs = RHS::rhs(); // zero RHS
  IPPtr ip = poissonForm.bf()->graphNorm();
  BCPtr bc = BC::bc();
  bc->addDirichlet(poissonForm.phi_hat(), SpatialFilter::allSpace(), Function::zero());
  
  SolutionPtr solution = Solution::solution(poissonForm.bf(), mesh, bc, rhs, ip);
  solution->solve();
  
  GDAMinimumRule* minRule = dynamic_cast<GDAMinimumRule*>(mesh->globalDofAssignment().get());
//  minRule->printGlobalDofInfo();
  
  Teuchos::RCP<Epetra_CrsMatrix> A = solution->getStiffnessMatrix();
  EpetraExt::RowMatrixToMatrixMarketFile("A.dat",*A, NULL, NULL, false);
  
  HDF5Exporter exporter(mesh);
  
  return 0;
}

FieldContainer<double> solutionData(FieldContainer<double> &points, SolutionPtr solution, VarPtr u1) {
  int numPoints = points.dimension(0);
  FieldContainer<double> values(numPoints);
  solution->solutionValues(values, u1->ID(), points);
  
  FieldContainer<double> xyzData(numPoints, 3);
  for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
    xyzData(ptIndex,0) = points(ptIndex,0);
    xyzData(ptIndex,1) = points(ptIndex,1);
    xyzData(ptIndex,2) = values(ptIndex);
  }
  return xyzData;
}

void writeSol(EnvPtr env, VarVector *orig_v,
              PresolverPtr pres, SolutionPtr sol, SolveStatus status,
              MINOTAUR_AMPL::AMPLInterface* iface)
{
  if (sol) {
    sol = pres->getPostSol(sol);
  }

  if (env->getOptions()->findFlag("AMPL")->getValue() ||
      true == env->getOptions()->findBool("write_sol_file")->getValue()) {
    iface->writeSolution(sol, status);
  } else if (sol && env->getLogger()->getMaxLevel()>=LogExtraInfo &&
             env->getOptions()->findBool("display_solution")->getValue()) {
    sol->writePrimal(env->getLogger()->msgStream(LogExtraInfo), orig_v);
  }
}

//.........这里部分代码省略.........
  VarPtr tau2 = vf->testVar("tau2", HDIV);
  VarPtr q = vf->testVar("q", HGRAD);

  BFPtr bf = Teuchos::rcp( new BF(vf) );

  bf->addTerm(1./mu*sigma1, tau1);
  bf->addTerm(1./mu*sigma2, tau2);
  bf->addTerm(u1, tau1->div());
  bf->addTerm(u2, tau2->div());
  bf->addTerm(-u1hat, tau1->dot_normal());
  bf->addTerm(-u2hat, tau2->dot_normal());

  bf->addTerm(sigma1, v1->grad());
  bf->addTerm(sigma2, v2->grad());
  bf->addTerm(-p, v1->dx());
  bf->addTerm(-p, v2->dy());
  bf->addTerm(t1c, v1);
  bf->addTerm(t2c, v2);
  bf->addTerm(mu*permInv*u1, v1);
  bf->addTerm(mu*permInv*u2, v2);

  bf->addTerm(-u1, q->dx());
  bf->addTerm(-u2, q->dy());
  bf->addTerm(u1hat, q->times_normal_x());
  bf->addTerm(u2hat, q->times_normal_y());

  RHSPtr rhs = RHS::rhs();

  BCPtr bc = BC::bc();

  SpatialFilterPtr y_equals_one = SpatialFilter::matchingY(1.0);
  SpatialFilterPtr y_equals_zero = SpatialFilter::matchingY(0);
  SpatialFilterPtr x_equals_one = SpatialFilter::matchingX(1.0);
  SpatialFilterPtr x_equals_zero = SpatialFilter::matchingX(0.0);
  bc->addDirichlet(u1hat, y_equals_zero, u1_exact);
  bc->addDirichlet(u2hat, y_equals_zero, u2_exact);
  bc->addDirichlet(u1hat, x_equals_zero, u1_exact);
  bc->addDirichlet(u2hat, x_equals_zero, u2_exact);
  bc->addDirichlet(u1hat, y_equals_one, u1_exact);
  bc->addDirichlet(u2hat, y_equals_one, u2_exact);
  bc->addDirichlet(u1hat, x_equals_one, u1_exact);
  bc->addDirichlet(u2hat, x_equals_one, u2_exact);
  bc->addZeroMeanConstraint(p);

  MeshPtr mesh = MeshFactory::quadMesh(bf, k+1, delta_k, 1, 1, 4, 4);

  map<string, IPPtr> brinkmanIPs;
  brinkmanIPs["Graph"] = bf->graphNorm();

  brinkmanIPs["Decoupled"] = Teuchos::rcp(new IP);
  brinkmanIPs["Decoupled"]->addTerm(tau1);
  brinkmanIPs["Decoupled"]->addTerm(tau2);
  brinkmanIPs["Decoupled"]->addTerm(tau1->div());
  brinkmanIPs["Decoupled"]->addTerm(tau2->div());
  brinkmanIPs["Decoupled"]->addTerm(permInv*v1);
  brinkmanIPs["Decoupled"]->addTerm(permInv*v2);
  brinkmanIPs["Decoupled"]->addTerm(v1->grad());
  brinkmanIPs["Decoupled"]->addTerm(v2->grad());
  brinkmanIPs["Decoupled"]->addTerm(q);
  brinkmanIPs["Decoupled"]->addTerm(q->grad());

  // brinkmanIPs["CoupledRobust"] = Teuchos::rcp(new IP);
  // brinkmanIPs["CoupledRobust"]->addTerm(tau->div()-beta*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(Function<double>::min(one/Function<double>::h(),Function<double>::constant(1./sqrt(epsilon)))*tau);
  // brinkmanIPs["CoupledRobust"]->addTerm(sqrt(epsilon)*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(beta*v->grad());
  // brinkmanIPs["CoupledRobust"]->addTerm(Function<double>::min(sqrt(epsilon)*one/Function<double>::h(),one)*v);

  IPPtr ip = brinkmanIPs[norm];

  SolutionPtr soln = TSolution<double>::solution(mesh, bc, rhs, ip);

  double threshold = 0.20;
  RefinementStrategy refStrategy(soln, threshold);

  ostringstream refName;
  refName << "brinkman";
  HDF5Exporter exporter(mesh,refName.str());

  for (int refIndex=0; refIndex <= numRefs; refIndex++)
  {
    soln->solve(false);

    double energyError = soln->energyErrorTotal();
    if (commRank == 0)
    {
      // if (refIndex > 0)
      // refStrategy.printRefinementStatistics(refIndex-1);
      cout << "Refinement:\t " << refIndex << " \tElements:\t " << mesh->numActiveElements()
           << " \tDOFs:\t " << mesh->numGlobalDofs() << " \tEnergy Error:\t " << energyError << endl;
    }

    exporter.exportSolution(soln, refIndex);

    if (refIndex != numRefs)
      refStrategy.refine();
  }

  return 0;
}

//.........这里部分代码省略.........
  {
    u_1 = varFactory.fieldVar("u_1");
    u_2 = varFactory.fieldVar("u_2");
    u_3 = varFactory.fieldVar("rho");
    u_4 = varFactory.fieldVar("T");
  }

  // test fxns
  VarPtr tau1 = varFactory.testVar("\\tau_1",HDIV);
  VarPtr tau2 = varFactory.testVar("\\tau_2",HDIV);
  VarPtr tau3 = varFactory.testVar("\\tau_3",HDIV);
  VarPtr v1 = varFactory.testVar("v_1",HGRAD);
  VarPtr v2 = varFactory.testVar("v_2",HGRAD);
  VarPtr v3 = varFactory.testVar("v_3",HGRAD);
  VarPtr v4 = varFactory.testVar("v_4",HGRAD);

  ////////////////////////////////////////////////////////////////////
  // CREATE BILINEAR FORM PTR AND MESH
  ////////////////////////////////////////////////////////////////////

  BFPtr bf = Teuchos::rcp( new BF(varFactory) ); // initialize bilinear form

  // create a pointer to a new mesh:
  //  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(domainPoints, horizontalCells, verticalCells, bf, H1Order, H1Order+pToAdd, useTriangles);
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildRampMesh(rampHeight,bf, H1Order, H1Order+pToAdd);
  //  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));
  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("REFTREE")));
  MeshInfo meshInfo(mesh); // gets info like cell measure, etc

  // for writing ref history to file
  Teuchos::RCP< RefinementHistory > refHistory = Teuchos::rcp( new RefinementHistory );
  mesh->registerObserver(refHistory);


  ////////////////////////////////////////////////////////////////////
  // CREATE SOLUTION OBJECT
  ////////////////////////////////////////////////////////////////////

  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  SolutionPtr prevTimeFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  Teuchos::RCP<Solution> solution = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP));
  int enrichDegree = 2; // just for kicks.
  if (rank==0)
    cout << "enriching cubature by " << enrichDegree << endl;
  solution->setCubatureEnrichmentDegree(enrichDegree); // double cubature enrichment
  if (reportTimingResults)
  {
    solution->setReportTimingResults(true);
  }

  mesh->registerSolution(solution);
  mesh->registerSolution(backgroundFlow); // u_t(i)
  mesh->registerSolution(prevTimeFlow); // u_t(i-1)

  // for loading refinement history
  if (replayFile.length() > 0)
  {
    RefinementHistory refHistory;
    replayFile = dir + replayFile;
    refHistory.loadFromFile(replayFile);
    refHistory.playback(mesh);
    int numElems = mesh->numActiveElements();
    if (rank==0)
    {
      double minSideLength = meshInfo.getMinCellSideLength() ;
      cout << "after replay, num elems = " << numElems << " and min side length = " << minSideLength << endl;
    }
  }
  if (solnLoadFile.length() > 0)
  {
    std::ostringstream ss;
    //    ss << dir <<  "solution_" << solnLoadFile;
    //    solution->readFromFile(ss.str());
    ss.str("");
    ss << dir << "backgroundFlow_" << solnLoadFile;
    backgroundFlow->readFromFile(ss.str());
    ss.str("");
    ss << dir << "prevTimeFlow_" << solnLoadFile;
    prevTimeFlow->readFromFile(ss.str());
  }

  ////////////////////////////////////////////////////////////////////
  // PSEUDO-TIME SOLVE STRATEGY
  ////////////////////////////////////////////////////////////////////

  VTKExporter exporter(solution, mesh, varFactory);
  VTKExporter backgroundFlowExporter(backgroundFlow, mesh, varFactory);


  if (rank==0)
  {
    exporter.exportSolution("dU");
    backgroundFlowExporter.exportSolution("U");
  }

  return 0;
}

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

  double minL2Increment = 1e-8;

  // get variable definitions:
  VarFactory varFactory = VGPStokesFormulation::vgpVarFactory();
  u1 = varFactory.fieldVar(VGP_U1_S);
  u2 = varFactory.fieldVar(VGP_U2_S);
  sigma11 = varFactory.fieldVar(VGP_SIGMA11_S);
  sigma12 = varFactory.fieldVar(VGP_SIGMA12_S);
  sigma21 = varFactory.fieldVar(VGP_SIGMA21_S);
  sigma22 = varFactory.fieldVar(VGP_SIGMA22_S);
  p = varFactory.fieldVar(VGP_P_S);

  u1hat = varFactory.traceVar(VGP_U1HAT_S);
  u2hat = varFactory.traceVar(VGP_U2HAT_S);
  t1n = varFactory.fluxVar(VGP_T1HAT_S);
  t2n = varFactory.fluxVar(VGP_T2HAT_S);

  v1 = varFactory.testVar(VGP_V1_S, HGRAD);
  v2 = varFactory.testVar(VGP_V2_S, HGRAD);
  tau1 = varFactory.testVar(VGP_TAU1_S, HDIV);
  tau2 = varFactory.testVar(VGP_TAU2_S, HDIV);
  q = varFactory.testVar(VGP_Q_S, HGRAD);

  FunctionPtr u1_0 = Teuchos::rcp( new U1_0(eps) );
  FunctionPtr u2_0 = Teuchos::rcp( new U2_0 );
  FunctionPtr zero = Function::zero();
  ParameterFunctionPtr Re_param = ParameterFunction::parameterFunction(1);
  VGPNavierStokesProblem problem = VGPNavierStokesProblem(Re_param,quadPoints,
                                   horizontalCells,verticalCells,
                                   H1Order, pToAdd,
                                   u1_0, u2_0,  // BC for u
                                   zero, zero); // zero forcing function
  SolutionPtr solution = problem.backgroundFlow();
  SolutionPtr solnIncrement = problem.solutionIncrement();

  Teuchos::RCP<Mesh> mesh = problem.mesh();
  mesh->registerSolution(solution);
  mesh->registerSolution(solnIncrement);

  ///////////////////////////////////////////////////////////////////////////

  // define bilinear form for stream function:
  VarFactory streamVarFactory;
  VarPtr phi_hat = streamVarFactory.traceVar("\\widehat{\\phi}");
  VarPtr psin_hat = streamVarFactory.fluxVar("\\widehat{\\psi}_n");
  VarPtr psi_1 = streamVarFactory.fieldVar("\\psi_1");
  VarPtr psi_2 = streamVarFactory.fieldVar("\\psi_2");
  VarPtr phi = streamVarFactory.fieldVar("\\phi");
  VarPtr q_s = streamVarFactory.testVar("q_s", HGRAD);
  VarPtr v_s = streamVarFactory.testVar("v_s", HDIV);
  BFPtr streamBF = Teuchos::rcp( new BF(streamVarFactory) );
  streamBF->addTerm(psi_1, q_s->dx());
  streamBF->addTerm(psi_2, q_s->dy());
  streamBF->addTerm(-psin_hat, q_s);

  streamBF->addTerm(psi_1, v_s->x());
  streamBF->addTerm(psi_2, v_s->y());
  streamBF->addTerm(phi, v_s->div());
  streamBF->addTerm(-phi_hat, v_s->dot_normal());

  Teuchos::RCP<Mesh> streamMesh, overkillMesh;

  streamMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
                                          streamBF, H1Order+pToAddForStreamFunction,
                                          H1Order+pToAdd+pToAddForStreamFunction, useTriangles);

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL); // initialize MPI
  
  Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options
  
  int numElements = 3;
  vector<vector<double>> domainDim(3,vector<double>{0.0,1.0}); // first index: spaceDim; second: 0/1 for x0, x1, etc.
  int polyOrder = 2, delta_k = 1;
  int spaceDim = 2;
  
  cmdp.setOption("numElements", &numElements );
  cmdp.setOption("polyOrder", &polyOrder );
  cmdp.setOption("delta_k", &delta_k );
  cmdp.setOption("x0", &domainDim[0][0] );
  cmdp.setOption("x1", &domainDim[0][1] );
  cmdp.setOption("y0", &domainDim[1][0] );
  cmdp.setOption("y1", &domainDim[1][1] );
  cmdp.setOption("z0", &domainDim[2][0] );
  cmdp.setOption("z1", &domainDim[2][1] );
  cmdp.setOption("spaceDim", &spaceDim);
  
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
  {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }
  
  vector<double> x0(spaceDim);
  vector<double> domainSize(spaceDim);
  vector<int> elementCounts(spaceDim);
  for (int d=0; d<spaceDim; d++)
  {
    x0[d] = domainDim[d][0];
    domainSize[d] = domainDim[d][1] - x0[d];
    elementCounts[d] = numElements;
  }
  
  bool conformingTraces = true; // no difference for primal/continuous formulations
  PoissonFormulation formCG(spaceDim, conformingTraces, PoissonFormulation::CONTINUOUS_GALERKIN);
  VarPtr q = formCG.q();
  VarPtr phi = formCG.phi();
  BFPtr bf = formCG.bf();
  
  MeshPtr bubnovMesh = MeshFactory::rectilinearMesh(bf, domainSize, elementCounts, polyOrder, 0, x0);

  // Right now, hanging nodes don't work with continuous field variables
  // there is a GDAMinimumRule test demonstrating the failure, SolvePoisson2DContinuousGalerkinHangingNode.
  // make a mesh with hanging nodes (when spaceDim > 1)
//  {
//    set<GlobalIndexType> cellsToRefine = {0};
//    bubnovMesh->hRefine(cellsToRefine);
//  }

  RHSPtr rhs = RHS::rhs();
  rhs->addTerm(1.0 * q); // unit forcing
  
  IPPtr ip = Teuchos::null; // will give Bubnov-Galerkin
  BCPtr bc = BC::bc();
  bc->addDirichlet(phi, SpatialFilter::allSpace(), Function::zero());
  
  SolutionPtr solution = Solution::solution(bf, bubnovMesh, bc, rhs, ip);
  solution->solve();
  
  HDF5Exporter exporter(bubnovMesh, "PoissonContinuousGalerkin");
  exporter.exportSolution(solution);

  /**** Sort-of-primal experiment ****/
  // an experiment: try doing "primal" DPG with IBP to the boundary
//  ip = IP::ip();
//  ip->addTerm(q->grad());
//  ip->addTerm(q);
//
//  solution = Solution::solution(bf, bubnovMesh, bc, rhs, ip);
//  solution->solve();
//  
//  HDF5Exporter primalNoFluxExporter(bubnovMesh, "PoissonPrimalNoFlux");
//  primalNoFluxExporter.exportSolution(solution);
  
  //*** Primal Formulation ***//
  PoissonFormulation form(spaceDim, conformingTraces, PoissonFormulation::PRIMAL);
  q = form.q();
  phi = form.phi();
  bf = form.bf();
  
  bc = BC::bc();
  bc->addDirichlet(phi, SpatialFilter::allSpace(), Function::zero());
  
  rhs = RHS::rhs();
  rhs->addTerm(1.0 * q); // unit forcing
  
  MeshPtr primalMesh = MeshFactory::rectilinearMesh(bf, domainSize, elementCounts, polyOrder, delta_k, x0);
  
  ip = IP::ip();
  ip->addTerm(q->grad());
  ip->addTerm(q);

  // Right now, hanging nodes don't work with continuous field variables
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif
  int polyOrder = 3;
  int pToAdd = 2; // for tests

  // define our manufactured solution or problem bilinear form:
  bool useTriangles = false;

  FieldContainer<double> meshPoints(4,2);

  meshPoints(0,0) = 0.0; // x1
  meshPoints(0,1) = 0.0; // y1
  meshPoints(1,0) = 1.0;
  meshPoints(1,1) = 0.0;
  meshPoints(2,0) = 1.0;
  meshPoints(2,1) = 1.0;
  meshPoints(3,0) = 0.0;
  meshPoints(3,1) = 1.0;

  int H1Order = polyOrder + 1;
  int horizontalCells = 4, verticalCells = 4;

  double energyThreshold = 0.2; // for mesh refinements
  double nonlinearStepSize = 0.5;
  double nonlinearRelativeEnergyTolerance = 1e-8; // used to determine convergence of the nonlinear solution

  ////////////////////////////////////////////////////////////////////
  // DEFINE VARIABLES
  ////////////////////////////////////////////////////////////////////

  // new-style bilinear form definition
  VarFactory varFactory;
  VarPtr fhat = varFactory.fluxVar("\\widehat{f}");
  VarPtr u = varFactory.fieldVar("u");

  VarPtr v = varFactory.testVar("v",HGRAD);
  BFPtr bf = Teuchos::rcp( new BF(varFactory) ); // initialize bilinear form

  ////////////////////////////////////////////////////////////////////
  // CREATE MESH
  ////////////////////////////////////////////////////////////////////

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshPoints, horizontalCells,
                            verticalCells, bf, H1Order,
                            H1Order+pToAdd, useTriangles);
  mesh->setPartitionPolicy(Teuchos::rcp(new ZoltanMeshPartitionPolicy("HSFC")));

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////
  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC,
                               nullRHS, nullIP) );

  vector<double> e1(2); // (1,0)
  e1[0] = 1;
  vector<double> e2(2); // (0,1)
  e2[1] = 1;

  FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
  FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );

  ////////////////////////////////////////////////////////////////////
  // DEFINE BILINEAR FORM
  ////////////////////////////////////////////////////////////////////

  // v:
  // (sigma, grad v)_K - (sigma_hat_n, v)_dK - (u, beta dot grad v) + (u_hat * n dot beta, v)_dK
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( fhat, v);

  // ==================== SET INITIAL GUESS ==========================
  mesh->registerSolution(backgroundFlow);
  FunctionPtr zero = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  FunctionPtr u0 = Teuchos::rcp( new U0 );

  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()] = u0;

  backgroundFlow->projectOntoMesh(functionMap);
  // ==================== END SET INITIAL GUESS ==========================

  ////////////////////////////////////////////////////////////////////
  // DEFINE INNER PRODUCT
  ////////////////////////////////////////////////////////////////////
  IPPtr ip = Teuchos::rcp( new IP );
  ip->addTerm( v );
  ip->addTerm( beta * v->grad() );

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

int main(int argc, char *argv[])
{
  // Process command line arguments
#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  int rank=mpiSession.getRank();
  int numProcs=mpiSession.getNProc();
#else
  int rank = 0;
  int numProcs = 1;
#endif

  ////////////////////   DECLARE VARIABLES   ///////////////////////
  // define test variables
  VarFactory varFactory;
  VarPtr v = varFactory.testVar("v", HGRAD);

  // define trial variables
  VarPtr beta_n_u_hat = varFactory.fluxVar("\\widehat{\\beta \\cdot n }");
  VarPtr u = varFactory.fieldVar("u");

  FunctionPtr beta = Teuchos::rcp(new Beta());

  ////////////////////   BUILD MESH   ///////////////////////
  BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
  // define nodes for mesh
  FieldContainer<double> meshBoundary(4,2);

  meshBoundary(0,0) = -1.0; // x1
  meshBoundary(0,1) = -1.0; // y1
  meshBoundary(1,0) =  1.0;
  meshBoundary(1,1) = -1.0;
  meshBoundary(2,0) =  1.0;
  meshBoundary(2,1) =  1.0;
  meshBoundary(3,0) = -1.0;
  meshBoundary(3,1) =  1.0;

  int horizontalCells = 32, verticalCells = 32;

  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = Mesh::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
                            confusionBF, H1Order, H1Order+pToAdd);

  ////////////////////////////////////////////////////////////////////
  // INITIALIZE FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////

  BCPtr nullBC = Teuchos::rcp((BC*)NULL);
  RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL);
  IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr prevTimeFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  SolutionPtr flowResidual = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );

  FunctionPtr u_prev_time = Teuchos::rcp( new PreviousSolutionFunction(prevTimeFlow, u) );

  ////////////////////   DEFINE BILINEAR FORM   ///////////////////////
  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr invDt = Teuchos::rcp(new ScalarParamFunction(1.0/dt));

  // v terms:
  confusionBF->addTerm( beta * u, - v->grad() );
  confusionBF->addTerm( beta_n_u_hat, v);

  confusionBF->addTerm( u, invDt*v );
  rhs->addTerm( u_prev_time * invDt * v );

  ////////////////////   SPECIFY RHS   ///////////////////////
  FunctionPtr f = Teuchos::rcp( new ConstantScalarFunction(0.0) );
  rhs->addTerm( f * v ); // obviously, with f = 0 adding this term is not necessary!

  ////////////////////   DEFINE INNER PRODUCT(S)   ///////////////////////
  // robust test norm
  IPPtr ip = confusionBF->graphNorm();
  // IPPtr ip = Teuchos::rcp(new IP);
  // ip->addTerm(v);
  // ip->addTerm(invDt*v - beta*v->grad());

  ////////////////////   CREATE BCs   ///////////////////////
  Teuchos::RCP<BCEasy> bc = Teuchos::rcp( new BCEasy );
  SpatialFilterPtr inflowBoundary = Teuchos::rcp( new InflowSquareBoundary(beta) );
  FunctionPtr u0 = Teuchos::rcp( new ConstantScalarFunction(0) );
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );

  bc->addDirichlet(beta_n_u_hat, inflowBoundary, beta*n*u0);

  Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );

  // ==================== Register Solutions ==========================
  mesh->registerSolution(solution);
  mesh->registerSolution(prevTimeFlow);
  mesh->registerSolution(flowResidual);

  // ==================== SET INITIAL GUESS ==========================
  FunctionPtr u_init = Teuchos::rcp(new InitialCondition());
  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()]      = u_init;

  prevTimeFlow->projectOntoMesh(functionMap);

  ////////////////////   SOLVE & REFINE   ///////////////////////
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
  Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
  choice::MpiArgs args( argc, argv );
#else
  choice::Args args( argc, argv );
#endif
  int rank = Teuchos::GlobalMPISession::getRank();
  int numProcs = Teuchos::GlobalMPISession::getNProc();

  int nCells = args.Input<int>("--nCells", "num cells",2);
  int numSteps = args.Input<int>("--numSteps", "num NR steps",20);

  int polyOrder = 0;
  
  // define our manufactured solution or problem bilinear form:
  bool useTriangles = false;
  
  int pToAdd = 1;

  args.Process();

  int H1Order = polyOrder + 1;
  
  ////////////////////////////////////////////////////////////////////
  // DEFINE VARIABLES 
  ////////////////////////////////////////////////////////////////////
  
  // new-style bilinear form definition
  VarFactory varFactory;
  VarPtr fn = varFactory.fluxVar("\\widehat{\\beta_n_u}");
  VarPtr u = varFactory.fieldVar("u");
  
  VarPtr v = varFactory.testVar("v",HGRAD);
  BFPtr bf = Teuchos::rcp( new BF(varFactory) ); // initialize bilinear form
  
  ////////////////////////////////////////////////////////////////////
  // CREATE MESH 
  ////////////////////////////////////////////////////////////////////
  
  // create a pointer to a new mesh:
  Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells , bf, H1Order, H1Order+pToAdd);
  
  ////////////////////////////////////////////////////////////////////
  // INITIALIZE BACKGROUND FLOW FUNCTIONS
  ////////////////////////////////////////////////////////////////////
  BCPtr nullBC = Teuchos::rcp((BC*)NULL); RHSPtr nullRHS = Teuchos::rcp((RHS*)NULL); IPPtr nullIP = Teuchos::rcp((IP*)NULL);
  SolutionPtr backgroundFlow = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  SolutionPtr solnPerturbation = Teuchos::rcp(new Solution(mesh, nullBC, nullRHS, nullIP) );
  
  vector<double> e1(2),e2(2);
  e1[0] = 1; e2[1] = 1;
  
  FunctionPtr u_prev = Teuchos::rcp( new PreviousSolutionFunction(backgroundFlow, u) );
  FunctionPtr beta = e1 * u_prev + Teuchos::rcp( new ConstantVectorFunction( e2 ) );
  
  ////////////////////////////////////////////////////////////////////
  // DEFINE BILINEAR FORM
  ////////////////////////////////////////////////////////////////////
  
  // v:
  bf->addTerm( -u, beta * v->grad());
  bf->addTerm( fn, v);

  ////////////////////////////////////////////////////////////////////
  // DEFINE RHS
  ////////////////////////////////////////////////////////////////////

  Teuchos::RCP<RHSEasy> rhs = Teuchos::rcp( new RHSEasy );
  FunctionPtr u_prev_squared_div2 = 0.5 * u_prev * u_prev;  
  rhs->addTerm((e1 * u_prev_squared_div2 + e2 * u_prev) * v->grad());

  // ==================== SET INITIAL GUESS ==========================

  mesh->registerSolution(backgroundFlow);
  FunctionPtr zero = Function::constant(0.0);
  FunctionPtr u0 = Teuchos::rcp( new U0 );
  FunctionPtr n = Teuchos::rcp( new UnitNormalFunction );
  //  FunctionPtr parity = Teuchos::rcp(new SideParityFunction);

  FunctionPtr u0_squared_div_2 = 0.5 * u0 * u0;

  map<int, Teuchos::RCP<Function> > functionMap;
  functionMap[u->ID()] = u0;
  //  functionMap[fn->ID()] = -(e1 * u0_squared_div_2 + e2 * u0) * n * parity;
  backgroundFlow->projectOntoMesh(functionMap);

  // ==================== END SET INITIAL GUESS ==========================

  ////////////////////////////////////////////////////////////////////
  // DEFINE INNER PRODUCT
  ////////////////////////////////////////////////////////////////////

  IPPtr ip = Teuchos::rcp( new IP );
  ip->addTerm( v );
  ip->addTerm(v->grad());
  //  ip->addTerm( beta * v->grad() ); // omitting term to make IP non-dependent on u

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

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

    /**************************************************************/
    // Given problem statistics .
    // Number of variables.
    UInt numvars = minlp->getNumVars();
    // number of constraints.
    UInt numcons = minlp->getNumCons();
    // linear constraints.
    UInt numlin = minlp->getNumLinCons();
    /*************************************************************/

    // set option for engine to resolve to solve NLP repeatedly.
    // Probbaly does nothing.
    e.setOptionsForRepeatedSolve();

    // load problem.
    e.load(minlp);
      
    // Solve problem.
    timer->start();
    status = e.solve();

    /********************************************************************/
    // Solution time of relaxation.
    Double timeinit = timer->query();
    timer->stop();
    // Solution objective value
    Double initobj = e.getSolutionValue();
    /********************************************************************/

    std::cout << "Relaxation objective value = " << initobj << std::endl; 
  
    // Get solution from engine.
    ConstSolutionPtr sol = e.getSolution();
  
    // Construct relaxation.
    RelaxationPtr rel = (RelaxationPtr) new Relaxation(minlp);
    
    // Time for cut generation.
    timer->start();
    // Generate kanpsack cover cuts.
    CoverCutGeneratorPtr knapgen = 
      (CoverCutGeneratorPtr) new CoverCutGenerator(rel, sol, env);

    /*******************************************************************/
    Double timecut = timer->query();
    timer->stop();
    /*******************************************************************/


    // Get statistics of cut generator.
    ConstCovCutGenStatsPtr knapstats = knapgen->getStats();

    /*******************************************************************/
    // Knapsack cut generator statistics.
    // knapsack constraints.
    UInt numknap = (knapgen->getKnapsackList())->getNumKnaps();
    // knapsacks that has cover sets.
    UInt numknapcov = knapgen->getNumCons();
    // knapsack subproblems solved, i.e number of lifting subproblems solved.
    UInt knaps = knapstats->knaps;
    // cover cuts including duplicates.
    UInt totalcuts = knapstats->totalcuts;
    // cuts without duplicates.
    UInt numknapcuts = knapstats->cuts;
    // violated cuts.

//.........这里部分代码省略.........
    if (loadSolution)
    {
      loadFilePrefix = rootDir + "/" + loadDir.str() + "/" + saveDir.str();
      if (commRank == 0) cout << "Loading previous solution " << loadFilePrefix << endl;
    }
    // ostringstream saveDir;
    // saveDir << problemName.str() << "_ref" << loadRef;
    string saveFilePrefix = rootDir + "/" + saveDir.str() + "/" + problemName.str();
    if (saveSolution && commRank == 0) cout << "Saving to " << saveFilePrefix << endl;

    Teuchos::ParameterList parameters;
    parameters.set("spaceDim", spaceDim);
    parameters.set("steady", steady);
    parameters.set("mu", 1./Re);
    parameters.set("useConformingTraces", useConformingTraces);
    parameters.set("fieldPolyOrder", p);
    parameters.set("delta_p", delta_p);
    parameters.set("numTElems", numTElems);
    parameters.set("norm", norm);
    parameters.set("savedSolutionAndMeshPrefix", loadFilePrefix);
    SpaceTimeIncompressibleFormulationPtr form = Teuchos::rcp(new SpaceTimeIncompressibleFormulation(problem, parameters));

    MeshPtr mesh = form->solutionUpdate()->mesh();
    vector<MeshPtr> meshesCoarseToFine;
    MeshPtr k0Mesh = Teuchos::rcp( new Mesh (mesh->getTopology()->deepCopy(), form->bf(), 1, delta_p) );
    meshesCoarseToFine.push_back(k0Mesh);
    meshesCoarseToFine.push_back(mesh);
    // mesh->registerObserver(k0Mesh);

    // Set up boundary conditions
    problem->setBCs(form);

    // Set up solution
    SolutionPtr solutionUpdate = form->solutionUpdate();
    SolutionPtr solutionBackground = form->solutionBackground();
    // dynamic_cast<AnalyticalIncompressibleProblem*>(problem.get())->projectExactSolution(solutionBackground);

    RefinementStrategyPtr refStrategy = form->getRefinementStrategy();
    Teuchos::RCP<HDF5Exporter> exporter;
    if (exportSolution)
      exporter = Teuchos::rcp(new HDF5Exporter(mesh,exportName, rootDir));

    Teuchos::RCP<Time> solverTime = Teuchos::TimeMonitor::getNewCounter("Solve Time");
    map<string, SolverPtr> solvers;
    solvers["KLU"] = Solver::getSolver(Solver::KLU, true);
#if defined(HAVE_AMESOS_SUPERLUDIST) || defined(HAVE_AMESOS2_SUPERLUDIST)
    solvers["SuperLUDist"] = Solver::getSolver(Solver::SuperLUDist, true);
#endif
#ifdef HAVE_AMESOS_MUMPS
    solvers["MUMPS"] = Solver::getSolver(Solver::MUMPS, true);
#endif
    bool useStaticCondensation = false;

    GMGOperator::MultigridStrategy multigridStrategy;
    if (multigridStrategyString == "Two-level")
    {
      multigridStrategy = GMGOperator::TWO_LEVEL;
    }
    else if (multigridStrategyString == "W-cycle")
    {
      multigridStrategy = GMGOperator::W_CYCLE;
    }
    else if (multigridStrategyString == "V-cycle")
    {
      multigridStrategy = GMGOperator::V_CYCLE;
    }