C++ MeshSource类代码示例

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();
}

	void MeshGroup::_build(int index)
	{
		const MultiMap<MeshSource *, Mesh *>::Node * i = mMeshBucket[index].Begin();
		while (i != mMeshBucket[index].End())
		{
			int iVertexCount = 0;
			int iIndexCount = 0;

			MeshSource * source = i->key;
			for (int j = 0; j < source->GetMeshBufferCount(); ++j)
			{
				MeshBuffer * mb = source->GetMeshBuffer(j);
				RenderOp * pRenderOp = mb->GetRenderOp();

				d_assert (
					!pRenderOp->vertexDeclarations[0].HasElement(eVertexSemantic::POSITION, eVertexType::FLOAT3) &&
					pRenderOp->vertexBuffers[0] != NULL &&
					pRenderOp->indexBuffer != NULL &&
					pRenderOp->primType != ePrimType::TRIANGLE_LIST);

				iVertexCount += pRenderOp->vertexBuffers[0]->GetCount();
				iIndexCount += pRenderOp->indexBuffer->GetCount();
			}

			int first = 0, last = 0;
			int vcount = 0, icount = 0;
			while (last < i->members.Size())
			{
				vcount += iVertexCount;
				icount += iIndexCount;
				last += 1;

				if (last == i->members.Size() ||
					iVertexCount > 4096 ||
					iIndexCount > 4096 * 3)
				{
					_genMesh(i->members, first, last, index != 0);

					first = last;
				}
			}

			i = i->Next();
		}
	}

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 {
        int nx = 128;
        double C = 4.0;
        std::string solverFile = "aztec-ml.xml";
        Sundance::setOption("nx", nx, "number of elements in x");
        Sundance::setOption("C", C, "Nitsche penalty");
        Sundance::setOption("solver", solverFile, "name of XML file for solver");

        Sundance::init( &argc , &argv );
        int np = MPIComm::world().getNProc();
        int npx = -1;
        int npy = -1;
        balanceXY(np, &npx, &npy);
        TEST_FOR_EXCEPT(npx < 1);
        TEST_FOR_EXCEPT(npy < 1);
        TEST_FOR_EXCEPT(npx * npy != np);

        VectorType<double> vecType = new EpetraVectorType();

        const int k = 1;
        const int splitBC = 1;

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

        BasisFamily L = new Lagrange( k );

        Expr u = new UnknownFunction( L , "u" );
        Expr v = new TestFunction( L , "v" );
        QuadratureFamily quad = new GaussianQuadrature( 2 * k );

        Expr h = new CellDiameterExpr();
        Expr alpha = C / h;
        Expr eqn = poissonEquationNitsche( splitBC, u , v , alpha , quad );
        Expr bc;


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

#ifdef HAVE_CONFIG_H
        ParameterXMLFileReader reader(searchForFile("SolverParameters/" + solverFile));
#else
        ParameterXMLFileReader reader(solverFile);
#endif
        ParameterList solverParams = reader.getParameters();
        LinearSolver<double> solver
            = LinearSolverBuilder::createSolver(solverParams);


        Expr soln = prob.solve( solver );


        FieldWriter w = new VTKWriter( "NitschePoisson2D" );
        w.addMesh( mesh );
        w.addField( "u" , new ExprFieldWrapper( soln ) );
        w.write();

        Expr x = new CoordExpr(0);
        Expr y = new CoordExpr(1);
        QuadratureFamily quad4 = new GaussianQuadrature(4);

        CellFilter interior = new MaximalCellFilter();
        const double pi = 4.0*atan(1.0);
        Expr exactSoln = sin(pi*x)*sin(pi*y);
        Expr err = exactSoln - soln;
        Expr errExpr = Integral(interior,
                                err*err,
                                quad4);
        FunctionalEvaluator errInt(mesh, errExpr);

        double errorSq = errInt.evaluate();
        cout << "error norm = " << sqrt(errorSq) << std::endl << std::endl;


        Sundance::passFailTest(sqrt(errorSq), 1.0e-4);
    }
    catch (std::exception &e)
    {
        Sundance::handleException(e);
    }
    Sundance::finalize();
    return Sundance::testStatus();
}

/** 
 * This example program sets up and solves the Laplace 
 * equation \f$-\nabla^2 u=0\f$. See the
 * document GettingStarted.pdf for more information.
 */
int main(int argc, char** argv)
{
  try
  {
    /* command-line options */
    std::string meshFile="plateWithHole3D-1";
    std::string solverFile = "aztec-ml.xml";
    Sundance::setOption("meshFile", meshFile, "mesh file");
    Sundance::setOption("solver", solverFile, 
      "name of XML file for solver");

    /* Initialize */
    Sundance::init(&argc, &argv);

    /* --- Specify vector representation to be used --- */
    VectorType<double> vecType = new EpetraVectorType();

    /* --- Read mesh --- */
    MeshType meshType = new BasicSimplicialMeshType();
    MeshSource meshSrc
      =  new ExodusMeshReader(meshFile, meshType);
    Mesh mesh = meshSrc.getMesh();

    /* --- Specification of geometric regions --- */

    /* Region "interior" consists of all maximal-dimension cells */
    CellFilter interior = new MaximalCellFilter();

    /* Identify boundary regions via labels in mesh */
    CellFilter edges = new DimensionalCellFilter(2);

    CellFilter south = edges.labeledSubset(1);
    CellFilter east = edges.labeledSubset(2);
    CellFilter north = edges.labeledSubset(3);
    CellFilter west = edges.labeledSubset(4);
    CellFilter hole = edges.labeledSubset(5);
    CellFilter down = edges.labeledSubset(6);
    CellFilter up = edges.labeledSubset(7);

    /* --- Symbolic equation definition --- */

    /* Test and unknown function */
    BasisFamily basis = new Lagrange(1);
    Expr u = new UnknownFunction(basis, "u");
    Expr v = new TestFunction(basis, "v");

    /* Gradient operator */
    Expr dx = new Derivative(0);
    Expr dy = new Derivative(1);
    Expr dz = new Derivative(2);
    Expr grad = List(dx, dy, dz);

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


    /** Write the weak form */
    Expr eqn 
      = Integral(interior, (grad*u)*(grad*v), quad1)
      + Integral(east, v, quad1);

    /* Write the essential boundary conditions */
    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(west, v*u/h, quad2);

    /* Set up linear problem */
    LinearProblem prob(mesh, eqn, bc, v, u, vecType);

    /* --- solve the problem --- */

    /* Create the solver as specified by parameters in 
     * an XML file */

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

    /* Solve! The solution is returned as an Expr containing a 
    * DiscreteFunction */
    Expr soln = prob.solve(solver);

    /* --- Postprocessing --- */

    /* Project the derivative onto the P1 basis */
    DiscreteSpace discSpace(mesh, List(basis, basis, basis), vecType);
    L2Projector proj(discSpace, grad*soln);
    Expr gradU = proj.project();

    /* Write the solution and its projected gradient to a VTK file */
    FieldWriter w = new VTKWriter("LaplaceDemo3D");
    w.addMesh(mesh);
    w.addField("soln", new ExprFieldWrapper(soln[0]));
    w.addField("du_dx", new ExprFieldWrapper(gradU[0]));
    w.addField("du_dy", new ExprFieldWrapper(gradU[1]));
    w.addField("du_dz", new ExprFieldWrapper(gradU[2]));
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  
  try
  {
    Sundance::init(&argc, &argv);
    int np = MPIComm::world().getNProc();
    TEUCHOS_TEST_FOR_EXCEPT(np > 1); // works only in serial for now

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

    /* Create a mesh. This will be our coarsest mesh */
    MeshType meshType = new BasicSimplicialMeshType();
    int nx = 3;

    MeshSource mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, 1,
      0.0, 1.0, nx, 1,
      meshType);
    Mesh coarse = mesher.getMesh();
    /* Label some random vertex. We'll check that the refinement preserves
     * vertex labels */
    coarse.setLabel(0, 3, 5);

    /* Refine the mesh. The URP object will store both meshes as well as
     * maps between cell indices at the two levels. */
    UniformRefinementPair urp(meshType, coarse);
    /* Run a consistency check on the refinement maps */
    int bad = urp.check();

    /* Get the fine mesh */
    const Mesh& fine = urp.fine();

    /* Create some expressions to be tested. The expression is linear
     * and so can be represented exactly with the P1 basis. Interpolation
     * to the fine mesh will also be exact. */
    Expr x = new CoordExpr(0);
    Expr y = new CoordExpr(1);
    Expr f = x + sqrt(2.0)*y;

    BasisFamily P1 = new Lagrange(1);
    
    /* Project a function onto the coarse space */
    DiscreteSpace coarseP1(coarse, P1, vecType);
    L2Projector projCoarseP1(coarseP1, f);
    Expr cP1 = projCoarseP1.project();
    Vector<double> cxP1 = DiscreteFunction::discFunc(cP1)->getVector();
    RCP<DOFMapBase> cDofMap = coarseP1.map();

    /* Project the same function onto the fine space */
    DiscreteSpace fineP1(fine, P1, vecType);
    L2Projector projFineP1(fineP1, f);
    Expr fP1 = projFineP1.project();
    Vector<double> fxP1 = DiscreteFunction::discFunc(fP1)->getVector();
    RCP<DOFMapBase> fDofMap = fineP1.map();

    /* Now we're going to interpolate the coarse function onto the fine
     * space. If the code is correct this will be equal to the function 
     * projected onto the fine mesh. */
    Vector<double> refinedX = fxP1.space().createMember();
    
    /* map the coarse space to the fine space */
    for (int f=0; f<fine.numCells(0); f++)
    {
      /* get the DOF for the node */
      int fDof = getNodeDof(fDofMap, f);
      /* If this fine vertex is on an edge of the coarse mesh, interpolate
       * from the facets of the edge. If the fine vertex is coincident
       * with a coarse vertex, use the value on that vertex */
      if (urp.newVertIsOnEdge()[f]) /* new vertex is on a coarse edge */
      {
        /* get the parent edge's LID in the coarse mesh */
        int cEdge = urp.newVertToOldLIDMap()[f];
        int ori; // dummy orientation, not needed here
        int cNode1 = coarse.facetLID(1, cEdge, 0, 0, ori);
        int cNode2 = coarse.facetLID(1, cEdge, 0, 1, ori);
        int cDof1 = getNodeDof(cDofMap, cNode1);
        int cDof2 = getNodeDof(cDofMap, cNode2);
        refinedX[fDof] = 0.5*(cxP1[cDof1]+cxP1[cDof2]);
      }
      else /* new vertex coincides with an old vertex */
      {
        int cNode = urp.newVertToOldLIDMap()[f];
        int cDof = getNodeDof(cDofMap, cNode);
        refinedX[fDof] = cxP1[cDof];
      }
    }
    
    /* Compare the R-mapped vector to the vector computed natively 
     * on the fine space */
    double err = (fxP1 - refinedX).norm2();
    cout << "error |fine - R*coarse| = " << err << endl;

    
    double finalTol = 0.5;
    Sundance::passFailTest(bad, finalTol);
  }
	catch(std::exception& e)
  {
    Sundance::handleException(e);
//.........这里部分代码省略.........

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

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

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

      //  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)
{
  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();
      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(); 
}

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

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