C++ TetrahedralMesh::GetBoundaryElementIteratorBegin方法代码示例

    void TestAnyNonZeroNeumannConditionsAndApplyNeumannToMeshBoundary()
    {
        // Load a 2D square mesh with 1 central non-boundary node
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_4_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        BoundaryConditionsContainer<2,2,1> bcc;
        BoundaryConditionsContainer<2,2,2> bcc_2unknowns;

        TS_ASSERT_EQUALS(bcc.AnyNonZeroNeumannConditions(), false);

        bcc.DefineZeroNeumannOnMeshBoundary(&mesh);

        TetrahedralMesh<2,2>::BoundaryElementIterator iter;
        iter = mesh.GetBoundaryElementIteratorBegin();
        while (iter != mesh.GetBoundaryElementIteratorEnd())
        {
            TS_ASSERT(bcc.HasNeumannBoundaryCondition(*iter));
            double value = bcc.GetNeumannBCValue(*iter, (*iter)->GetNode(0)->GetPoint());
            TS_ASSERT_DELTA(value, 0.0, 1e-8);

            iter++;
        }
        TS_ASSERT_EQUALS(bcc.AnyNonZeroNeumannConditions(), false);

        iter = mesh.GetBoundaryElementIteratorBegin();

        ConstBoundaryCondition<2>* p_boundary_condition2 = new ConstBoundaryCondition<2>(-1);

        bcc_2unknowns.AddNeumannBoundaryCondition(*iter, p_boundary_condition2);
        TS_ASSERT_EQUALS(bcc_2unknowns.AnyNonZeroNeumannConditions(), true);
    }

    /* Define a particular test. Note the {{{throw(Exception)}}} at the end of the
     * declaration. This causes {{{Exception}}} messages to be printed out if an
     * {{{Exception}}} is thrown, rather than just getting the message "terminate
     * called after throwing an instance of 'Exception' " */
    void TestSolvingNonlinearEllipticPde() throw(Exception)
    {
        /* As usual, first create a mesh. */
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        /* Next, instantiate the PDE to be solved. */
        MyNonlinearPde pde;

        /*
         * Then we have to define the boundary conditions. First, the Dirichlet boundary
         * condition, u=0 on x=0, using the boundary node iterator.
         */
        BoundaryConditionsContainer<2,2,1> bcc;
        ConstBoundaryCondition<2>* p_zero_bc = new ConstBoundaryCondition<2>(0.0);
        for (TetrahedralMesh<2,2>::BoundaryNodeIterator node_iter = mesh.GetBoundaryNodeIteratorBegin();
             node_iter != mesh.GetBoundaryNodeIteratorEnd();
             node_iter++)
        {
            if (fabs((*node_iter)->GetPoint()[1]) < 1e-12)
            {
                bcc.AddDirichletBoundaryCondition(*node_iter, p_zero_bc);
            }
        }

        /* And then the Neumman conditions. Neumann boundary condition are defined on
         * surface elements, and for this problem, the Neumman boundary value depends
         * on the position in space, so we make use of the {{{FunctionalBoundaryCondition}}}
         * object, which contains a pointer to a function, and just returns the value
         * of that function for the required point when the {{{GetValue}}} method is called.
         */
        FunctionalBoundaryCondition<2>* p_functional_bc = new FunctionalBoundaryCondition<2>(&MyNeummanFunction);
        /* Loop over surface elements. */
        for (TetrahedralMesh<2,2>::BoundaryElementIterator elt_iter = mesh.GetBoundaryElementIteratorBegin();
             elt_iter != mesh.GetBoundaryElementIteratorEnd();
             elt_iter++)
        {
            /* Get the y value of any node (here, the zero-th). */
            double y = (*elt_iter)->GetNodeLocation(0,1);
            /* If y=1... */
            if (fabs(y-1.0) < 1e-12)
            {
                /* ... then associate the functional boundary condition, (Dgradu).n = y,
                 *  with the surface element... */
                bcc.AddNeumannBoundaryCondition(*elt_iter, p_functional_bc);
            }
            else
            {
                /* ...else associate the zero boundary condition (i.e. zero flux) with this
                 * element. */
                bcc.AddNeumannBoundaryCondition(*elt_iter, p_zero_bc);
            }
        }
        /* Note that in the above loop, the zero Neumman boundary condition was applied
         * to all surface elements for which y!=1, which included the Dirichlet surface
         * y=0. This is OK, as Dirichlet boundary conditions are applied to the finite
         * element matrix after Neumman boundary conditions, where the appropriate rows
         * in the matrix are overwritten.
         *
         * This is the solver for solving nonlinear problems, which, as usual,
         * takes in the mesh, the PDE, and the boundary conditions. */
        SimpleNonlinearEllipticSolver<2,2> solver(&mesh, &pde, &bcc);

        /* The solver also needs to be given an initial guess, which will be
         * a PETSc vector. We can make use of a helper method to create it.
         */
        Vec initial_guess = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 0.25);

        /* '''Optional:''' To use Chaste's Newton solver to solve nonlinear vector equations that are
         * assembled, rather than the default PETSc nonlinear solvers, we can
         * do the following: */
        SimpleNewtonNonlinearSolver newton_solver;
        solver.SetNonlinearSolver(&newton_solver);
        /* '''Optional:''' We can also manually set tolerances, and whether to print statistics, with
         * this nonlinear vector equation solver */
        newton_solver.SetTolerance(1e-10);
        newton_solver.SetWriteStats();

        /* Now call {{{Solve}}}, passing in the initial guess */
        Vec answer = solver.Solve(initial_guess);

        /* Note that we could have got the solver to not use an analytical Jacobian
         * and use a numerically-calculated Jacobian instead, by passing in false as a second
         * parameter:
         */
        //Vec answer = solver.Solve(initial_guess, false);

        /* Once solved, we can check the obtained solution against the analytical
         * solution. */
        ReplicatableVector answer_repl(answer);
        for (unsigned i=0; i<answer_repl.GetSize(); i++)
        {
            double y = mesh.GetNode(i)->GetPoint()[1];
            double exact_u = sqrt(y*(4-y));
            TS_ASSERT_DELTA(answer_repl[i], exact_u, 0.15);
//.........这里部分代码省略.........

    void TestHeatEquationWithSourceWithCoupledOdeSystemIn1dWithZeroNeumann()
    {
        // Create mesh of domain [0,1]
        TrianglesMeshReader<1,1> mesh_reader("mesh/test/data/1D_0_to_1_100_elements");
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Create PDE system object
        HeatEquationWithSourceForCoupledOdeSystem<1> pde;

        // Define zero Neumann boundary conditions
        BoundaryConditionsContainer<1,1,1> bcc;
        ConstBoundaryCondition<1>* p_boundary_condition = new ConstBoundaryCondition<1>(0.0);
        TetrahedralMesh<1,1>::BoundaryElementIterator iter = mesh.GetBoundaryElementIteratorBegin();
        bcc.AddNeumannBoundaryCondition(*iter, p_boundary_condition);
        iter = mesh.GetBoundaryElementIteratorEnd();
        iter--;
        bcc.AddNeumannBoundaryCondition(*iter, p_boundary_condition);

        // Create the correct number of ODE systems
        double a = 5.0;
        std::vector<AbstractOdeSystemForCoupledPdeSystem*> ode_systems;
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            ode_systems.push_back(new OdeSystemForCoupledHeatEquationWithSource(a));
        }

        // Create PDE system solver
        LinearParabolicPdeSystemWithCoupledOdeSystemSolver<1,1,1> solver(&mesh, &pde, &bcc, ode_systems);

        // Test setting end time and timestep
        TS_ASSERT_THROWS_THIS(solver.SetTimes(1.0, 0.0), "Start time has to be less than end time");
        TS_ASSERT_THROWS_THIS(solver.SetTimeStep(0.0), "Time step has to be greater than zero");

        // Set end time and timestep
        double t_end = 0.1;
        solver.SetTimes(0, t_end);
        solver.SetTimeStep(0.001);

        // Set initial condition u(x,0) = 1 + cos(pi*x)
        std::vector<double> init_cond(mesh.GetNumNodes());
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];
            init_cond[i] = 1 + cos(M_PI*x);
        }
        Vec initial_condition = PetscTools::CreateVec(init_cond);
        solver.SetInitialCondition(initial_condition);

        // Solve PDE system and store result
        Vec result = solver.Solve();
        ReplicatableVector result_repl(result);

        /*
         * Test that solution is given by
         *
         * u(x,t) = 1 + (1 - exp(-a*t))/a + exp(-pi*pi*t)*cos(pi*x),
         * v(x,t) = exp(-a*t),
         *
         * with t = t_end.
         */
        for (unsigned i=0; i<result_repl.GetSize(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];

            double u = 1 + (1 - exp(-a*t_end))/a + exp(-M_PI*M_PI*t_end)*cos(M_PI*x);
            TS_ASSERT_DELTA(result_repl[i], u, 0.1);

            double u_from_v = solver.GetOdeSystemAtNode(i)->rGetPdeSolution()[0];
            TS_ASSERT_DELTA(result_repl[i], u_from_v, 0.1);

            double v = exp(-a*t_end);
            TS_ASSERT_DELTA(ode_systems[i]->rGetStateVariables()[0], v, 0.1);
        }

        // Test the method GetOdeSystemAtNode()
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            TS_ASSERT(solver.GetOdeSystemAtNode(i) != NULL);
            TS_ASSERT_DELTA(static_cast<OdeSystemForCoupledHeatEquationWithSource*>(solver.GetOdeSystemAtNode(i))->GetA(), 5.0, 1e-6);
        }

        // Tidy up
        PetscTools::Destroy(initial_condition);
        PetscTools::Destroy(result);
    }

    void TestSolvingEllipticPde() throw(Exception)
    {
        /* First we declare a mesh reader which reads mesh data files of the 'Triangle'
         * format. The path given is relative to the main Chaste directory. As we are in 2d,
         * the reader will look for three datafiles, [name].nodes, [name].ele and [name].edge.
         * Note that the first template argument here is the spatial dimension of the
         * elements in the mesh ({{{ELEMENT_DIM}}}), and the second is the dimension of the nodes,
         * i.e. the dimension of the space the mesh lives in ({{{SPACE_DIM}}}). Usually
         * {{{ELEMENT_DIM}}} and {{{SPACE_DIM}}} will be equal. */
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
        /* Now declare a tetrahedral mesh with the same dimensions... */
        TetrahedralMesh<2,2> mesh;
        /* ... and construct the mesh using the mesh reader. */
        mesh.ConstructFromMeshReader(mesh_reader);

        /* Next we instantiate an instance of our PDE we wish to solve. */
        MyPde pde;

        /* A set of boundary conditions are stored in a {{{BoundaryConditionsContainer}}}. The
         * three template arguments are ELEMENT_DIM, SPACE_DIM and PROBLEM_DIM, the latter being
         * the number of unknowns we are solving for. We have one unknown (ie u is a scalar, not
         * a vector), so in this case {{{PROBLEM_DIM}}}=1. */
        BoundaryConditionsContainer<2,2,1> bcc;

        /* Defining the boundary conditions is the only particularly fiddly part of solving PDEs,
         * unless they are very simple, such as u=0 on the boundary, which could be done
         * as follows: */
        //bcc.DefineZeroDirichletOnMeshBoundary(&mesh);

        /* We want to specify u=0 on x=0 and y=0.  To do this, we first create the boundary condition
         * object saying what the value of the condition is at any particular point in space.  Here
         * we use the class `ConstBoundaryCondition`, a subclass of `AbstractBoundaryCondition` that
         * yields the same constant value (0.0 here) everywhere it is used.
         *
         * Note that the object is allocated with `new`.  The `BoundaryConditionsContainer` object deals
         * with deleting its associated boundary condition objects.  Note too that we could allocate a
         * separate condition object for each boundary node, but using the same object where possible is
         * more memory efficient.
         */
        ConstBoundaryCondition<2>* p_zero_boundary_condition = new ConstBoundaryCondition<2>(0.0);
        /* We then get a boundary node iterator from the mesh... */
        TetrahedralMesh<2,2>::BoundaryNodeIterator iter = mesh.GetBoundaryNodeIteratorBegin();
        /* ...and loop over the boundary nodes, getting the x and y values. */
        while (iter < mesh.GetBoundaryNodeIteratorEnd())
        {
            double x = (*iter)->GetPoint()[0];
            double y = (*iter)->GetPoint()[1];
            /* If x=0 or y=0... */
            if ((x==0) || (y==0))
            {
                /* ...associate the zero boundary condition created above with this boundary node
                 * ({{{*iter}}} being a pointer to a {{{Node<2>}}}).
                 */
                bcc.AddDirichletBoundaryCondition(*iter, p_zero_boundary_condition);
            }
            iter++;
        }

        /* Now we create Neumann boundary conditions for the ''surface elements'' on x=1 and y=1. Note that
         * Dirichlet boundary conditions are defined on nodes, whereas Neumann boundary conditions are
         * defined on surface elements. Note also that the natural boundary condition statement for this
         * PDE is (D grad u).n = g(x) (where n is the outward-facing surface normal), and g(x) is a prescribed
         * function, ''not'' something like du/dn=g(x). Hence the boundary condition we are specifying is
         * (D grad u).n = 0.
         *
         * '''Important note for 1D:''' This means that if we were solving 2u,,xx,,=f(x) in 1D, and
         * wanted to specify du/dx=1 on the LHS boundary, the Neumann boundary value we have to specify is
         * -2, as n=-1 (outward facing normal) so (D gradu).n = -2 when du/dx=1.
         *
         * To define Neumann bcs, we reuse the zero boundary condition object defined above, but apply it
         * at surface elements.  We loop over these using another iterator provided by the mesh class.
         */
        TetrahedralMesh<2,2>::BoundaryElementIterator surf_iter
            = mesh.GetBoundaryElementIteratorBegin();
        while (surf_iter != mesh.GetBoundaryElementIteratorEnd())
        {
            /* Get the x and y values of any node (here, the 0th) in the element. */
            unsigned node_index = (*surf_iter)->GetNodeGlobalIndex(0);
            double x = mesh.GetNode(node_index)->GetPoint()[0];
            double y = mesh.GetNode(node_index)->GetPoint()[1];

            /* If x=1 or y=1... */
            if ( (fabs(x-1.0) < 1e-6) || (fabs(y-1.0) < 1e-6) )
            {
                /* ...associate the boundary condition with the surface element. */
                bcc.AddNeumannBoundaryCondition(*surf_iter, p_zero_boundary_condition);
            }

            /* Finally increment the iterator. */
            surf_iter++;
        }

        /* Next we define the solver of the PDE.
         * To solve an {{{AbstractLinearEllipticPde}}} (which is the type of PDE {{{MyPde}}} is),
         * we use a {{{SimpleLinearEllipticSolver}}}. The solver, again templated over
         * {{{ELEMENT_DIM}}} and {{{SPACE_DIM}}}, needs to be given (pointers to) the mesh,
         * pde and boundary conditions.
         */
        SimpleLinearEllipticSolver<2,2> solver(&mesh, &pde, &bcc);

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

    // In this test we have no cardiac tissue, so that the equations are just sigma * phi_e''=0
    // throughout the domain (with a Neumann boundary condition on x=1 and a dirichlet boundary
    // condition (ie grounding) on x=0), so the exact solution can be calculated and compared
    // against.
    void Test1dProblemOnlyBathGroundedOneSide() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(0.5);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainBathOnlyBath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_bath");

        c_vector<double,1> centre;
        centre(0) = 0.5;
        BathCellFactory<1> cell_factory(-1e6, centre);

        TrianglesMeshReader<1,1> reader("mesh/test/data/1D_0_to_1_10_elements");
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructFromMeshReader(reader);

        for(unsigned i=0; i<mesh.GetNumElements(); i++)
        {
            mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidBathId());
        }

        // create boundary conditions container
        double boundary_val = 1.0;
        boost::shared_ptr<BoundaryConditionsContainer<1,1,2> > p_bcc(new BoundaryConditionsContainer<1,1,2>);
        ConstBoundaryCondition<1>* p_bc_stim = new ConstBoundaryCondition<1>(boundary_val);
        ConstBoundaryCondition<1>* p_zero_stim = new ConstBoundaryCondition<1>(0.0);

        // loop over boundary elements and set (sigma\gradphi).n = 1.0 on RHS edge
        for(TetrahedralMesh<1,1>::BoundaryElementIterator iter
              = mesh.GetBoundaryElementIteratorBegin();
           iter != mesh.GetBoundaryElementIteratorEnd();
           iter++)
        {
            if (((*iter)->GetNodeLocation(0))[0]==1.0)
            {
                /// \todo: I think you need to provide a boundary condition for unknown#1 if you are gonig to provide one for unknown#2?
                p_bcc->AddNeumannBoundaryCondition(*iter, p_zero_stim, 0);
                p_bcc->AddNeumannBoundaryCondition(*iter, p_bc_stim,   1);
            }
        }

        BidomainWithBathProblem<1> bidomain_problem( &cell_factory );

        bidomain_problem.SetBoundaryConditionsContainer(p_bcc);
        bidomain_problem.SetMesh(&mesh);
        bidomain_problem.Initialise();

        // fix phi=0 on LHS edge
        std::vector<unsigned> fixed_nodes;
        fixed_nodes.push_back(0);
        bidomain_problem.SetFixedExtracellularPotentialNodes(fixed_nodes);

        bidomain_problem.Solve();

        Vec sol = bidomain_problem.GetSolution();
        ReplicatableVector sol_repl(sol);

        // test phi = x*boundary_val/sigma (solution of phi''=0, phi(0)=0, sigma*phi'(1)=boundary_val
        for(unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double bath_cond = HeartConfig::Instance()->GetBathConductivity();
            double x = mesh.GetNode(i)->rGetLocation()[0];
            TS_ASSERT_DELTA(sol_repl[2*i],   0.0,   1e-12);               // V
            TS_ASSERT_DELTA(sol_repl[2*i+1], x*boundary_val/bath_cond, 1e-4);   // phi_e
        }
    }

    // test 2D problem - takes a long time to run.
    // solution is incorrect to specified tolerance.
    void xTestSimpleLinearParabolicSolver2DNeumannWithSmallTimeStepAndFineMesh()
    {
        // Create mesh from mesh reader
        FemlabMeshReader<2,2> mesh_reader("mesh/test/data/",
                                          "femlab_fine_square_nodes.dat",
                                          "femlab_fine_square_elements.dat",
                                          "femlab_fine_square_edges.dat");

        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Instantiate PDE object
        HeatEquation<2> pde;

        // Boundary conditions - zero dirichlet on boundary;
        BoundaryConditionsContainer<2,2,1> bcc;
        TetrahedralMesh<2,2>::BoundaryNodeIterator iter = mesh.GetBoundaryNodeIteratorBegin();

        while (iter != mesh.GetBoundaryNodeIteratorEnd())
        {
            double x = (*iter)->GetPoint()[0];
            double y = (*iter)->GetPoint()[1];

            ConstBoundaryCondition<2>* p_dirichlet_boundary_condition =
                new ConstBoundaryCondition<2>(x);

            if (fabs(y) < 0.01)
            {
                bcc.AddDirichletBoundaryCondition(*iter, p_dirichlet_boundary_condition);
            }

            if (fabs(y - 1.0) < 0.01)
            {
                bcc.AddDirichletBoundaryCondition(*iter, p_dirichlet_boundary_condition);
            }

            if (fabs(x) < 0.01)
            {
                bcc.AddDirichletBoundaryCondition(*iter, p_dirichlet_boundary_condition);
            }

            iter++;
        }

        TetrahedralMesh<2,2>::BoundaryElementIterator surf_iter = mesh.GetBoundaryElementIteratorBegin();
        ConstBoundaryCondition<2>* p_neumann_boundary_condition =
            new ConstBoundaryCondition<2>(1.0);

        while (surf_iter != mesh.GetBoundaryElementIteratorEnd())
        {
            int node = (*surf_iter)->GetNodeGlobalIndex(0);
            double x = mesh.GetNode(node)->GetPoint()[0];

            if (fabs(x - 1.0) < 0.01)
            {
                bcc.AddNeumannBoundaryCondition(*surf_iter, p_neumann_boundary_condition);
            }

            surf_iter++;
        }

        // Solver
        SimpleLinearParabolicSolver<2,2> solver(&mesh,&pde,&bcc);

        // Initial condition u(0,x,y) = sin(0.5*M_PI*x)*sin(M_PI*y)+x
        std::vector<double> init_cond(mesh.GetNumNodes());
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];
            double y = mesh.GetNode(i)->GetPoint()[1];
            init_cond[i] = sin(0.5*M_PI*x)*sin(M_PI*y)+x;
        }
        Vec initial_condition = PetscTools::CreateVec(init_cond);

        double t_end = 0.1;
        solver.SetTimes(0, t_end);
        solver.SetTimeStep(0.001);

        solver.SetInitialCondition(initial_condition);

        Vec result = solver.Solve();
        ReplicatableVector result_repl(result);

        // Check solution is u = e^{-5/4*M_PI*M_PI*t} sin(0.5*M_PI*x)*sin(M_PI*y)+x, t=0.1
        for (unsigned i=0; i<result_repl.GetSize(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];
            double y = mesh.GetNode(i)->GetPoint()[1];
            double u = exp((-5/4)*M_PI*M_PI*t_end) * sin(0.5*M_PI*x) * sin(M_PI*y) + x;
            TS_ASSERT_DELTA(result_repl[i], u, 0.001);
        }

        PetscTools::Destroy(result);
        PetscTools::Destroy(initial_condition);
    }

    /**
     * Simple Parabolic PDE u' = del squared u
     *
     * With u = x on 5 boundaries of the unit cube, and
     * u_n = 1 on the x face of the cube.
     *
     * Subject to the initial condition
     * u(0,x,y,z)=sin( PI x)sin( PI y)sin( PI z) + x
     */
    void TestSimpleLinearParabolicSolver3DNeumannOnCoarseMesh()
    {
        // Create mesh from mesh reader
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");

        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Instantiate PDE object
        HeatEquation<3> pde;

        // Boundary conditions
        BoundaryConditionsContainer<3,3,1> bcc;
        TetrahedralMesh<3,3>::BoundaryNodeIterator iter = mesh.GetBoundaryNodeIteratorBegin();

        while (iter != mesh.GetBoundaryNodeIteratorEnd())
        {
            double x = (*iter)->GetPoint()[0];
            double y = (*iter)->GetPoint()[1];
            double z = (*iter)->GetPoint()[2];


            if ((fabs(y) < 0.01) || (fabs(y - 1.0) < 0.01) ||
                (fabs(x) < 0.01) ||
                (fabs(z) < 0.01) || (fabs(z - 1.0) < 0.01) )
            {
                ConstBoundaryCondition<3>* p_dirichlet_boundary_condition =
                    new ConstBoundaryCondition<3>(x);
                bcc.AddDirichletBoundaryCondition(*iter, p_dirichlet_boundary_condition);
            }

            iter++;
        }

        TetrahedralMesh<3,3>::BoundaryElementIterator surf_iter = mesh.GetBoundaryElementIteratorBegin();
        ConstBoundaryCondition<3>* p_neumann_boundary_condition =
            new ConstBoundaryCondition<3>(1.0);

        while (surf_iter != mesh.GetBoundaryElementIteratorEnd())
        {
            int node = (*surf_iter)->GetNodeGlobalIndex(0);
            double x = mesh.GetNode(node)->GetPoint()[0];

            if (fabs(x - 1.0) < 0.01)
            {
                bcc.AddNeumannBoundaryCondition(*surf_iter, p_neumann_boundary_condition);
            }

            surf_iter++;
        }

        // Solver
        SimpleLinearParabolicSolver<3,3> solver(&mesh,&pde,&bcc);

        // Initial condition u(0,x,y) = sin(0.5*PI*x)*sin(PI*y)+x
        Vec initial_condition = PetscTools::CreateVec(mesh.GetNumNodes());

        double* p_initial_condition;
        VecGetArray(initial_condition, &p_initial_condition);

        int lo, hi;
        VecGetOwnershipRange(initial_condition, &lo, &hi);
        for (int global_index = lo; global_index < hi; global_index++)
        {
            int local_index = global_index - lo;
            double x = mesh.GetNode(global_index)->GetPoint()[0];
            double y = mesh.GetNode(global_index)->GetPoint()[1];
            double z = mesh.GetNode(global_index)->GetPoint()[2];

            p_initial_condition[local_index] = sin(0.5*M_PI*x)*sin(M_PI*y)*sin(M_PI*z)+x;
        }
        VecRestoreArray(initial_condition, &p_initial_condition);

        solver.SetTimes(0, 0.1);
        solver.SetTimeStep(0.01);

        solver.SetInitialCondition(initial_condition);
        Vec result = solver.Solve();

        // Check result
        double* p_result;
        VecGetArray(result, &p_result);

        // Solution should be u = e^{-5/2*PI*PI*t} sin(0.5*PI*x)*sin(PI*y)*sin(PI*z)+x, t=0.1
        for (int global_index = lo; global_index < hi; global_index++)
        {
            int local_index = global_index - lo;
            double x = mesh.GetNode(global_index)->GetPoint()[0];
            double y = mesh.GetNode(global_index)->GetPoint()[1];
            double z = mesh.GetNode(global_index)->GetPoint()[2];

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