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

    void TestScalingWithMethod()
    {
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        double mesh_volume = mesh.GetVolume();

        mesh.Scale(1.0);
        TS_ASSERT_DELTA(mesh_volume, mesh.GetVolume(), 1e-6);

        mesh.Scale(2.0, 3.0, 4.0);
        TS_ASSERT_DELTA(24.0*mesh_volume, mesh.GetVolume(), 1e-6);

        ChastePoint<3> corner_after = mesh.GetNode(6)->GetPoint();
        TS_ASSERT_DELTA(corner_after[0], 2.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[1], 3.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[2], 4.0, 1e-7);

        mesh.Scale(0.5, 1.0/3.0, 0.25);

        TS_ASSERT_DELTA(mesh_volume,mesh.GetVolume(),1e-6);

        corner_after = mesh.GetNode(6)->GetPoint();
        TS_ASSERT_DELTA(corner_after[0], 1.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[1], 1.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[2], 1.0, 1e-7);
    }

    /*
     * 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 Test2DSimulations() throw(Exception)
    {
        double conductivity_scale = 1;
        double h = 0.01; // cm
        double ode_time_step = 0.005; //ms
        double pde_time_step = 0.01; //ms
        unsigned num_stims = 1;

        TetrahedralMesh<2,2> mesh;
        unsigned num_elem_x = (unsigned)(0.5/h);  // num elements to make 5mm
        unsigned num_elem_y = (unsigned)(0.5/h);  // num elements to make 5mm
        //unsigned num_elem_z = (unsigned)(0.15/h);// Num elements to make 0.3cm
        double pacing_cycle_length = 350;
        double stim_mag = -500000;
        double stim_dur = 3;
        double area = 0.005;

        mesh.ConstructRectangularMesh(num_elem_x, num_elem_y);
        mesh.Scale(h,h); // Get mesh into units of cm.

        std::string archive_dir_base("LongPostprocessing_archives/archive");
        std::string archive_dir_current;

        // Setup
        HeartConfig::Instance()->SetSimulationDuration(pacing_cycle_length); //ms
        HeartConfig::Instance()->SetOutputDirectory("LongPostprocessing");
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");

        // These lines make postprocessing fast or slow.
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(ode_time_step, pde_time_step, 10); // Leads to 10MB VTK file
        //HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(ode_time_step, pde_time_step, 0.01); // Leads to 1GB VTK file

        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.4*conductivity_scale*1.171, 1.4*conductivity_scale*1.171));
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1400.0); // 1/cm
        HeartConfig::Instance()->SetCapacitance(1.0); // uF/cm^2
        HeartConfig::Instance()->SetVisualizeWithMeshalyzer();
#ifdef CHASTE_VTK
        HeartConfig::Instance()->SetVisualizeWithVtk();
#endif

        std::vector<std::pair<double,double> > apds_requested;
        apds_requested.push_back(std::pair<double, double>(90,-30)); //repolarisation percentage and threshold
        HeartConfig::Instance()->SetApdMaps(apds_requested);
//        std::vector<double> excitation_threshold;
//        excitation_threshold.push_back(-30.0);
//        HeartConfig::Instance()->SetUpstrokeTimeMaps(excitation_threshold);
//        HeartConfig::Instance()->SetMaxUpstrokeVelocityMaps(excitation_threshold);

        for (unsigned stim_counter=0; stim_counter < num_stims; stim_counter++ )
        {
            // Load problem
            MonodomainProblem<2> *p_monodomain_problem;
            if (stim_counter==0)
            {
                PointStimulusCellFactory<2> cell_factory(stim_mag, stim_dur, pacing_cycle_length, area);
                p_monodomain_problem = new MonodomainProblem<2>( &cell_factory );
                p_monodomain_problem->SetMesh(&mesh);
                p_monodomain_problem->Initialise();
            }
            else
            {
                p_monodomain_problem = CardiacSimulationArchiver<MonodomainProblem<2> >::Load(archive_dir_current);
            }

            HeartConfig::Instance()->SetSimulationDuration((double) (stim_counter+1)*pacing_cycle_length); //ms

            // set new directories to work from
            std::stringstream stringoutput;
            stringoutput << stim_counter;
            std::string stim_counter_string = stringoutput.str();

            archive_dir_current = archive_dir_base + "_" + stim_counter_string;
            OutputFileHandler archive_directory(archive_dir_current, true); // Clean a folder for new results
            HeartConfig::Instance()->SetOutputFilenamePrefix("results_" + stim_counter_string);

            // Solve problem (this does the postprocessing too when HeartConfig options are set).
            p_monodomain_problem->Solve();

            HeartEventHandler::Headings();
            HeartEventHandler::Report();

            // Save problem to archive
            CardiacSimulationArchiver<MonodomainProblem<2> >::Save(*p_monodomain_problem, archive_dir_current, false);
            std::cout << "Archived to " << archive_dir_current << "\n" << std::flush;

            // Copy the postprocessing results into the archive folders so they aren't wiped.
            std::vector<std::string> files;
            files.push_back("Apd_90_minus_30_Map");
//                files.push_back("MaxUpstrokeVelocityMap_-30");
//                files.push_back("UpstrokeTimeMap_-30");

            for (unsigned i=0; i<files.size(); i++)
            {
                FileFinder file_to_copy(HeartConfig::Instance()->GetOutputDirectory() + "/output/" + files[i] + ".dat", RelativeTo::ChasteTestOutput);
                TS_ASSERT(file_to_copy.IsFile());
                archive_directory.CopyFileTo(file_to_copy);
            }
        }// close for loop
    }//close void Test2dSimulations