C++ DistributedTetrahedralMesh::ConstructFromMeshReader方法代码示例

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

        c_matrix<double, 3, 3> dummy;
        double jacobian_det = 0.0;

        double scaled_volume = 0.0;
        for (AbstractTetrahedralMesh<3, 3>::ElementIterator iter = mesh.GetElementIteratorBegin();
                iter != mesh.GetElementIteratorEnd();
                ++iter)
        {
            iter->CalculateJacobian(dummy, jacobian_det);
            scaled_volume += jacobian_det;
        }
        TS_ASSERT_DELTA(scaled_volume, 1.0*6, 1e-6);

        mesh.RotateX(M_PI);
        //mesh.Translate(100.0, 0.0, 0.0);

        double scaled_volume_after = 0.0;
        for (AbstractTetrahedralMesh<3, 3>::ElementIterator iter = mesh.GetElementIteratorBegin();
                iter != mesh.GetElementIteratorEnd();
                ++iter)
        {
            iter->CalculateJacobian(dummy, jacobian_det);
            scaled_volume_after += jacobian_det;
        }
        TS_ASSERT_DELTA(scaled_volume_after, 1.0*6, 1e-6);

    }

    void TestCheckForBathElementsNoDeadlock() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("bidomain_bath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory;
        BidomainWithBathProblem<1> bidomain_problem( &bidomain_cell_factory );

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

        try
        {
            mesh.GetElement(0)->SetAttribute(HeartRegionCode::GetValidBathId());
        }
        catch(Exception&)
        {
            // I don't own element 0
        }

        bidomain_problem.SetMesh(&mesh);

        // Fails because no bath
        TS_ASSERT_THROWS_NOTHING(bidomain_problem.Initialise());

        // Prevent an EventHandling exception in later tests
        HeartEventHandler::EndEvent(HeartEventHandler::EVERYTHING);
    }

    void TestSetLogInfo() throw (Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");
        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";

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

        StreeterFibreGenerator<3> fibre_generator(mesh);
        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        fibre_generator.SetApexToBase(0);

        OutputFileHandler handler("shorter_streeter_loginfo");

        fibre_generator.WriteData(handler, "box_heart.ortho");

        FileFinder node_regions_file = handler.FindFile("node_regions.data");
        FileFinder wall_thickness_file = handler.FindFile("wall_thickness.data");
        FileFinder averaged_thickness_file = handler.FindFile("averaged_thickness.data");

        TS_ASSERT_EQUALS(node_regions_file.IsFile(), false);
        TS_ASSERT_EQUALS(wall_thickness_file.IsFile(), false);
        TS_ASSERT_EQUALS(averaged_thickness_file.IsFile(), false);

        fibre_generator.SetLogInfo(true);

        fibre_generator.WriteData(handler, "box_heart.ortho");

        TS_ASSERT_EQUALS(node_regions_file.IsFile(), true);
        TS_ASSERT_EQUALS(wall_thickness_file.IsFile(), true);
        TS_ASSERT_EQUALS(averaged_thickness_file.IsFile(), true);
    }

    void TestConstructStreeterOnRightWedge() throw(Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/human_wedge_mesh/HumanWedgeMesh");
        std::string epi_face_file = "heart/test/data/human_wedge_mesh/epi.tri";
        std::string endo_face_file = "heart/test/data/human_wedge_mesh/endo.tri";

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

        StreeterFibreGenerator<3> fibre_generator(mesh);

        //Assume we are in the left ventricle
        fibre_generator.SetSurfaceFiles(epi_face_file, endo_face_file, "", true);

        fibre_generator.SetApexToBase(0);
        fibre_generator.SetWriteFileAsBinary();

        OutputFileHandler handler("human_wedge_mesh", false);

        fibre_generator.WriteData(handler, "HumanWedgeMeshRight.ortho");

        FileFinder fibre_file1 = handler.FindFile("HumanWedgeMeshRight.ortho");
        FileFinder fibre_file2("heart/test/data/human_wedge_mesh/HumanWedgeMeshRight.ortho", RelativeTo::ChasteSourceRoot);

        CompareGeneratedWithReferenceFile(fibre_file1, ORTHO, fibre_file2, ORTHO);
    }

    void TestSimpleOrthotropic() throw (Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");
        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";

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

        StreeterFibreGenerator<3> fibre_generator(mesh);
        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        fibre_generator.SetApexToBase(0);

        OutputFileHandler handler("shorter_streeter", false);
        fibre_generator.WriteData(handler, "box_heart.ortho");

        FileFinder fibre_file_ascii = handler.FindFile("box_heart.ortho");
        FileFinder fibre_file_reference("heart/test/data/box_shaped_heart/box_heart.ortho", RelativeTo::ChasteSourceRoot);

        CompareGeneratedWithReferenceFile(fibre_file_ascii, ORTHO, fibre_file_reference, ORTHO);

        fibre_generator.SetWriteFileAsBinary();
        fibre_generator.WriteData(handler, "box_heart_binary.ortho");

        FileFinder fibre_file_binary = handler.FindFile("box_heart_binary.ortho");

        CompareGeneratedWithReferenceFile(fibre_file_binary, ORTHO, fibre_file_reference, ORTHO);
    }

    void TestExceptions()
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");

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

        StreeterFibreGenerator<3> fibre_generator(mesh);

        // No surfaces defined
        OutputFileHandler handler("streeter", false);
        TS_ASSERT_THROWS_THIS(fibre_generator.WriteData(handler, "file.fibres"),
                "Files defining the heart surfaces not set");

        // Wrong surface filename
        TS_ASSERT_THROWS_THIS(fibre_generator.SetSurfaceFiles("wrong_name", "wrong_name", "wrong_name", false),
                "Wrong surface definition file name wrong_name");

        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";


        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        OutputFileHandler shorter_handler("shorter_streeter", false);
        TS_ASSERT_THROWS_THIS(fibre_generator.WriteData(shorter_handler, "vector_not_set.ortho"),
            "Apex to base vector has not been set");
        TS_ASSERT_THROWS_THIS(fibre_generator.SetApexToBase(999),
            "Apex to base coordinate axis was out of range");

        c_vector<double, 3> axis;
        axis[0] = 0.0;
        axis[1] = 0.0;
        axis[2] = 0.0;
        TS_ASSERT_THROWS_THIS(fibre_generator.SetApexToBase(axis),
            "Apex to base vector should be non-zero");
        axis[1] = 42.0; //Will be normalised
        fibre_generator.SetApexToBase(axis);
    }

    /**
     * Tests archiving of the tissue object.
     * It creates one, changes the default values of some member variables and saves.
     * Then it tries to load from the archive and checks that the member variables are with the right values.
     */
    void TestSaveAndLoadExtendedBidomainTissue() throw (Exception)
    {
        HeartConfig::Instance()->Reset();
        // Archive settings
        FileFinder archive_dir("extended_tissue_archive", RelativeTo::ChasteTestOutput);
        std::string archive_file = "extended_bidomain_tissue.arch";

        bool cache_replication_saved = false;
        double saved_printing_timestep = 2.0;
        double default_printing_timestep = HeartConfig::Instance()->GetPrintingTimeStep();

        c_matrix<double, 3, 3> intra_tensor_before_archiving;
        c_matrix<double, 3, 3> intra_tensor_second_cell_before_archiving;
        c_matrix<double, 3, 3> extra_tensor_before_archiving;

        //creation and save
        {
            // This call is required to set the appropriate conductivity media and to make sure that HeartConfig
            // knows the mesh filename despite we use our own mesh reader.
            HeartConfig::Instance()->SetMeshFileName("mesh/test/data/cube_136_elements");

            TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
            DistributedTetrahedralMesh<3,3> mesh;
            mesh.ConstructFromMeshReader(mesh_reader);

            UnStimulatedCellFactory first_cell;
            StimulatedCellFactory second_cell;
            ExtracellularStimulusFactory extra_factory;
            first_cell.SetMesh(&mesh);
            second_cell.SetMesh(&mesh);
            extra_factory.SetMesh(&mesh);
            ExtendedBidomainTissue<3> extended_tissue( &first_cell, &second_cell , &extra_factory);

            //set a value different from default for the conductivities of the second cell
            extended_tissue.SetIntracellularConductivitiesSecondCell(Create_c_vector(25.0,26.0,27.0));
            //this is normally done by the problem class, but in this test we do it manually
            extended_tissue.CreateIntracellularConductivityTensorSecondCell();
            extended_tissue.SetCacheReplication(cache_replication_saved); // Not the default to check it is archived...

            //shuffle default values to check if they get archived properly
            extended_tissue.SetAmFirstCell(11.0);
            extended_tissue.SetAmSecondCell(22.0);
            extended_tissue.SetAmGap(33.0);
            extended_tissue.SetCmFirstCell(44.0);
            extended_tissue.SetCmSecondCell(55.0);
            extended_tissue.SetGGap(66.0);

            //again, away from default value to check for archiving
            extended_tissue.SetUserSuppliedExtracellularStimulus(true);

            //set some heterogeneities in Ggap
            std::vector<boost::shared_ptr<AbstractChasteRegion<3> > > heterogeneity_areas;
            std::vector<double> Ggap_values;
            ChastePoint<3> cornerA(-1, -1, 0);
            ChastePoint<3> cornerB(0.001, 0.001, 0.001);
            boost::shared_ptr<ChasteCuboid<3> > p_cuboid_1(new ChasteCuboid<3>(cornerA, cornerB));
            heterogeneity_areas.push_back(p_cuboid_1);
            //within the first area
            Ggap_values.push_back(143.0);
            extended_tissue.SetGgapHeterogeneities(heterogeneity_areas, Ggap_values);
            extended_tissue.CreateGGapConductivities();

            // Some checks to make sure HeartConfig is being saved and loaded by this too.
            HeartConfig::Instance()->SetPrintingTimeStep(saved_printing_timestep);
            TS_ASSERT_DELTA(HeartConfig::Instance()->GetPrintingTimeStep(), saved_printing_timestep, 1e-9);

            intra_tensor_before_archiving = extended_tissue.rGetIntracellularConductivityTensor(0);
            intra_tensor_second_cell_before_archiving = extended_tissue.rGetIntracellularConductivityTensorSecondCell(0);
            extra_tensor_before_archiving = extended_tissue.rGetExtracellularConductivityTensor(0);

            // Save
            ArchiveOpener<boost::archive::text_oarchive, std::ofstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_oarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacTissue<3>* const p_archive_bidomain_tissue = &extended_tissue;
            (*p_arch) << p_archive_bidomain_tissue;

            HeartConfig::Reset();
            TS_ASSERT_DELTA(HeartConfig::Instance()->GetPrintingTimeStep(), default_printing_timestep, 1e-9);
            TS_ASSERT_DIFFERS(saved_printing_timestep, default_printing_timestep);
        }
        //load
        {
            ArchiveOpener<boost::archive::text_iarchive, std::ifstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_iarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacTissue<3>* p_abstract_tissue;
            (*p_arch) >> p_abstract_tissue;
            assert(p_abstract_tissue!=NULL);

            //dynamic cast so we are able to test specific variables of ExtendedBidomainTissue
            ExtendedBidomainTissue<3>* p_extended_tissue = dynamic_cast<ExtendedBidomainTissue<3>*>(p_abstract_tissue);
            assert(p_extended_tissue != NULL);

            const c_matrix<double, 3, 3>& intra_tensor_after_archiving = p_extended_tissue->rGetIntracellularConductivityTensor(0);
//.........这里部分代码省略.........

    void TestWithBathAndElectrodes()
    {
        /* First, set the end time and output info. In this simulation
         * we'll explicitly read the mesh, alter it, then pass it
         * to the problem class, so we don't set the mesh file name.
         */
        HeartConfig::Instance()->SetSimulationDuration(3.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainTutorialWithBath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");

        /* Bath problems seem to require decreased ODE timesteps.
         */
        HeartConfig::Instance()->SetOdeTimeStep(0.001);  //ms

        /* Use the {{{PlaneStimulusCellFactory}}} to define a set
         * of Luo-Rudy cells. We pass the stimulus magnitude as 0.0
         * as we don't want any stimulated cells.
         */
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellMLBackwardEuler,2> cell_factory(0.0);

        /*
         * Now, we load up a rectangular mesh (in triangle/tetgen format), done as follows,
         * using {{{TrianglesMeshReader}}}.  Note that we use a distributed mesh, so the data
         * is shared among processes if run in parallel.
         */
        TrianglesMeshReader<2,2> reader("mesh/test/data/2D_0_to_1mm_400_elements");
        DistributedTetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(reader);

        /*
         * In most simulations there is one valid tissue identifier and one valid bath identifier
         * (for elements).
         * One of these can be assigned to an element with
         *  * {{{mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidTissueId());}}}
         *  * {{{mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidBathId());}}}
         *
         * If we want heterogeneous conductivities outside the heart (for example for torso and blood)
         * then we will need different identifiers:
         */
        std::set<unsigned> tissue_ids;
        static unsigned tissue_id=0;
        tissue_ids.insert(tissue_id);

        std::set<unsigned> bath_ids;
        static unsigned bath_id1=1;
        bath_ids.insert(bath_id1);
        static unsigned bath_id2=2;
        bath_ids.insert(bath_id2);

        HeartConfig::Instance()->SetTissueAndBathIdentifiers(tissue_ids, bath_ids);

        /* In bath problems, each element has an attribute which must be set
         * to 0 (cardiac tissue) or 1 (bath). This can be done by having an
         * extra column in the element file (see the file formats documentation,
         * or for example
         * mesh/test/data/1D_0_to_1_10_elements_with_two_attributes.ele,
         * and note that the header in this file has 1 at the end to indicate that
         * the file defines an attribute for each element). We have read in a mesh
         * without this type of information set up, so we set it up manually,
         * by looping over elements and setting those more than 2mm from the centre
         * as bath elements (by default, the others are cardiac elements).
         */
        for (AbstractTetrahedralMesh<2,2>::ElementIterator iter = mesh.GetElementIteratorBegin();
             iter != mesh.GetElementIteratorEnd();
             ++iter)
        {
            double x = iter->CalculateCentroid()[0];
            double y = iter->CalculateCentroid()[1];
            if (sqrt((x-0.05)*(x-0.05) + (y-0.05)*(y-0.05)) > 0.02)
            {
                if (y<0.05)
                {
                    //Outside circle on the bottom
                    iter->SetAttribute(bath_id1);
                }
                else
                {
                    //Outside circle on the top
                    iter->SetAttribute(bath_id2);
                }
            }
            else
            {
                //IDs default to 0, but we want to be safe
                iter->SetAttribute(tissue_id);
            }
        }

        /* HOW_TO_TAG Cardiac/Problem definition
         * Tell Chaste that a mesh has been modified
         *
         * Since we have modified the mesh by setting element attributes, we need to inform Chaste of this fact.
         * If we do not, problems will arise when [wiki:UserTutorials/CardiacCheckpointingAndRestarting checkpointing],
         * since the code that saves the simulation state will assume that it can just reuse the original mesh files,
         * and thus won't save the new element attributes.
         *
         * (Some mesh modifications, that use methods on the mesh class directly, will automatically record that
         * the mesh has been modified.  Since we're just modifying elements, this information isn't propagated at
         * present.)
         */
//.........这里部分代码省略.........