C++ CellFilter类代码示例

CellFilter CellFilter::operator-(const CellFilter& other) const 
{
  if (other.isNull())
  {
    return *this;
  }
  else if (isKnownDisjointWith(other) || other.isKnownDisjointWith(*this))
  {
    return *this;
  }
  else if (isKnownSubsetOf(other))
  {
    CellFilter rtn;
    return rtn;
  }
  else if (*this == other)
  {
    CellFilter rtn;
    return rtn;
  }
  else
  {
    CellFilter rtn 
      = new BinaryCellFilter(*this, other, BinaryCellFilter::Difference);
    rtn.registerDisjoint(other);
    this->registerSubset(rtn);
    return rtn;
  }
}

HomogeneousDOFMap::HomogeneousDOFMap(const Mesh& mesh, 
  const BasisFamily& basis,
  int numFuncs, 
  int setupVerb)
  : DOFMapBase(mesh, setupVerb), 
    dim_(mesh.spatialDim()),
    dofs_(mesh.spatialDim()+1),
    maximalDofs_(),
    haveMaximalDofs_(false),
    localNodePtrs_(mesh.spatialDim()+1),
    nNodesPerCell_(mesh.spatialDim()+1),
    totalNNodesPerCell_(mesh.spatialDim()+1, 0),
    numFacets_(mesh.spatialDim()+1),
    originalFacetOrientation_(2),
    basisIsContinuous_(false)
{
  verbosity() = DOFMapBase::classVerbosity();
  
  CellFilter maximalCells = new MaximalCellFilter();
  cellSets().append(maximalCells.getCells(mesh));
  cellDimOnCellSets().append(mesh.spatialDim());

  allocate(mesh, basis, numFuncs);
  initMap();
}

void CellFilter::registerDisjoint(const CellFilter& sub) const
{
  SubsetManager::registerDisjoint(*this, sub);
  
  for (Set<CellFilter>::const_iterator 
         i=sub.knownDisjoints().begin(); i!=sub.knownDisjoints().end(); i++)
  {
    SubsetManager::registerDisjoint(*this, *i);
  }
}

/* weak form of poisson with Nitsche-type weak BC's */
Expr poissonEquationNitsche( bool splitBC,
                             Expr u ,
                             Expr v ,
                             Expr alpha ,
                             QuadratureFamily quad )
{
    CellFilter interior = new MaximalCellFilter();
    CellFilter boundary = new BoundaryCellFilter();
    CellFilter left = boundary.subset( new LeftPointTest() );
    CellFilter right = boundary.subset( new RightPointTest() );
    CellFilter top = boundary.subset( new TopPointTest() );
    CellFilter bottom = boundary.subset( new BottomPointTest() );

    CellFilter allBdry = left+right+top+bottom;

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

    Expr uvTerm;
    if (splitBC)
    {
        Out::os() << "BC expressions split over domains" << std::endl;
        uvTerm = Integral( left , alpha*u * v , quad )
                 + Integral( right , alpha*u * v , quad )
                 + Integral( top , alpha*u * v , quad )
                 + Integral( bottom , alpha*u * v , quad );
    }
    else
    {
        Out::os() << "BC expressions not split over domains" << std::endl;
        uvTerm = Integral( allBdry , alpha*u * v , quad );
    }


    const double pi = 4.0*atan(1.0);
    Expr force = 2.0*pi*pi*sin(pi*x)*sin(pi*y);
    return Integral( interior , (grad*v) * (grad*u) - force * v , quad )
           /* du/dn term */
           - Integral( left , -(dx*u)*v , quad )
           - Integral( top , (dy*u)*v , quad )
           - Integral( right , (dx*u)*v , quad )
           - Integral( bottom , -(dy*u)*v , quad )
           /* dv/dn term */
           - Integral( left , -(dx*v)*u , quad )
           - Integral( top , (dy*v)*u , quad )
           - Integral( right , (dx*v)*u , quad )
           - Integral( bottom , -(dy*v)*u , quad )
           /* u,v term  -- alpha = C / h */
           + uvTerm;
}

DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
  const RCP<FunctionSupportResolver>& fsr,
  const VectorType<double>& vecType,
  int setupVerb)
  : setupVerb_(setupVerb),
    map_(), 
    mesh_(mesh), 
    subdomains_(),
    basis_(basis),
    vecSpace_(), 
    vecType_(vecType),
    ghostImporter_()
  ,transformationBuilder_(new DiscreteSpaceTransfBuilder())
{
  bool partitionBCs = false;
  DOFMapBuilder builder(mesh, fsr, partitionBCs, setupVerb);

  map_ = builder.colMap()[0];
  Array<Set<CellFilter> > cf = builder.unkCellFilters()[0];

  for (int i=0; i<cf.size(); i++)
  {
    Array<Array<CellFilter> > dimCF(mesh.spatialDim()+1);
    for (Set<CellFilter>::const_iterator 
           iter=cf[i].begin(); iter != cf[i].end(); iter++)
    {
      CellFilter f = *iter;
      int dim = f.dimension(mesh);
      dimCF[dim].append(f);
    }
    for (int d=mesh.spatialDim(); d>=0; d--)
    {
      if (dimCF[d].size() == 0) continue;
      CellFilter f = dimCF[d][0];
      for (int j=1; j<dimCF[d].size(); j++)
      {
        f = f + dimCF[d][j];
      }
      subdomains_.append(f);
      break;
    }
  }
  RCP<Array<int> > dummyBCIndices;
  
  // set up the transformation
  transformationBuilder_ = rcp(new DiscreteSpaceTransfBuilder( mesh , basis , map_ ));

  initVectorSpace(dummyBCIndices, partitionBCs);
  initImporter();
}

bool CellFilter::isKnownDisjointWith(const CellFilter& other) const
{
  if (other.knownDisjoints().contains(*this)) return true;
  if (this->knownDisjoints().contains(other)) return true;

  return false;
}

CellSet connectedNodeSet(const CellFilter& f, const Mesh& mesh)
{
  CellSet cells = f.getCells(mesh);
  int dim = cells.dimension();
  if (dim==0) return cells;


  Array<int> cellLID;

  for (CellIterator i=cells.begin(); i!=cells.end(); i++)
  {
    cellLID.append(*i);
  }

  Array<int> nodes;
  Array<int> fo;

  mesh.getFacetLIDs(dim, cellLID, 0, nodes, fo);

  Set<int> nodeSet;

  for (int i=0; i<nodes.size(); i++)
  {
    nodeSet.put(nodes[i]);
  }
  
  return CellSet(mesh, 0, PointCell, nodeSet);
}

bool CellFilter::isSubsetOf(const CellFilter& other,
  const Mesh& mesh) const
{
  if (isKnownSubsetOf(other)) 
  {
    return true;
  }
  else
  {
    CellSet myCells = getCells(mesh);
    CellSet yourCells = other.getCells(mesh);
    CellSet inter = myCells.setIntersection(yourCells);
    if (inter.begin() == inter.end()) return false;
    CellSet diff = myCells.setDifference(inter);
    return (diff.begin() == diff.end());
  }
}

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

int main(int argc, char** argv)
{
  try
  {
    /*
     * Initialization code
     */
    std::string meshFile="plateWithHole2D-1";
    std::string solverFile = "nox-aztec.xml";
    Sundance::setOption("meshFile", meshFile, "mesh file");
    Sundance::setOption("solver", solverFile, 
      "name of XML file for solver");

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

    // This next line is just a hack to deal with some 
    // transitional code in the
    // element integration logic. 
    Sundance::ElementIntegral::alwaysUseCofacets() = false;

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

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

    /* 
     * Specification of cell filters
     */
    CellFilter interior = new MaximalCellFilter();
    CellFilter edges = new DimensionalCellFilter(1);

    CellFilter south = edges.labeledSubset(1);
    CellFilter east = edges.labeledSubset(2);
    CellFilter north = edges.labeledSubset(3);
    CellFilter west = edges.labeledSubset(4);

    /* 
     * <Header level="subsubsection" name="symb_setup">
     * Setup of symbolic problem description
     * </Header>
     * 
     * Create unknown and test functions discretized on the space
     * first-order Lagrange polynomials. 
     */
    BasisFamily basis = new Lagrange(2);
    Expr u = new UnknownFunction(basis, "u");
    Expr v = new TestFunction(basis, "v");

    /* 
     * Create differential operators and coordinate functions. Directions
     * are indexed starting from zero. The \verb+List()+ function can 
     * collect expressions into a vector. 
     */
    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);

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

    /* 
     * Create the weak form and the BCs
     */
    Expr source=exp(u);
    Expr eqn 
      = Integral(interior, (grad*u)*(grad*v)+v*source, quad4);

    Expr h = new CellDiameterExpr();
    Expr bc = EssentialBC(west+east, v*(u-1.0)/h, quad2);

    /* 
     * <Header level="subsubsection" name="lin_prob">
     * Creation of initial guess
     * </Header>
     *
     * So far the setup has been almost identical to that for the linear
     * problem, the only difference being the nonlinear term in the
     * equation set. 
     */
    DiscreteSpace discSpace(mesh, basis, vecType);
    L2Projector proj(discSpace, 1.0);
    Expr u0 = proj.project();

    /* 
     * <Header level="subsubsection" name="lin_prob">
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
  try
		{
      int depth = 0;
      bool useCCode = false;
      Sundance::ElementIntegral::alwaysUseCofacets() = true;
      Sundance::clp().setOption("depth", &depth, "expression depth");
      Sundance::clp().setOption("C", "symb", &useCCode, "Code type (C or symbolic)");
      Sundance::init(&argc, &argv);

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

      /* Read the mesh */
      MeshType meshType = new BasicSimplicialMeshType();

      MeshSource mesher 
        = new ExodusMeshReader("cube-0.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();
      CellFilter faces = new DimensionalCellFilter(2);
      CellFilter side1 = faces.labeledSubset(1);
      CellFilter side2 = faces.labeledSubset(2);
      CellFilter side3 = faces.labeledSubset(3);
      CellFilter side4 = faces.labeledSubset(4);
      CellFilter side5 = faces.labeledSubset(5);
      CellFilter side6 = faces.labeledSubset(6);

      
      /* Create unknown and test functions, discretized using second-order
       * Lagrange interpolants */
      Expr u = new UnknownFunction(new Lagrange(1), "u");
      Expr v = new TestFunction(new Lagrange(1), "v");

      /* Create differential operator and coordinate functions */
      Expr dx = new Derivative(0);
      Expr dy = new Derivative(1);
      Expr dz = new Derivative(2);
      Expr grad = List(dx, dy, dz);
      Expr x = new CoordExpr(0);
      Expr y = new CoordExpr(1);
      Expr z = new CoordExpr(2);

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

      /* Define the weak form */
      //Expr eqn = Integral(interior, (grad*v)*(grad*u) + v, quad);
      
      Expr coeff = 1.0;
#ifdef FOR_TIMING
      if (useCCode)
      {
        coeff = Poly(depth, x);
      }
      else
      {
        for (int i=0; i<depth; i++)
        {
          Expr t = 1.0;
          for (int j=0; j<depth; j++) t = t*x;
          coeff = coeff + 2.0*t - t - t;
        }
      }
#endif
      Expr eqn = Integral(interior, coeff*(grad*v)*(grad*u) /*+ 2.0*v*/, quad2);

      /* Define the Dirichlet BC */
      Expr exactSoln = x;//(x + 1.0)*x - 1.0/4.0;
      Expr h = new CellDiameterExpr();

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

      Expr bc = EssentialBC(side4, v*(u-exactSoln), quad4)
        + EssentialBC(side6, v*(u-exactSoln), quad4, watchBC);

      /* We can now set up the linear problem! */
      LinearProblem prob(mesh, eqn, bc, v, u, vecType);

#ifdef HAVE_CONFIG_H
      ParameterXMLFileReader reader(searchForFile("SolverParameters/aztec-ml.xml"));
#else
      ParameterXMLFileReader reader("aztec-ml.xml");
#endif
      ParameterList solverParams = reader.getParameters();
      std::cerr << "params = " << solverParams << std::endl;


      LinearSolver<double> solver 
        = LinearSolverBuilder::createSolver(solverParams);
//.........这里部分代码省略.........

SubmaximalNodalDOFMap
::SubmaximalNodalDOFMap(const Mesh& mesh, 
  const CellFilter& cf,
  int nFuncs,
  int setupVerb)
  : DOFMapBase(mesh, setupVerb),
    dim_(0),
    nTotalFuncs_(nFuncs),
    domain_(cf),
    domains_(tuple(cf)),
    nodeLIDs_(),
    nodeDOFs_(),
    lidToPtrMap_(),
    mapStructure_()
{
  Tabs tab0(0);
  SUNDANCE_MSG1(setupVerb, tab0 << "in SubmaximalNodalDOFMap ctor");
  Tabs tab1;
  SUNDANCE_MSG2(setupVerb, tab1 << "domain " << domain_);
  SUNDANCE_MSG2(setupVerb, tab1 << "N funcs " << nFuncs);

  const MPIComm& comm = mesh.comm();
  int rank = comm.getRank();
  int nProc = comm.getNProc();
  
  dim_ = cf.dimension(mesh);  
  TEUCHOS_TEST_FOR_EXCEPT(dim_ != 0);

  CellSet nodes = cf.getCells(mesh);
  int nc = nodes.numCells();
  nodeLIDs_.reserve(nc);
  nodeDOFs_.reserve(nc);

  Array<Array<int> > remoteNodes(nProc);
  
  int nextDOF = 0;
  int k=0; 
  for (CellIterator c=nodes.begin(); c!=nodes.end(); c++, k++)
  {
    int nodeLID = *c;
    lidToPtrMap_.put(nodeLID, k);
    nodeLIDs_.append(nodeLID);
    int remoteOwner = rank;
    if (isRemote(0, nodeLID, remoteOwner))
    {
      int GID = mesh.mapLIDToGID(0, nodeLID);
      remoteNodes[remoteOwner].append(GID);
      for (int f=0; f<nFuncs; f++) nodeDOFs_.append(-1);
    }
    else
    {
      for (int f=0; f<nFuncs; f++) nodeDOFs_.append(nextDOF++);
    }
  }

  /* Compute offsets for each processor */
  int localCount = nextDOF;
  computeOffsets(localCount);
  
  /* Resolve remote DOF numbers */
  shareRemoteDOFs(remoteNodes);

  BasisFamily basis = new Lagrange(1);
  mapStructure_ = rcp(new MapStructure(nTotalFuncs_, basis.ptr()));
}

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