C++ GEOObjects::getPointVec方法代码示例

bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects &geo_objects, double eps)
{
    if (mesh.getDimension() != 2)
    {
        ERR ("Mesh to geometry conversion is only working for 2D meshes.");
        return false;
    }

    // nodes to points conversion
    std::string mesh_name(mesh.getName());
    {
        auto points = std::make_unique<std::vector<GeoLib::Point*>>();
        points->reserve(mesh.getNumberOfNodes());

        for (auto node_ptr : mesh.getNodes())
            points->push_back(new GeoLib::Point(*node_ptr, node_ptr->getID()));

        geo_objects.addPointVec(std::move(points), mesh_name, nullptr, eps);
    }
    const std::vector<std::size_t> id_map (geo_objects.getPointVecObj(mesh_name)->getIDMap());

    // elements to surface triangles conversion
    std::string const mat_name ("MaterialIDs");
    auto bounds (MeshInformation::getValueBounds<int>(mesh, mat_name));
    const unsigned nMatGroups(bounds.second-bounds.first+1);
    auto sfcs = std::make_unique<std::vector<GeoLib::Surface*>>();
    sfcs->reserve(nMatGroups);
    auto const& points = *geo_objects.getPointVec(mesh_name);
    for (unsigned i=0; i<nMatGroups; ++i)
        sfcs->push_back(new GeoLib::Surface(points));

    const std::vector<MeshLib::Element*> &elements = mesh.getElements();
    const std::size_t nElems (mesh.getNumberOfElements());

    MeshLib::PropertyVector<int> const*const materialIds =
        mesh.getProperties().existsPropertyVector<int>("MaterialIDs")
            ? mesh.getProperties().getPropertyVector<int>("MaterialIDs")
            : nullptr;

    for (unsigned i=0; i<nElems; ++i)
    {
        auto surfaceId = !materialIds ? 0 : ((*materialIds)[i] - bounds.first);
        MeshLib::Element* e (elements[i]);
        if (e->getGeomType() == MeshElemType::TRIANGLE)
            (*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
        if (e->getGeomType() == MeshElemType::QUAD)
        {
            (*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
            (*sfcs)[surfaceId]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(2)], id_map[e->getNodeIndex(3)]);
        }
        // all other element types are ignored (i.e. lines)
    }

    std::for_each(sfcs->begin(), sfcs->end(), [](GeoLib::Surface* sfc){ if (sfc->getNumberOfTriangles()==0) delete sfc; sfc = nullptr;});
    auto sfcs_end = std::remove(sfcs->begin(), sfcs->end(), nullptr);
    sfcs->erase(sfcs_end, sfcs->end());

    geo_objects.addSurfaceVec(std::move(sfcs), mesh_name);
    return true;
}

	void createSurfaces()
	{
		auto const points = geo_objects.getPointVec(geo_name);
		const std::vector<std::size_t> pnt_id_map(
			geo_objects.getPointVecObj(geo_name)->getIDMap());

		auto sfcs = std::unique_ptr<std::vector<GeoLib::Surface*>>(
			new std::vector<GeoLib::Surface*>(2));
		std::map<std::string, std::size_t>* sfc_names =
			new std::map<std::string, std::size_t>;
		(*sfcs)[0] = new GeoLib::Surface(*points);
		(*sfcs)[0]->addTriangle(pnt_id_map[0], pnt_id_map[3], pnt_id_map[1]);
		(*sfcs)[0]->addTriangle(pnt_id_map[1], pnt_id_map[3], pnt_id_map[4]);
		(*sfcs)[0]->addTriangle(pnt_id_map[1], pnt_id_map[4], pnt_id_map[2]);
		(*sfcs)[0]->addTriangle(pnt_id_map[2], pnt_id_map[4], pnt_id_map[5]);
		(*sfc_names)["FirstSurface"] = 0;
		(*sfcs)[1] = new GeoLib::Surface(*points);
		(*sfcs)[1]->addTriangle(pnt_id_map[3], pnt_id_map[6], pnt_id_map[8]);
		(*sfcs)[1]->addTriangle(pnt_id_map[3], pnt_id_map[8], pnt_id_map[5]);
		(*sfc_names)["SecondSurface"] = 1;
		geo_objects.addSurfaceVec(std::move(sfcs), geo_name, sfc_names);
	}

TEST(GeoLib, SurfaceIsPointInSurface)
{
    std::vector<std::function<double(double, double)>> surface_functions;
    surface_functions.push_back(constant);
    surface_functions.push_back(coscos);

    for (auto f : surface_functions) {
        std::random_device rd;

        std::string name("Surface");
        // generate ll and ur in random way
        std::mt19937 random_engine_mt19937(rd());
        std::normal_distribution<> normal_dist_ll(-10, 2);
        std::normal_distribution<> normal_dist_ur(10, 2);
        MathLib::Point3d ll(std::array<double,3>({{
                normal_dist_ll(random_engine_mt19937),
                normal_dist_ll(random_engine_mt19937),
                0.0
            }
        }));
        MathLib::Point3d ur(std::array<double,3>({{
                normal_dist_ur(random_engine_mt19937),
                normal_dist_ur(random_engine_mt19937),
                0.0
            }
        }));
        for (std::size_t k(0); k<3; ++k)
            if (ll[k] > ur[k])
                std::swap(ll[k], ur[k]);

        // random discretization of the domain
        std::default_random_engine re(rd());
        std::uniform_int_distribution<std::size_t> uniform_dist(2, 25);
        std::array<std::size_t,2> n_steps = {{uniform_dist(re),uniform_dist(re)}};

        std::unique_ptr<MeshLib::Mesh> sfc_mesh(
            MeshLib::MeshGenerator::createSurfaceMesh(
                name, ll, ur, n_steps, f
            )
        );

        // random rotation angles
        std::normal_distribution<> normal_dist_angles(
            0, boost::math::double_constants::two_pi);
        std::array<double,3> euler_angles = {{
                normal_dist_angles(random_engine_mt19937),
                normal_dist_angles(random_engine_mt19937),
                normal_dist_angles(random_engine_mt19937)
            }
        };

        MathLib::DenseMatrix<double, std::size_t> rot_mat(getRotMat(
                    euler_angles[0], euler_angles[1], euler_angles[2]));

        std::vector<MeshLib::Node*> const& nodes(sfc_mesh->getNodes());
        GeoLib::rotatePoints<MeshLib::Node>(rot_mat, nodes);

        MathLib::Vector3 const normal(0,0,1.0);
        MathLib::Vector3 const surface_normal(rot_mat * normal);
        double const eps(1e-6);
        MathLib::Vector3 const displacement(eps * surface_normal);

        GeoLib::GEOObjects geometries;
        MeshLib::convertMeshToGeo(*sfc_mesh, geometries);

        std::vector<GeoLib::Surface*> const& sfcs(*geometries.getSurfaceVec(name));
        GeoLib::Surface const*const sfc(sfcs.front());
        std::vector<GeoLib::Point*> const& pnts(*geometries.getPointVec(name));
        // test triangle edge point of the surface triangles
        for (auto const p : pnts) {
            EXPECT_TRUE(sfc->isPntInSfc(*p));
            MathLib::Point3d q(*p);
            for (std::size_t k(0); k<3; ++k)
                q[k] += displacement[k];
            EXPECT_FALSE(sfc->isPntInSfc(q));
        }
        // test edge middle points of the triangles
        for (std::size_t k(0); k<sfc->getNTriangles(); ++k) {
            MathLib::Point3d p, q, r;
            std::tie(p,q,r) = getEdgeMiddlePoints(*(*sfc)[k]);
            EXPECT_TRUE(sfc->isPntInSfc(p));
            EXPECT_TRUE(sfc->isPntInSfc(q));
            EXPECT_TRUE(sfc->isPntInSfc(r));
        }
    }
}