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

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

    /*
     * Define a particular test.
     */
    void TestSchnackenbergSystemOnButterflyMesh() throw (Exception)
    {
        /* As usual, we first create a mesh. Here we are using a 2d mesh of a butterfly-shaped domain. */
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/butterfly");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        /* We scale the mesh to an appropriate size. */
        mesh.Scale(0.2, 0.2);

        /* Next, we instantiate the PDE system to be solved. We pass the parameter values into the
         * constructor.  (The order is D,,1,,  D,,2,,  k,,1,,  k,,-1,,  k,,2,,  k,,3,,) */
        SchnackenbergCoupledPdeSystem<2> pde(1e-4, 1e-2, 0.1, 0.2, 0.3, 0.1);

        /*
         * Then we have to define the boundary conditions. As we are in 2d, {{{SPACE_DIM}}}=2 and
         * {{{ELEMENT_DIM}}}=2. We also have two unknowns u and v,
         * so in this case {{{PROBLEM_DIM}}}=2. The value of each boundary condition is
         * given by the spatially uniform steady state solution of the Schnackenberg system,
         * given by u = (k,,1,, + k,,2,,)/k,,-1,,, v = k,,2,,k,,-1,,^2^/k,,3,,(k,,1,, + k,,2,,)^2^.
         */
        BoundaryConditionsContainer<2,2,2> bcc;
        ConstBoundaryCondition<2>* p_bc_for_u = new ConstBoundaryCondition<2>(2.0);
        ConstBoundaryCondition<2>* p_bc_for_v = new ConstBoundaryCondition<2>(0.75);
        for (TetrahedralMesh<2,2>::BoundaryNodeIterator node_iter = mesh.GetBoundaryNodeIteratorBegin();
             node_iter != mesh.GetBoundaryNodeIteratorEnd();
             ++node_iter)
        {
            bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_u, 0);
            bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_v, 1);
        }

        /* This is the solver for solving coupled systems of linear parabolic PDEs and ODEs,
         * which takes in the mesh, the PDE system, the boundary conditions and optionally
         * a vector of ODE systems (one for each node in the mesh). Since in this example
         * we are solving a system of coupled PDEs only, we do not supply this last argument. */
        LinearParabolicPdeSystemWithCoupledOdeSystemSolver<2,2,2> solver(&mesh, &pde, &bcc);

        /* Then we set the end time and time step and the output directory to which results will be written. */
        double t_end = 10;
        solver.SetTimes(0, t_end);
        solver.SetTimeStep(1e-1);
        solver.SetSamplingTimeStep(1);
        solver.SetOutputDirectory("TestSchnackenbergSystemOnButterflyMesh");

        /* We create a vector of initial conditions for u and v that are random perturbations
         * of the spatially uniform steady state and pass this to the solver. */
        std::vector<double> init_conds(2*mesh.GetNumNodes());
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            init_conds[2*i] = fabs(2.0 + RandomNumberGenerator::Instance()->ranf());
            init_conds[2*i + 1] = fabs(0.75 + RandomNumberGenerator::Instance()->ranf());
        }
        Vec initial_condition = PetscTools::CreateVec(init_conds);
        solver.SetInitialCondition(initial_condition);

        /* We now solve the PDE system and write results to VTK files, for
         * visualization using Paraview.  Results will be written to CHASTE_TEST_OUTPUT/TestSchnackenbergSystemOnButterflyMesh
         * as a results.pvd file and several results_[time].vtu files.
         * You should see something like [[Image(u.png, 350px)]] for u and [[Image(v.png, 350px)]] for v.
         */
        solver.SolveAndWriteResultsToFile();

        /*
         * All PETSc {{{Vec}}}s should be destroyed when they are no longer needed.
         */
        PetscTools::Destroy(initial_condition);
    }

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

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

    // 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];

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