C++ MeshSource::getMesh方法代码示例

Mesh MeshBuilder::createMesh(const ParameterList& params)
{
  TEST_FOR_EXCEPTION(!params.isParameter("type"), RuntimeError,
                     "field name 'type' expected but not found in MeshBuilder "
                     "input parameter list: " << params);

  std::string type = params.get<string>("type");

  MeshSource mesher;

  if (type=="Rectangle")
    {
      mesher = new PartitionedRectangleMesher(params);
    }
  else if (type=="Line")
    {
      mesher = new PartitionedLineMesher(params);
    }
  else if (type=="Exodus")
    {
      mesher = new ExodusNetCDFMeshReader(params);
    }
  else if (type=="Triangle")
    {
      mesher = new TriangleMeshReader(params);
    }

  TEST_FOR_EXCEPTION(mesher.ptr().get()==0, RuntimeError,
                     "null mesh source in MeshBuilder::createMesh()");

  return mesher.getMesh();
}

int main(int argc, char** argv)
{
  
  try
		{
      GlobalMPISession session(&argc, &argv);

      TimeMonitor t(totalTimer());

      MeshType meshType = new BasicSimplicialMeshType();

      MeshSource mesher = new PartitionedRectangleMesher(0.0, 1.0, 32, 1,
                                                         0.0, 1.0, 32, 1,
                                                         meshType);

      Mesh mesh2D = mesher.getMesh();

      MeshTransformation extruder = new ExtrusionMeshTransformation(0.0, 1.0, 32, meshType);

      Mesh mesh3D = extruder.apply(mesh2D);

      FieldWriter w3 = new VTKWriter("test3d");

      w3.addMesh(mesh3D);

      w3.write();

      std::cout << "num elements = " << mesh3D.numCells(3) << std::endl;
      std::cout << "num nodes = " << mesh3D.numCells(0) << std::endl;

      TimeMonitor::summarize();
    }
	catch(std::exception& e)
		{
      std::cerr << e.what() << std::endl;
		}
}

int main(int argc, char** argv)
{
  try
		{
      Sundance::init(&argc, &argv);
      int np = MPIComm::world().getNProc();

      /* We will do our linear algebra using Epetra */
      VectorType<double> vecType = new EpetraVectorType();

      /* Create a mesh. It will be of type BasisSimplicialMesh, and will
       * be built using a PartitionedLineMesher. */
      MeshType meshType = new BasicSimplicialMeshType();
      MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, 1*np, meshType);
      Mesh mesh = mesher.getMesh();

      /* Create a cell filter that will identify the maximal cells
       * in the interior of the domain */
      CellFilter interior = new MaximalCellFilter();

      /* Create the Spectral Basis */
      int ndim = 1;
      int order = 2;
      SpectralBasis sbasis = new HermiteSpectralBasis(ndim, order); 
      
      /* Create unknown and test functions, discretized using first-order
       * Lagrange interpolants */
      Expr u = new UnknownFunction(new Lagrange(1), sbasis, "u");
      Expr v = new TestFunction(new Lagrange(1), sbasis, "v");

      /* Create the stochastic input function. */
      Expr a0 = new Sundance::Parameter(1.0);
      Expr a1 = new Sundance::Parameter(0.1);
      Expr a2 = new Sundance::Parameter(0.01);
      Expr alpha = new SpectralExpr(sbasis, tuple(a0, a1, a2));

      /* Create a discrete space, and discretize the function 1.0 on it */
      cout << "forming discrete space" << std::endl;
      DiscreteSpace discSpace(mesh, new Lagrange(1), sbasis, vecType);
      cout << "forming discrete func" << std::endl;
      Expr u0 = new DiscreteFunction(discSpace, 0.5, "u0");

      /* We need a quadrature rule for doing the integrations */
      QuadratureFamily quad = new GaussianQuadrature(2);

      /* Now we set up the weak form of our equation. */
      Expr eqn = Integral(interior, v*(u*u-alpha), quad);

      cout << "equation = " << eqn << std::endl;

      /* There are no boundary conditions for this problem, so the
       * BC expression is empty */
      Expr bc;

      /* We can now set up the nonlinear problem! */

      NonlinearProblem prob(mesh, eqn, bc, v, u, u0, vecType);



#ifdef HAVE_CONFIG_H
      ParameterXMLFileReader reader(searchForFile("SolverParameters/nox.xml"));
#else
      ParameterXMLFileReader reader("nox.xml");
#endif
      ParameterList noxParams = reader.getParameters();

      std::cerr << "solver params = " << noxParams << std::endl;

      NOXSolver solver(noxParams);

      prob.solve(solver);

      /* Inspect solution values. The solution is constant in space,
       * so we can simply take the first NTerms entries in the vector */
      Vector<double> vec = DiscreteFunction::discFunc(u0)->getVector();

      int k=0;
      for (SequentialIterator<double> i=vec.space().begin(); i!=vec.space().end(); i++, k++)
        {
          cout << "u[" << k << "] = " << vec[i] << std::endl;
        }

      double tol = 1.0e-12;
      double errorSq = 0.0;
      Sundance::passFailTest(errorSq, tol);
    }
	catch(std::exception& e)
		{
      std::cerr << e.what() << std::endl;
		}
  Sundance::finalize(); return Sundance::testStatus(); 
}

int main(int argc, char** argv)
{
  try
  {
    const double pi = 4.0*atan(1.0);
    double lambda = 1.25*pi*pi;

    int nx = 32;
    int nt = 10;
    double tFinal = 1.0/lambda;

    Sundance::setOption("nx", nx, "Number of elements");
    Sundance::setOption("nt", nt, "Number of timesteps");
    Sundance::setOption("tFinal", tFinal, "Final time");
    
    Sundance::init(&argc, &argv);

    /* Creation of vector type */
    VectorType<double> vecType = new EpetraVectorType();

    /* Set up mesh */
    MeshType meshType = new BasicSimplicialMeshType();
      
    MeshSource meshSrc = new PartitionedRectangleMesher(
      0.0, 1.0, nx,
      0.0, 1.0, nx,
      meshType);
    Mesh mesh = meshSrc.getMesh();

    /* 
     * Specification of cell filters
     */
    CellFilter interior = new MaximalCellFilter();
    CellFilter edges = new DimensionalCellFilter(1);
    CellFilter west = edges.coordSubset(0, 0.0);
    CellFilter east = edges.coordSubset(0, 1.0);
    CellFilter south = edges.coordSubset(1, 0.0);
    CellFilter north = edges.coordSubset(1, 1.0);

    /* set up test and unknown functions */
    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "u");
    Expr v = new TestFunction(basis, "v");

    /* set up differential operators */
    Expr grad = gradient(2);

    Expr x = new CoordExpr(0);
    Expr y = new CoordExpr(1);

    Expr t = new Sundance::Parameter(0.0);
    Expr tPrev = new Sundance::Parameter(0.0);


    DiscreteSpace discSpace(mesh, basis, vecType);
    Expr uExact = cos(0.5*pi*y)*sin(pi*x)*exp(-lambda*t);
    L2Projector proj(discSpace, uExact);
    Expr uPrev = proj.project();


    /* 
     * We need a quadrature rule for doing the integrations 
     */
    QuadratureFamily quad = new GaussianQuadrature(2);

    double deltaT = tFinal/nt;

    Expr gWest = -pi*exp(-lambda*t)*cos(0.5*pi*y);
    Expr gWestPrev = -pi*exp(-lambda*tPrev)*cos(0.5*pi*y);
    
    /* Create the weak form */
    Expr eqn = Integral(interior, v*(u-uPrev)/deltaT
      + 0.5*(grad*v)*(grad*u + grad*uPrev), quad)
      + Integral(west, -0.5*v*(gWest+gWestPrev), quad);

    Expr bc = EssentialBC(east + north, v*u, quad);

    
    LinearProblem prob(mesh, eqn, bc, v, u, vecType);

    
    LinearSolver<double> solver 
      = LinearSolverBuilder::createSolver("amesos.xml");

    FieldWriter w0 = new VTKWriter("TransientHeat2D-0");
    w0.addMesh(mesh);
    w0.addField("T", new ExprFieldWrapper(uPrev[0]));
    w0.write();

    for (int i=0; i<nt; i++)
    {
      t.setParameterValue((i+1)*deltaT);
      tPrev.setParameterValue(i*deltaT);
      Out::root() << "t=" << (i+1)*deltaT << endl;
      Expr uNext = prob.solve(solver);
      
      ostringstream oss;
      oss << "TransientHeat2D-" << i+1;
      FieldWriter w = new VTKWriter(oss.str());
      w.addMesh(mesh);
//.........这里部分代码省略.........

bool NonlinReducedIntegration()
{
  int np = MPIComm::world().getNProc();

  int n = 4;
  bool increaseProbSize = true;
  if ( (np % 4)==0 ) increaseProbSize = false;

  Array<double> h;
  Array<double> errQuad;
  Array<double> errReduced;

  for (int i=0; i<4; i++)
  {
    n *= 2;
    int nx = n;
    int ny = n;

    VectorType<double> vecType = new EpetraVectorType();

    MeshType meshType = new BasicSimplicialMeshType();
      
    int npx = -1;
    int npy = -1;
    PartitionedRectangleMesher::balanceXY(np, &npx, &npy);
    TEUCHOS_TEST_FOR_EXCEPT(npx < 1);
    TEUCHOS_TEST_FOR_EXCEPT(npy < 1);
    TEUCHOS_TEST_FOR_EXCEPT(npx * npy != np);
    if (increaseProbSize)
    {
      nx = nx*npx;
      ny = ny*npy;
    }
    MeshSource mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, npx, 
      0.0,  1.0, ny, npy, meshType);
    Mesh mesh = mesher.getMesh();


    WatchFlag watchMe("watch eqn");
    watchMe.setParam("integration setup", 0);
    watchMe.setParam("integration", 0);
    watchMe.setParam("fill", 0);
    watchMe.setParam("evaluation", 0);
    watchMe.deactivate();

    WatchFlag watchBC("watch BCs");
    watchBC.setParam("integration setup", 0);
    watchBC.setParam("integration", 0);
    watchBC.setParam("fill", 0);
    watchBC.setParam("evaluation", 0);
    watchBC.deactivate();
    


    CellFilter interior = new MaximalCellFilter();
    CellFilter edges = new DimensionalCellFilter(1);

    CellFilter left = edges.subset(new CoordinateValueCellPredicate(0,0.0));
    CellFilter right = edges.subset(new CoordinateValueCellPredicate(0,1.0));
    CellFilter top = edges.subset(new CoordinateValueCellPredicate(1,1.0));
    CellFilter bottom = edges.subset(new CoordinateValueCellPredicate(1,0.0));

    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "u");
    Expr v = new TestFunction(basis, "v");

    Expr dx = new Derivative(0);
    Expr dy = new Derivative(1);
    Expr grad = List(dx, dy);
    Expr x = new CoordExpr(0);
    Expr y = new CoordExpr(1);

    QuadratureFamily quad = new ReducedQuadrature();
    QuadratureFamily quad2 = new GaussianQuadrature(2);

    /* Define the weak form */
    const double pi = 4.0*atan(1.0);

    Expr c = cos(pi*x);
    Expr s = sin(pi*x);
    Expr ch = cosh(y);
    Expr sh = sinh(y);
    Expr s2 = s*s; 
    Expr c2 = c*c;
    Expr sh2 = sh*sh;
    Expr ch2 = ch*ch;
    Expr pi2 = pi*pi;
    Expr uEx = s*ch;
    Expr eu = exp(uEx);
    Expr f = -(ch*eu*(-1 + pi2)*s) + ch2*(c2*eu*pi2 - s2) + eu*s2*sh2;

    Expr eqn = Integral(interior, exp(u)*(grad*u)*(grad*v)
      + v*f + v*u*u, quad, watchMe)
      + Integral(right, v*exp(u)*pi*cosh(y), quad,watchBC);
    /* Define the Dirichlet BC */
    Expr bc = EssentialBC(left+top, v*(u-uEx), quad, watchBC);

    Expr eqn2 = Integral(interior, exp(u)*(grad*u)*(grad*v)
      + v*f + v*u*u, quad2, watchMe)
      + Integral(right, v*exp(u)*pi*cosh(y), quad2,watchBC);
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  try
		{
			Sundance::init(&argc, &argv);
      int np = MPIComm::world().getNProc();
      
      int nx = 48;
      int ny = 48;
      int npx = -1;
      int npy = -1;
      PartitionedRectangleMesher::balanceXY(np, &npx, &npy);
      TEUCHOS_TEST_FOR_EXCEPT(npx < 1);
      TEUCHOS_TEST_FOR_EXCEPT(npy < 1);
      TEUCHOS_TEST_FOR_EXCEPT(npx * npy != np);
      MeshType meshType = new BasicSimplicialMeshType();
      MeshSource mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, npx, 
        0.0,  1.0, ny, npy, meshType);

      Mesh mesh = mesher.getMesh();
      CellFilter interior = new MaximalCellFilter();
      CellFilter bdry = new BoundaryCellFilter();
      
      /* Create a vector space factory, used to 
       * specify the low-level linear algebra representation */
      VectorType<double> vecType = new EpetraVectorType();
  
      /* create a discrete space on the mesh */
      DiscreteSpace discreteSpace(mesh, new Lagrange(1), vecType);

      /* initialize the design, state, and multiplier vectors */
      Expr alpha0 = new DiscreteFunction(discreteSpace, 1.0, "alpha0");
      Expr u0 = new DiscreteFunction(discreteSpace, 1.0, "u0");
      Expr lambda0 = new DiscreteFunction(discreteSpace, 1.0, "lambda0");

      /* create symbolic objects for test and unknown functions */
      Expr u = new UnknownFunction(new Lagrange(1), "u");
      Expr lambda = new UnknownFunction(new Lagrange(1), "lambda");
      Expr alpha = new UnknownFunction(new Lagrange(1), "alpha");

      /* create symbolic differential operators */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr grad = List(dx, dy);

      /* create symbolic coordinate functions */
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);

      /* create target function */
      const double pi = 4.0*atan(1.0);
      Expr uStar = sin(pi*x)*sin(pi*y);
      
      /* create quadrature rules of different orders */
      QuadratureFamily q1 = new GaussianQuadrature(1);
      QuadratureFamily q2 = new GaussianQuadrature(2);
      QuadratureFamily q4 = new GaussianQuadrature(4);

      /* Regularization weight */
      double R = 0.001;
      double U0 = 1.0/(1.0 + 4.0*pow(pi,4.0)*R);
      double A0 = -2.0*pi*pi*U0;

      /* Form objective function */
      Expr reg = Integral(interior, 0.5 * R * alpha*alpha, q2);
      Expr fit = Integral(interior, 0.5 * pow(u-uStar, 2.0), q4);

      Expr constraintEqn = Integral(interior, 
        (grad*lambda)*(grad*u) + lambda*alpha, q2);
      Expr L = reg + fit + constraintEqn;

      Expr constraintBC = EssentialBC(bdry, lambda*u, q2);
      Functional Lagrangian(mesh, L, constraintBC, vecType);
      
      LinearSolver<double> solver 
        = LinearSolverBuilder::createSolver("amesos.xml");

      RCP<ObjectiveBase> obj = rcp(new LinearPDEConstrainedObj(
        Lagrangian, u, u0, lambda, lambda0, alpha, alpha0,
        solver));

      Vector<double> xInit = obj->getInit();

      bool doFDCheck = false;
      if (doFDCheck)
      {
        Out::root() << "Doing FD check of gradient..." << endl;
        bool fdOK = obj->fdCheck(xInit, 1.0e-6, 2);
        if (fdOK) 
        {
          Out::root() << "FD check OK" << endl;
        }
        else
        {
          Out::root() << "FD check FAILED" << endl;
        }
      }

      RCP<UnconstrainedOptimizerBase> opt 
          = OptBuilder::createOptimizer("basicLMBFGS.xml");
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  int stat = 0;
  int verb=1;
  try
  {
    GlobalMPISession session(&argc, &argv);

    TimeMonitor t(totalTimer());

    int pMax = 2;
    int dim=2;

    bool isInternalBdry = false;

    Utils::setChopVal(1.0e-14);

    CellType cellType = TriangleCell;

    //       Point a = Point(1.0, 1.0);
    //       Point b = Point(1.2, 1.6);
    //       Point c = Point(0.8, 1.3);

    Point a = Point(0.0, 0.0);
    Point b = Point(1.0, 0.0);
    Point c = Point(0.0, 1.0);

    Point d = Point(0.1, 0.1);
    Point e = Point(1.0, 0.0);
    Point f = Point(0.0, 1.0);

    int nCells = 2;

    CellJacobianBatch JBatch;
    JBatch.resize(nCells, 2, 2);
    double* J = JBatch.jVals(0);
    J[0] = b[0] - a[0];
    J[1] = c[0] - a[0];
    J[2] = b[1] - a[1];
    J[3] = c[1] - a[1];

    J[4] = e[0] - d[0];
    J[5] = f[0] - d[0];
    J[6] = e[1] - d[1];
    J[7] = f[1] - d[1];


      
    Array<int> dummy;
    double coeff = 1.0;
    RCP<Array<double> > A = rcp(new Array<double>());
    RCP<Array<double> > B = rcp(new Array<double>());

    QuadratureFamily q4 = new GaussianQuadrature(4);

    int nErrors = 0;

    std::cerr << std::endl << std::endl 
         << "---------------- One-forms --------------------" 
         << std::endl << std::endl;
    for (int p=0; p<=pMax; p++)
    {
      BasisFamily P = new Lagrange(p);
      for (int dp=0; dp<=1; dp++)
      {
        if (dp > p) continue;

        int numTestDir = 1;
        if (dp==1) numTestDir = dim;
        for (int t=0; t<numTestDir; t++)
        {
          int alpha = t;
          Tabs tab;

          ParametrizedCurve curve = new DummyParametrizedCurve();
          MeshType meshType = new BasicSimplicialMeshType();
          MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, 10, meshType);
          Mesh mesh = mesher.getMesh();
          RCP<Array<int> > cellLIDs;

          RefIntegral ref(dim, cellType, dim, cellType, P, alpha, dp, q4 , isInternalBdry, curve, mesh ,verb);
          A->resize(JBatch.numCells() * ref.nNodes());
          for (int ai=0; ai<A->size(); ai++) (*A)[ai]=0.0;
          ref.transformOneForm(JBatch, JBatch, dummy, cellLIDs , coeff, A);
          std::cerr << tab << "transformed reference element" << std::endl;
          if (dp>0) std::cerr << tab << "test diff direction=" << t << std::endl;
          for (int cell=0; cell<nCells; cell++)
          {
            std::cerr << tab << "{";
            for (int r=0; r<ref.nNodesTest(); r++)
            {
              if (r!=0) std::cerr << ", ";
              std::cerr << Utils::chop((*A)[cell*ref.nNodesTest()+r]);
            }
            std::cerr << "}" << std::endl;
          }
          QuadratureIntegral quad(dim, cellType, dim, cellType, P, alpha, dp, q4, isInternalBdry, curve, mesh, verb);
          Array<double> quadCoeff(2*quad.nQuad(), 1.0);
          B->resize(JBatch.numCells() * quad.nNodes());
          for (int ai=0; ai<B->size(); ai++) (*B)[ai]=0.0;
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  
  try
		{
      Sundance::init(&argc, &argv);
      int np = MPIComm::world().getNProc();

      //  DOFMapBase::classVerbosity() = VerbExtreme;

      const double density = 1000.0; // kg/m^3
      const double porosity = 0.442; // dimensionless %
      const double A = 175.5; // dimensionless fit parameter
      const double B = 1.83;  // dimensionless fit parameter
      const double criticalRe = 36.73;  // dimensionless fit parameter
      const double dynvisc = 1.31;  // kg/(m-s)
      const double graindia = 1.9996e-4;  // m 
      const double charvel = 1.0;  // m/s

      
      double Reynolds = density*graindia*charvel/(dynvisc*porosity);

      Expr Re = new Parameter(Reynolds);

      /* We will do our linear algebra using Epetra */
      VectorType<double> vecType = new EpetraVectorType();

      /* Create a mesh. It will be of type BasisSimplicialMesh, and will
       * be built using a PartitionedLineMesher. */
      MeshType meshType = new BasicSimplicialMeshType();
      MeshSource mesher = new PartitionedLineMesher(0.0, 1000.0, 100*np, meshType);
      Mesh mesh = mesher.getMesh();

      /* Create a cell filter that will identify the maximal cells
       * in the interior of the domain */
      CellFilter interior = new MaximalCellFilter();
      CellFilter points = new DimensionalCellFilter(0);
      CellPredicate leftPointFunc = new PositionalCellPredicate(leftPointTest);
      CellPredicate rightPointFunc = new PositionalCellPredicate(rightPointTest);
      CellFilter leftPoint = points.subset(leftPointFunc);
      CellFilter rightPoint = points.subset(rightPointFunc);
      
      /* Create unknown and test functions, discretized using first-order
       * Lagrange interpolants */
      
      Expr p = new UnknownFunction(new Lagrange(2), "p");
      Expr q = new UnknownFunction(new Lagrange(2), "q");
 
      Expr u = new TestFunction(new Lagrange(2), "u");
      Expr v = new TestFunction(new Lagrange(2), "v");

      /* Create differential operator and coordinate function */
      Expr dx = new Derivative(0);
      Expr x = new CoordExpr(0);

      /* We need a quadrature rule for doing the integrations */
      QuadratureFamily quad = new GaussianQuadrature(4);

      /* Define the weak form */
      Expr MassEqn = Integral(interior, q*(dx*u), quad)
	+ Integral(leftPoint, - q*u,quad)
	+ Integral(rightPoint,  - q*u,quad);
      Expr MomEqn = Integral(interior, (density/porosity)*q*q*(dx*v) + porosity*p*(dx*v) - porosity*q*v*A - (porosity*q*v*B*Re*Re)/((Re+criticalRe)*(1-porosity)), quad)
	+ Integral(leftPoint, - density*q*q*v/porosity - porosity*p*v,quad)
	+ Integral(rightPoint,- density*q*q*v/porosity - porosity*p*v,quad);

      /* Define the Dirichlet BC */
      Expr leftbc = EssentialBC(leftPoint, v*(q-charvel), quad);
      Expr rightbc = EssentialBC(rightPoint, v*(q-charvel), quad);

      /* Create a discrete space, and discretize the function 1.0 on it */
      BasisFamily L2 = new Lagrange(2);
      Array<BasisFamily> basis = tuple(L2, L2);
      DiscreteSpace discSpace(mesh, basis, vecType);
      Expr u0 = new DiscreteFunction(discSpace, 1.0, "u0");
      Expr p0 = u0[0];
      Expr q0 = u0[1];
     
 
/* Create a TSF NonlinearOperator object */
      std::cerr << "about to make nonlinear object" << std::endl;
      std::cerr.flush();

      NonlinearOperator<double> F 
        = new NonlinearProblem(mesh, MassEqn+MomEqn, leftbc+rightbc, Sundance::List(u,v),Sundance::List(p,q) , u0, vecType);
    
      //      F.verbosity() = VerbExtreme;
      /* Get the initial guess */
  
      Vector<double> x0 = F.getInitialGuess();
   
      
      /* Create an Aztec solver for solving the linear subproblems */
      std::map<int,int> azOptions;
      std::map<int,double> azParams;
      
      azOptions[AZ_solver] = AZ_gmres;
      azOptions[AZ_precond] = AZ_dom_decomp;
      azOptions[AZ_subdomain_solve] = AZ_ilu;
      azOptions[AZ_graph_fill] = 1;
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  try
  {
    int nx = 32;
    double convTol = 1.0e-8;
    double lambda = 0.5;
    Sundance::setOption("nx", nx, "Number of elements");
    Sundance::setOption("tol", convTol, "Convergence tolerance");
    Sundance::setOption("lambda", lambda, "Lambda (parameter in Bratu's equation)");

    Sundance::init(&argc, &argv);

    Out::root() << "Bratu problem (lambda=" << lambda << ")" << endl;
    Out::root() << "Newton's method, linearized by hand" << endl << endl;

    VectorType<double> vecType = new EpetraVectorType();

    MeshType meshType = new BasicSimplicialMeshType();
    MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
    Mesh mesh = mesher.getMesh();

    CellFilter interior = new MaximalCellFilter();
    CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
    CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
    CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
    
    BasisFamily basis = new Lagrange(1);
    Expr w = new UnknownFunction(basis, "w");
    Expr v = new TestFunction(basis, "v");

    Expr grad = gradient(1);

    Expr x = new CoordExpr(0);



    const double pi = 4.0*atan(1.0);
    Expr uExact = sin(pi*x);
    Expr R = pi*pi*uExact - lambda*exp(uExact);

    QuadratureFamily quad4 = new GaussianQuadrature(4);
    QuadratureFamily quad2 = new GaussianQuadrature(2);

    DiscreteSpace discSpace(mesh, basis, vecType);
    Expr uPrev = new DiscreteFunction(discSpace, 0.5);
    Expr stepVal = copyDiscreteFunction(uPrev);

    Expr eqn 
      = Integral(interior, (grad*v)*(grad*w) + (grad*v)*(grad*uPrev) 
        - v*lambda*exp(uPrev)*(1.0+w) - v*R, quad4);

    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(left+right, v*(uPrev+w)/h, quad2); 

    LinearProblem prob(mesh, eqn, bc, v, w, vecType);

    LinearSolver<double> linSolver 
      = LinearSolverBuilder::createSolver("amesos.xml");

    Out::root() << "Newton iteration" << endl;
    int maxIters = 20;
    Expr soln ;
    bool converged = false;

    for (int i=0; i<maxIters; i++)
    {
      /* solve for the next u */
      prob.solve(linSolver, stepVal);
      Vector<double> stepVec = getDiscreteFunctionVector(stepVal);
      double deltaU = stepVec.norm2();
      Out::root() << "Iter=" << setw(3) << i << " ||Delta u||=" << setw(20)
                  << deltaU << endl;
      addVecToDiscreteFunction(uPrev, stepVec);
      if (deltaU < convTol) 
      {
        soln = uPrev;
        converged = true;
        break;
      }
    } 
    TEUCHOS_TEST_FOR_EXCEPTION(!converged, std::runtime_error, 
      "Newton iteration did not converge after " 
      << maxIters << " iterations");
    
    FieldWriter writer = new DSVWriter("HandCodedBratu.dat");
    writer.addMesh(mesh);
    writer.addField("soln", new ExprFieldWrapper(soln[0]));
    writer.write();

    Out::root() << "Converged!" << endl << endl;

    double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
    Out::root() << "L2 Norm of error: " << L2Err << endl;
    
    Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
  }
	catch(std::exception& e) 
  {
    Sundance::handleException(e);
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  
  try
		{
      Sundance::init(&argc, &argv);
      int np = MPIComm::world().getNProc();
      TEST_FOR_EXCEPT(np != 1);

      /* We will do our linear algebra using Epetra */
      VectorType<double> vecType = new EpetraVectorType();

      /* Create a periodic mesh */
      int nx = 1000;
      const double pi = 4.0*atan(1.0);
      MeshType meshType = new PeriodicMeshType1D();
      MeshSource mesher = new PeriodicLineMesher(0.0, 2.0*pi, nx, meshType);
      Mesh mesh = mesher.getMesh();

      /* Create a cell filter that will identify the maximal cells
       * in the interior of the domain */
      CellFilter interior = new MaximalCellFilter();
      
      /* Create unknown and test functions, discretized using first-order
       * Lagrange interpolants */
      Expr u = new UnknownFunction(new Lagrange(1), "u");
      Expr v = new TestFunction(new Lagrange(1), "v");

      /* Create differential operator and coordinate function */
      Expr dx = new Derivative(0);
      Expr x = new CoordExpr(0);

      /* We need a quadrature rule for doing the integrations */
      QuadratureFamily quad = new GaussianQuadrature(4);

      
      /* Define the weak form */
      Expr eqn = Integral(interior, 
        (dx*v)*(dx*u) - 2.0*v*(dx*u) - v*u + v*sin(2*x),
                          quad);
      Expr bc ; // no explicit BC needed

      /* We can now set up the linear problem! */

      LinearProblem prob(mesh, eqn, bc, v, u, vecType);


      ParameterXMLFileReader reader("amesos.xml");
      ParameterList solverParams = reader.getParameters();

      LinearSolver<double> solver 
        = LinearSolverBuilder::createSolver(solverParams);


      Out::os() << "solving problem " << std::endl;
      Expr soln = prob.solve(solver);

      Expr uExact = -1.0/25.0 * (4.0*cos(2.0*x) + 3.0*sin(2.0*x));

      Expr uErr = uExact - soln;
      
      Expr uErrExpr = Integral(interior, 
                              uErr*uErr,
                              new GaussianQuadrature(6));
      
      FunctionalEvaluator uErrInt(mesh, uErrExpr);

      double uErrorSq = uErrInt.evaluate();
      std::cerr << "u error norm = " << sqrt(uErrorSq) << std::endl << std::endl;


      /* make sure the unfolded solution is also correct */

      Out::os() << "unfolding " << std::endl;
      Expr unfoldedSoln = unfoldPeriodicDiscreteFunction(soln);
      

      Expr ufErr = uExact - unfoldedSoln;
      
      Expr ufErrExpr = Integral(interior, 
                              ufErr*ufErr,
                              new GaussianQuadrature(6));
      
      Mesh unfoldedMesh = DiscreteFunction::discFunc(unfoldedSoln)->mesh();
      FunctionalEvaluator ufErrInt(unfoldedMesh, ufErrExpr);

      double ufErrorSq = ufErrInt.evaluate();
      std::cerr << "unfolded error norm = " << sqrt(ufErrorSq) << std::endl << std::endl;

      double tol = 1.0e-3;
      Sundance::passFailTest(sqrt(uErrorSq + ufErrorSq), tol);

    }
	catch(std::exception& e)
		{
      Sundance::handleException(e);
		}
  Sundance::finalize(); return Sundance::testStatus(); 
}

bool DuffingFloquet()
{
  int np = MPIComm::world().getNProc();
  TEUCHOS_TEST_FOR_EXCEPT(np != 1);

  const double pi = 4.0*atan(1.0);

  /* We will do our linear algebra using Epetra */
  VectorType<double> vecType = new EpetraVectorType();

  /* Create a periodic mesh */
  int nx = 128;

  MeshType meshType = new PeriodicMeshType1D();
  MeshSource mesher = new PeriodicLineMesher(0.0, 2.0*pi, nx, meshType);
  Mesh mesh = mesher.getMesh();

  /* Create a cell filter that will identify the maximal cells
   * in the interior of the domain */
  CellFilter interior = new MaximalCellFilter();
  CellFilter pts = new DimensionalCellFilter(0);
      
  CellFilter left = pts.subset(new CoordinateValueCellPredicate(0,0.0));
  CellFilter right = pts.subset(new CoordinateValueCellPredicate(0,2.0*pi));
      
  /* Create unknown and test functions, discretized using first-order
   * Lagrange interpolants */
  Expr u1 = new UnknownFunction(new Lagrange(1), "u1");
  Expr u2 = new UnknownFunction(new Lagrange(1), "u2");
  Expr v1 = new TestFunction(new Lagrange(1), "v1");
  Expr v2 = new TestFunction(new Lagrange(1), "v2");

  /* Create differential operator and coordinate function */
  Expr dx = new Derivative(0);
  Expr x = new CoordExpr(0);

  /* We need a quadrature rule for doing the integrations */
  QuadratureFamily quad = new GaussianQuadrature(4);

  double F0 = 0.5;
  double gamma = 2.0/3.0;
  double a0 = 1.0;
  double w0 = 1.0;
  double eps = 0.5;

  Expr u1Guess = -0.75*cos(x) + 0.237*sin(x);
  Expr u2Guess = 0.237*cos(x) + 0.75*sin(x);

  DiscreteSpace discSpace(mesh, 
    List(new Lagrange(1), new Lagrange(1)),
    vecType);
  L2Projector proj(discSpace, List(u1Guess, u2Guess));
  Expr u0 = proj.project();


  Expr rhs1 = u2;
  Expr rhs2 = -w0*w0*u1 - gamma*u2 - eps*w0*w0*pow(u1,3.0)/a0/a0 
    + F0*w0*w0*sin(x);

  /* Define the weak form */
  Expr eqn = Integral(interior, 
    v1*(dx*u1 - rhs1) + v2*(dx*u2 - rhs2),
    quad);
  Expr dummyBC ; 

  NonlinearProblem prob(mesh, eqn, dummyBC, List(v1,v2), List(u1,u2), 
    u0, vecType);


  ParameterXMLFileReader reader("nox.xml");
  ParameterList solverParams = reader.getParameters();

  Out::root() << "finding periodic solution" << endl;
  NOXSolver solver(solverParams);
  prob.solve(solver);

  /* unfold the solution onto a non-periodic mesh */
      
  Expr uP = unfoldPeriodicDiscreteFunction(u0, "u_p");
  Out::root() << "uP=" << uP << endl;
      
  Mesh unfoldedMesh = DiscreteFunction::discFunc(uP)->mesh();
  DiscreteSpace unfDiscSpace = DiscreteFunction::discFunc(uP)->discreteSpace();

  FieldWriter writer = new MatlabWriter("Floquet.dat");
  writer.addMesh(unfoldedMesh);
  writer.addField("u_p[0]", new ExprFieldWrapper(uP[0]));
  writer.addField("u_p[1]", new ExprFieldWrapper(uP[1]));

  Array<Expr> a(2);
  a[0] = new Sundance::Parameter(0.0, "a1");
  a[1] = new Sundance::Parameter(0.0, "a2");


  Expr bc = EssentialBC(left, v1*(u1-uP[0]-a[0]) + v2*(u2-uP[1]-a[1]), quad);

  NonlinearProblem unfProb(unfoldedMesh, eqn, bc, 
    List(v1,v2), List(u1,u2), uP, vecType);

  unfProb.setEvalPoint(uP);
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  try
  {
    int nx = 32;
    double convTol = 1.0e-8;
    double lambda = 0.5;
    Sundance::setOption("nx", nx, "Number of elements");
    Sundance::setOption("tol", convTol, "Convergence tolerance");
    Sundance::setOption("lambda", lambda, "Lambda (parameter in Bratu's equation)");

    Sundance::init(&argc, &argv);

    Out::root() << "Bratu problem (lambda=" << lambda << ")" << endl;
    Out::root() << "Fixed-point iteration" << endl << endl;

    VectorType<double> vecType = new EpetraVectorType();

    MeshType meshType = new BasicSimplicialMeshType();
    MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
    Mesh mesh = mesher.getMesh();

    CellFilter interior = new MaximalCellFilter();
    CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
    CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
    CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
    
    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "u");
    Expr v = new TestFunction(basis, "v");

    Expr grad = gradient(1);

    Expr x = new CoordExpr(0);



    const double pi = 4.0*atan(1.0);
    Expr uExact = sin(pi*x);
    Expr R = pi*pi*uExact - lambda*exp(uExact);

    QuadratureFamily quad4 = new GaussianQuadrature(4);
    QuadratureFamily quad2 = new GaussianQuadrature(2);

    DiscreteSpace discSpace(mesh, basis, vecType);
    Expr uPrev = new DiscreteFunction(discSpace, 0.5);
    Expr uCur = copyDiscreteFunction(uPrev);

    Expr eqn 
      = Integral(interior, (grad*u)*(grad*v) - v*lambda*exp(uPrev) - v*R, quad4);

    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(left+right, v*u/h, quad4); 

    LinearProblem prob(mesh, eqn, bc, v, u, vecType);

    Expr normSqExpr = Integral(interior, pow(u-uPrev, 2.0), quad2);
    Functional normSqFunc(mesh, normSqExpr, vecType);
    FunctionalEvaluator normSqEval = normSqFunc.evaluator(u, uCur);

    LinearSolver<double> linSolver 
      = LinearSolverBuilder::createSolver("amesos.xml");

    Out::root() << "Fixed-point iteration" << endl;
    int maxIters = 20;
    Expr soln ;
    bool converged = false;

    for (int i=0; i<maxIters; i++)
    {
      /* solve for the next u */
      prob.solve(linSolver, uCur);
      /* evaluate the norm of (uCur-uPrev) using 
       * the FunctionalEvaluator defined above */
      double deltaU = sqrt(normSqEval.evaluate());
      Out::root() << "Iter=" << setw(3) << i << " ||Delta u||=" << setw(20)
                  << deltaU << endl; 
      /* check for convergence */  
      if (deltaU < convTol) 
      {
        soln = uCur;
        converged = true;
        break;
      }
      /* get the vector from the current discrete function */
      Vector<double> uVec = getDiscreteFunctionVector(uCur);
      /* copy the vector into the previous discrete function */ 
      setDiscreteFunctionVector(uPrev, uVec);
    } 
    TEUCHOS_TEST_FOR_EXCEPTION(!converged, std::runtime_error, 
      "Fixed point iteration did not converge after " 
      << maxIters << " iterations");
    
    FieldWriter writer = new DSVWriter("FixedPointBratu.dat");
    writer.addMesh(mesh);
    writer.addField("soln", new ExprFieldWrapper(soln[0]));
    writer.write();

    Out::root() << "Converged!" << endl << endl;

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

int main(int argc, char** argv)
{
  try
  {
    int nx = 32;
    double convTol = 1.0e-8;
    double lambda = 0.5;
    Sundance::setOption("nx", nx, "Number of elements");
    Sundance::setOption("tol", convTol, "Convergence tolerance");
    Sundance::setOption("lambda", lambda, "Lambda (parameter in Bratu's equation)");

    Sundance::init(&argc, &argv);

    Out::root() << "Bratu problem (lambda=" << lambda << ")" << endl;
    Out::root() << "Newton's method with automated linearization" 
                << endl << endl;

    VectorType<double> vecType = new EpetraVectorType();

    MeshType meshType = new BasicSimplicialMeshType();
    MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, meshType);
    Mesh mesh = mesher.getMesh();

    CellFilter interior = new MaximalCellFilter();
    CellFilter sides = new DimensionalCellFilter(mesh.spatialDim()-1);
    CellFilter left = sides.subset(new CoordinateValueCellPredicate(0, 0.0));
    CellFilter right = sides.subset(new CoordinateValueCellPredicate(0, 1.0));
    
    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "w");
    Expr v = new TestFunction(basis, "v");

    Expr grad = gradient(1);

    Expr x = new CoordExpr(0);

    const double pi = 4.0*atan(1.0);
    Expr uExact = sin(pi*x);
    Expr R = pi*pi*uExact - lambda*exp(uExact);

    QuadratureFamily quad4 = new GaussianQuadrature(4);
    QuadratureFamily quad2 = new GaussianQuadrature(2);

    DiscreteSpace discSpace(mesh, basis, vecType);
    Expr uPrev = new DiscreteFunction(discSpace, 0.5);

    Expr eqn 
      = Integral(interior, (grad*v)*(grad*u) - v*lambda*exp(u) - v*R, quad4);

    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(left+right, v*u/h, quad2); 

    NonlinearProblem prob(mesh, eqn, bc, v, u, uPrev, vecType);

    NonlinearSolver<double> solver 
      = NonlinearSolverBuilder::createSolver("playa-newton-amesos.xml");

    Out::root() << "Newton solve" << endl;

    SolverState<double> state = prob.solve(solver);
    
    TEUCHOS_TEST_FOR_EXCEPTION(state.finalState() != SolveConverged,
      std::runtime_error,
      "Nonlinear solve failed to converge: message=" << state.finalMsg());
    
    Expr soln = uPrev;
    FieldWriter writer = new DSVWriter("AutoLinearizedBratu.dat");
    writer.addMesh(mesh);
    writer.addField("soln", new ExprFieldWrapper(soln[0]));
    writer.write();

    Out::root() << "Converged!" << endl << endl;

    double L2Err = L2Norm(mesh, interior, soln-uExact, quad4);
    Out::root() << "L2 Norm of error: " << L2Err << endl;
    
    Sundance::passFailTest(L2Err, 1.5/((double) nx*nx));
  }
	catch(std::exception& e) 
  {
    Sundance::handleException(e);
  }
  Sundance::finalize(); 
  return Sundance::testStatus();
}

bool AToCDensitySample()
{

  /* We will do our linear algebra using Epetra */
  VectorType<double> vecType = new EpetraVectorType();

  /* Create a mesh. It will be of type BasisSimplicialMesh, and will
   * be built using a PartitionedLineMesher. */
  MeshType meshType = new BasicSimplicialMeshType();
  MeshSource mesher = new PartitionedRectangleMesher(-1.0, 1.0, 32, 1,
    -1.0, 1.0, 32, 1,
    meshType);
  Mesh mesh = mesher.getMesh();

  /* Create a cell filter that will identify the maximal cells
   * in the interior of the domain */
  CellFilter interior = new MaximalCellFilter();

  Expr x = new CoordExpr(0);
  Expr y = new CoordExpr(1);

  BasisFamily L1 = new Lagrange(1);
  DiscreteSpace discSpace(mesh, List(L1, L1), vecType);
      
  /* Discretize some expression for the force. We'll pick a linear function
   * so that it can be interpolated exactly, letting us check the 
   * validity of our interpolations. */
  L2Projector proj(discSpace, List(x, y));
  Expr F = proj.project();

      
  /* create a sampler */
  cout << "making grid" << std::endl;
  AToCPointLocator locator(mesh, interior, createVector(tuple(200, 200)));
      
  AToCDensitySampler sampler(locator, vecType);

  CToAInterpolator forceInterpolator(locator, F);

  cout << "making points" << std::endl;
  /* create a bunch of particles */
  int nCells = mesh.numCells(2);
  int nPts = 15000;

  Array<double> pos(2*nPts);
  Array<double> f(F.size() * nPts);
  Array<Point> physPts;

  /* We'll generate random sample points in a way that lets us make an exact check
   * of the density recovery. We pick random cells, then random local coordinates
   * within each cell. This way, we can compute the density exactly as we
   * go, giving us something to check the recovered density against. */
  Array<int> counts(nCells);
  for (int i=0; i<nPts; i++)
  {
    /* pick a random cell */
    int cell = (int) floor(nCells * drand48());
    counts[cell]++;
    /* generate a point in local coordinates */
    double s = drand48();
    double t = drand48() * (1.0-s);
    Point refPt(s, t);
    /* map to physical coordinates */
    mesh.pushForward(2, tuple(cell), tuple(refPt), physPts);
    Point X = physPts[0];
    pos[2*i] = X[0];
    pos[2*i+1] = X[1];
  }

  cout << "sampling..." << std::endl;
  Expr density = sampler.sample(createVector(pos), 1.0);

  cout << "computing forces..." << std::endl;
  forceInterpolator.interpolate(pos, f);

  double maxForceErr = 0.0;
  for (int i=0; i<nPts; i++)
  {
    double x0 = pos[2*i];
    double y0 = pos[2*i+1];
    double fx = x0;
    double fy = y0;
    double df = ::fabs(fx - f[2*i]) + ::fabs(fy - f[2*i+1]);
    maxForceErr = max(maxForceErr, df);
  }
  cout << "max force error = " << maxForceErr << std::endl;

  cout << "writing..." << std::endl;

  /* Write the field in VTK format */
  FieldWriter w = new VTKWriter("Density2d");
  w.addMesh(mesh);
  w.addField("rho", new ExprFieldWrapper(density));
  w.write();

  double errorSq = 0.0;

  double tol = 1.0e-6;
  return SundanceGlobal::passFailTest(::sqrt(errorSq), tol);
}

bool BlockStochPoissonTest1D()
{
  /* We will do our linear algebra using Epetra */
  VectorType<double> vecType = new EpetraVectorType();

  /* Read a mesh */
  MeshType meshType = new BasicSimplicialMeshType();
  int nx = 32;
  MeshSource mesher = new PartitionedLineMesher(0.0, 1.0, nx, 
    meshType);
  Mesh mesh = mesher.getMesh();

  /* Create a cell filter that will identify the maximal cells
   * in the interior of the domain */
  CellFilter interior = new MaximalCellFilter();
  CellFilter pts = new DimensionalCellFilter(0);
  CellFilter left = pts.subset(new CoordinateValueCellPredicate(0,0.0));
  CellFilter right = pts.subset(new CoordinateValueCellPredicate(0,1.0));

  Expr x = new CoordExpr(0);

  /* Create the stochastic coefficients */
  int nDim = 1;
  int order = 6;
#ifdef HAVE_SUNDANCE_STOKHOS
  Out::root() << "using Stokhos hermite basis" << std::endl;
  SpectralBasis pcBasis = new Stokhos::HermiteBasis<int,double>(order);
#else
  Out::root() << "using George's hermite basis" << std::endl;
  SpectralBasis pcBasis = new HermiteSpectralBasis(nDim, order);
#endif
    
  Array<Expr> q(pcBasis.nterms());
  Array<Expr> kappa(pcBasis.nterms());
  Array<Expr> uEx(pcBasis.nterms());

  double a = 0.1;

  q[0] = -2 + pow(a,2)*(4 - 9*x)*x - 2*pow(a,3)*(-1 + x)*(1 + 3*x*(-3 + 4*x));
  q[1] = -(a*(-3 + 10*x + 2*a*(-1 + x*(8 - 9*x +
          a*(-4 + 3*(5 - 4*x)*x + 12*a*(-1 + x)*(1 + 5*(-1 + x)*x))))));
  q[2] = a*(-4 + 6*x + a*(1 - x*(2 + 3*x) + a*(4 - 28*x + 30*pow(x,2))));
  q[3] = -(pow(a,2)*(-3 + x*(20 - 21*x +
        a*(-4 + 3*(5 - 4*x)*x + 24*a*(-1 + x)*(1 + 5*(-1 + x)*x)))));
  q[4] = pow(a,3)*(1 + x*(-6 + x*(3 + 4*x)));
  q[5] = -4*pow(a,4)*(-1 + x)*x*(1 + 5*(-1 + x)*x);
  q[6] = 0.0;

  uEx[0] = -((-1 + x)*x);
  uEx[1] = -(a*(-1 + x)*pow(x,2));
  uEx[2] = a*pow(-1 + x,2)*x;
  uEx[3] = pow(a,2)*pow(-1 + x,2)*pow(x,2);
  uEx[4] = 0.0;
  uEx[5] = 0.0;
  uEx[6] = 0.0;

  kappa[0] = 1.0;
  kappa[1] = a*x;
  kappa[2] = -(pow(a,2)*(-1 + x)*x);

  kappa[3] = 1.0; // unused
  kappa[4] = 1.0; // unused
  kappa[5] = 1.0; // unused
  kappa[6] = 1.0; // unused


  Array<Expr> uBC(pcBasis.nterms());
  for (int i=0; i<pcBasis.nterms(); i++) uBC[i] = 0.0;

  int L = nDim+2;
  int P = pcBasis.nterms();
  Out::os() << "L = " << L << std::endl;
  Out::os() << "P = " << P << std::endl;
    
  /* Create the unknown and test functions. Do NOT use the spectral
   * basis here */
  Expr u = new UnknownFunction(new Lagrange(4), "u");
  Expr v = new TestFunction(new Lagrange(4), "v");

  /* Create differential operator and coordinate function */
  Expr dx = new Derivative(0);
  Expr grad = dx;


  /* We need a quadrature rule for doing the integrations */
  QuadratureFamily quad = new GaussianQuadrature(12);

  /* Now we create problem objects to build each $K_j$ and $f_j$.
   * There will be L matrix-vector pairs */
  Array<Expr> eqn(P);
  Array<Expr> bc(P);
  Array<LinearProblem> prob(P);
  Array<LinearOperator<double> > KBlock(L);
  Array<Vector<double> > fBlock(P);
  Array<Vector<double> > solnBlock;

  for (int j=0; j<P; j++)
  {
    eqn[j] = Integral(interior, kappa[j]*(grad*v)*(grad*u) + v*q[j], quad);
    bc[j] = EssentialBC(left+right, v*(u-uBC[j]), quad);
//.........这里部分代码省略.........