C++ geolib::GEOObjects类代码示例

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

int XmlStnInterface::readFile(const QString &fileName)
{
	GEOLIB::GEOObjects* geoObjects = _project->getGEOObjects();
	QFile* file = new QFile(fileName);
	if (!file->open(QIODevice::ReadOnly | QIODevice::Text))
	{
		std::cout << "XmlStnInterface::readFile() - Can't open xml-file." << "\n";
		delete file;
		return 0;
	}
	if (!checkHash(fileName))
	{
		delete file;
		return 0;
	}

	QDomDocument doc("OGS-STN-DOM");
	doc.setContent(file);
	QDomElement docElement = doc.documentElement(); //root element, used for identifying file-type
	if (docElement.nodeName().compare("OpenGeoSysSTN"))
	{
		std::cout << "XmlStnInterface::readFile() - Unexpected XML root." << "\n";
		delete file;
		return 0;
	}

	QDomNodeList lists = docElement.childNodes();
	for (int i = 0; i < lists.count(); i++)
	{
		// read all the station lists
		QDomNodeList stationList = lists.at(i).childNodes();
		std::vector<GEOLIB::Point*>* stations = new std::vector<GEOLIB::Point*>;
		std::string stnName("[NN]");

		for (int j = 0; j < stationList.count(); j++)
		{
			const QDomNode station_node(stationList.at(j));
			const QString station_type(station_node.nodeName());
			if (station_type.compare("name") == 0)
				stnName = station_node.toElement().text().toStdString();
			else if (station_type.compare("stations") == 0)
				this->readStations(station_node, stations, fileName.toStdString());
			else if (station_type.compare("boreholes") == 0)
				this->readStations(station_node, stations, fileName.toStdString());
		}

		if (!stations->empty())
			geoObjects->addStationVec(stations, stnName);
		else
			delete stations;
	}

	delete file;

	return 1;
}

GMSHPrefsDialog::GMSHPrefsDialog(GeoLib::GEOObjects const& geoObjects, QDialog* parent)
    : QDialog(parent), _allGeo(new QStringListModel), _selGeo(new QStringListModel)
{
    setupUi(this);

    // default parameters
    this->param1->setText("2");
    this->param2->setText("0.3");
    this->param3->setText("0.05");
    this->param4->setText("0");

    // object will be deleted by Qt
    auto* max_number_of_points_in_quadtree_leaf_validator(
        new StrictIntValidator(1, 1000, this->param1));
    param1->setValidator (max_number_of_points_in_quadtree_leaf_validator);
    // object will be deleted by Qt
    auto* mesh_density_scaling_pnts_validator(
        new StrictDoubleValidator(0, 1, 5, this->param2));
    param2->setValidator (mesh_density_scaling_pnts_validator);
    // object will be deleted by Qt#
    auto* mesh_density_scaling_stations_validator(
        new StrictDoubleValidator(0, 1, 5, this->param3));
    param3->setValidator (mesh_density_scaling_stations_validator);

    std::vector<std::string> geoNames;
    geoObjects.getGeometryNames(geoNames);

    // get station names
    std::vector<std::string> geo_station_names;
    geoObjects.getStationVectorNames(geo_station_names);

    for (auto& geo_station_name : geo_station_names)
        geoNames.push_back(geo_station_name);

    std::size_t nGeoObjects(geoNames.size());

    QStringList list;
    for (unsigned i = 0; i < nGeoObjects; ++i)
        list.append(QString::fromStdString(geoNames[i]));

    if (list.empty())
    {
        this->selectGeoButton->setDisabled(true);
        this->deselectGeoButton->setDisabled(true);
        list.append("[No geometry available.]");
    }
    _allGeo->setStringList(list);
    this->allGeoView->setModel(_allGeo);
    this->selectedGeoView->setModel(_selGeo);
    this->radioAdaptive->toggle(); // default is adaptive meshing
    this->on_radioAdaptive_toggled(true);
}

// geometry_sets contains the geometric points the boundary conditions will be
// set on, geo_name is the name the geometry can be accessed with, out_fname is
// the base file name the gli and bc as well as the gml file will be written to.
void writeBCsAndGeometry(GeoLib::GEOObjects& geometry_sets,
                         std::string& geo_name, std::string const& out_fname,
                         std::string const& bc_type, bool write_gml)
{
	if (write_gml) {
		INFO("write points to \"%s.gml\".", geo_name.c_str());
		FileIO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gml");
	}
	FileIO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gli");

	bool liquid_flow(false);
	if (bc_type == "LIQUID_FLOW")
		liquid_flow = true;


	GeoLib::PointVec const* pnt_vec_objs(geometry_sets.getPointVecObj(geo_name));
	std::vector<GeoLib::Point*> const& pnts(*(pnt_vec_objs->getVector()));
	std::ofstream bc_out (out_fname+".bc");
	for (std::size_t k(0); k<pnts.size(); k++) {
		std::string const& pnt_name(pnt_vec_objs->getItemNameByID(k));
		if (!pnt_name.empty()) {
			if (liquid_flow)
				writeLiquidFlowPointBC(bc_out, pnt_name);
			else
				writeGroundwaterFlowPointBC(bc_out, pnt_name, (*pnts[k])[2]);
		}
	}
	bc_out << "#STOP\n";
	bc_out.close();
}

TEST_F(OGSIOVer4InterfaceTest, InvalidTIN_ZeroAreaTri)
{
    std::string tin_fname(BaseLib::BuildInfo::tests_tmp_path+"Surface.tin");
    std::ofstream tin_out (tin_fname);
    tin_out << "0 0.0 0.0 0.0 1.0 0.0.0 0.0 0.0 0.0\n";
    tin_out.close();

    // read geometry
    GeoLib::GEOObjects geometries;
    std::vector<std::string> errors;
    std::string geometry_name("TestGeometry");
    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
                                  errors, "dummy_for_gmsh_path");

    std::vector<GeoLib::Surface*> const*
        sfcs(geometries.getSurfaceVec(geometry_name));
    ASSERT_TRUE(sfcs == nullptr);

    std::remove(tin_fname.c_str());
}

void convertMeshNodesToGeometry(std::vector<MeshLib::Node*> const& nodes,
	std::vector<std::size_t> const& node_ids,
	std::string & geo_name,
	GeoLib::GEOObjects & geometry_sets)
{
	// copy data
	auto pnts = std::unique_ptr<std::vector<GeoLib::Point*>>(
	    new std::vector<GeoLib::Point*>);
	std::map<std::string, std::size_t>* pnt_names(
		new std::map<std::string, std::size_t>);
	std::size_t cnt(0);
	for (std::size_t id: node_ids) {
		pnts->push_back(new GeoLib::Point(*(nodes[id]), cnt));
		pnt_names->insert(std::pair<std::string, std::size_t>(
			geo_name+"-PNT-"+std::to_string(cnt), cnt));
		cnt++;
	}

	// create data structures for geometry
	geometry_sets.addPointVec(std::move(pnts), geo_name, pnt_names);
}

void createSetOfTestPointsAndAssociatedNames(GeoLib::GEOObjects & geo_objs, std::string &name, GeoLib::Point const& shift)
{
	std::vector<GeoLib::Point*> *pnts(new std::vector<GeoLib::Point*>);
	std::map<std::string, std::size_t>* pnt_name_map(new std::map< std::string, std::size_t>);

	const std::size_t pnts_per_edge(8);
	for (std::size_t k(0); k < pnts_per_edge; k++) {
		const std::size_t k_offset(k * pnts_per_edge * pnts_per_edge);
		for (std::size_t j(0); j < pnts_per_edge; j++) {
			const std::size_t offset(j * pnts_per_edge + k_offset);
			for (std::size_t i(0); i < pnts_per_edge; i++) {
				pnts->push_back(new GeoLib::Point(i+shift[0], j+shift[1], k+shift[2]));
				std::string pnt_name(
						name + "-" + BaseLib::number2str(i) + "-" + BaseLib::number2str(j) + "-"
								+ BaseLib::number2str(k));
				pnt_name_map->insert(std::pair< std::string, std::size_t>(pnt_name, i + offset));
			}
		}
	}

	geo_objs.addPointVec(pnts, name, pnt_name_map);
}

ProcessVariable::ProcessVariable(
    ConfigTree const& config,
    MeshLib::Mesh const& mesh,
    GeoLib::GEOObjects const& geometries)
    : _name(config.get<std::string>("name")),
      _mesh(mesh)
{
    DBUG("Constructing process variable %s", this->_name.c_str());

    // Initial condition
    {
        auto const& ic_config = config.find("initial_condition");
        if (ic_config == config.not_found())
            INFO("No initial condition found.");


        std::string const type =
            config.get<std::string>("initial_condition.type");
        if (type == "Uniform")
        {
            _initial_condition.reset(new UniformInitialCondition(ic_config->second));
        }
        else
        {
            ERR("Unknown type of the initial condition.");
        }
    }

    // Boundary conditions
    {
        auto const& bcs_config = config.find("boundary_conditions");
        if (bcs_config == config.not_found())
            INFO("No boundary conditions found.");

        for (auto const& bc_iterator : bcs_config->second)
        {
            ConfigTree const& bc_config = bc_iterator.second;

            // Find corresponding GeoObject
            std::string const geometrical_set_name =
                bc_config.get<std::string>("geometrical_set");
            std::string const geometry_name =
                bc_config.get<std::string>("geometry");

            GeoLib::GeoObject const* const geometry = geometries.getGeoObject(
                    geometrical_set_name, geometry_name);
            DBUG("Found geometry type \"%s\"",
                GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str());

            // Construct type dependent boundary condition
            std::string const type = bc_config.get<std::string>("type");

            if (type == "UniformDirichlet")
            {
                _dirichlet_bcs.emplace_back(
                    new UniformDirichletBoundaryCondition(
                        geometry, bc_config));
            }
            else if (type == "UniformNeumann")
            {
                _neumann_bc_configs.emplace_back(
                    new NeumannBcConfig(geometry, bc_config));
            }
            else
            {
                ERR("Unknown type \'%s\' of the boundary condition.",
                        type.c_str());
            }
        }

    }
}

int XmlStnInterface::write(std::ostream& stream)
{
	if (this->_exportName.empty())
	{
		std::cout << "Error in XmlStnInterface::write() - No station list specified..." << "\n";
		return 0;
	}

	GEOLIB::GEOObjects* geoObjects = _project->getGEOObjects();

	stream << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>\n"; // xml definition
	stream << "<?xml-stylesheet type=\"text/xsl\" href=\"OpenGeoSysSTN.xsl\"?>\n\n"; // stylefile definition

	QDomDocument doc("OGS-STN-DOM");
	QDomElement root = doc.createElement("OpenGeoSysSTN");
	root.setAttribute( "xmlns:ogs", "http://www.opengeosys.net" );
	root.setAttribute( "xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance" );
	root.setAttribute( "xsi:noNamespaceSchemaLocation", "http://141.65.34.25/OpenGeoSysSTN.xsd" );

	const std::vector<GEOLIB::Point*>* stations (geoObjects->getStationVec(_exportName));
	bool isBorehole =
	        (static_cast<GEOLIB::Station*>((*stations)[0])->type() ==
	         GEOLIB::Station::BOREHOLE) ? true : false;

	doc.appendChild(root);
	QDomElement stationListTag = doc.createElement("stationlist");
	root.appendChild(stationListTag);

	QDomElement listNameTag = doc.createElement("name");
	stationListTag.appendChild(listNameTag);
	QDomText stationListNameText = doc.createTextNode(QString::fromStdString(_exportName));
	listNameTag.appendChild(stationListNameText);
	QString listType = (isBorehole) ? "boreholes" : "stations";
	QDomElement stationsTag = doc.createElement(listType);
	stationListTag.appendChild(stationsTag);

	bool useStationValue(false);
	double sValue=static_cast<GEOLIB::Station*>((*stations)[0])->getStationValue();
	size_t nStations(stations->size());
	for (size_t i = 1; i < nStations; i++)
		if ((static_cast<GEOLIB::Station*>((*stations)[i])->getStationValue() - sValue) < std::numeric_limits<double>::min())
		{
			useStationValue = true;
			break;
		}

	for (size_t i = 0; i < nStations; i++)
	{
		QString stationType =  (isBorehole) ? "borehole" : "station";
		QDomElement stationTag = doc.createElement(stationType);
		stationTag.setAttribute( "id", QString::number(i) );
		stationTag.setAttribute( "x",  QString::number((*(*stations)[i])[0], 'f') );
		stationTag.setAttribute( "y",  QString::number((*(*stations)[i])[1], 'f') );
		stationTag.setAttribute( "z",  QString::number((*(*stations)[i])[2], 'f') );
		stationsTag.appendChild(stationTag);

		QDomElement stationNameTag = doc.createElement("name");
		stationTag.appendChild(stationNameTag);
		QDomText stationNameText =
		        doc.createTextNode(QString::fromStdString(static_cast<GEOLIB::Station*>((*stations)[i])->getName()));
		stationNameTag.appendChild(stationNameText);

		if (useStationValue)
		{
			QDomElement stationValueTag = doc.createElement("value");
			stationTag.appendChild(stationValueTag);
			QDomText stationValueText =
					doc.createTextNode(QString::number(static_cast<GEOLIB::Station*>((*stations)[i])->getStationValue()));
			stationValueTag.appendChild(stationValueText);
		}

		if (isBorehole)
			writeBoreholeData(doc, stationTag,
			                  static_cast<GEOLIB::StationBorehole*>((*stations)[i]));
	}

	std::string xml = doc.toString().toStdString();
	stream << xml;
	return 1;
}

int main (int argc, char* argv[])
{
    ApplicationsLib::LogogSetup logog_setup;

    TCLAP::CmdLine cmd("Maps geometric objects to the surface of a given mesh."
        "The documentation is available at https://docs.opengeosys.org/docs/tools/model-preparation/map-geometric-object-to-the-surface-of-a-mesh",
        ' ',
        "0.1");
    TCLAP::ValueArg<std::string> mesh_in("m", "mesh-file",
        "the name of the file containing the mesh", true,
        "", "file name");
    cmd.add(mesh_in);
    TCLAP::ValueArg<std::string> input_geometry_fname("i", "input-geometry",
        "the name of the file containing the input geometry", true,
        "", "file name");
    cmd.add(input_geometry_fname);
    TCLAP::ValueArg<bool> additional_insert_mapping("a", "additional-insert-mapping",
        "if true advanced mapping algorithm will be applied, i.e. a new "
        "geometry will be created and possibly new points will be inserted.", false,
        true, "boolean value");
    cmd.add(additional_insert_mapping);
    TCLAP::ValueArg<std::string> output_geometry_fname("o", "output-geometry",
        "the name of the file containing the input geometry", true,
        "", "file name");
    cmd.add(output_geometry_fname);
    cmd.parse(argc, argv);

    // *** read geometry
    GeoLib::GEOObjects geometries;
    {
        GeoLib::IO::BoostXmlGmlInterface xml_io(geometries);
        if (xml_io.readFile(input_geometry_fname.getValue())) {
            INFO("Read geometry from file \"%s\".",
                input_geometry_fname.getValue().c_str());
        } else {
            return EXIT_FAILURE;
        }
    }

    std::string geo_name;
    {
        std::vector<std::string> geo_names;
        geometries.getGeometryNames(geo_names);
        geo_name = geo_names[0];
    }

    MeshGeoToolsLib::GeoMapper geo_mapper(geometries, geo_name);

    // *** read mesh
    std::unique_ptr<MeshLib::Mesh> mesh(
        MeshLib::IO::readMeshFromFile(mesh_in.getValue()));

    if (additional_insert_mapping.getValue()) {
        geo_mapper.advancedMapOnMesh(*mesh);
    } else {
        geo_mapper.mapOnMesh(mesh.get());
    }

    {
        GeoLib::IO::BoostXmlGmlInterface xml_io(geometries);
        xml_io.setNameForExport(geo_name);
        xml_io.writeToFile(output_geometry_fname.getValue());
    }
    return EXIT_SUCCESS;
}

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

int main (int argc, char* argv[])
{
	ApplicationsLib::LogogSetup logog_setup;

	TCLAP::CmdLine cmd(
		"Creates boundary conditions for mesh nodes along polylines."
		"The documentation is available at https://docs.opengeosys.org/docs/tools/model-preparation/create-boundary-conditions-along-a-polyline",
		' ',
		"0.1");
	TCLAP::ValueArg<bool> gml_arg("", "gml",
		"if switched on write found nodes to file in gml format", false, 0, "bool");
	cmd.add(gml_arg);

	TCLAP::ValueArg<std::string> output_base_fname("o", "output-base-file-name",
		"the base name of the file the output (geometry (gli) and boundary"\
		"condition (bc)) will be written to", true,
		"", "file name");
	cmd.add(output_base_fname);

	TCLAP::ValueArg<std::string> bc_type("t", "type",
		"the process type the boundary condition will be written for "\
		"currently LIQUID_FLOW (primary variable PRESSURE1) and "\
		"GROUNDWATER_FLOW (primary variable HEAD, default) are supported", true,
		"",
		"process type as string (LIQUID_FLOW or GROUNDWATER_FLOW (default))");
	cmd.add(bc_type);

	TCLAP::ValueArg<double> search_length_arg("s", "search-length",
		"The size of the search length. The default value is "
		"std::numeric_limits<double>::epsilon()", false,
		std::numeric_limits<double>::epsilon(), "floating point number");
	cmd.add(search_length_arg);

	TCLAP::ValueArg<std::string> geometry_fname("i", "input-geometry",
		"the name of the file containing the input geometry", true,
		"", "file name");
	cmd.add(geometry_fname);

	TCLAP::ValueArg<std::string> mesh_arg("m", "mesh-file",
		"the name of the file containing the mesh", true,
		"", "file name");
	cmd.add(mesh_arg);

	cmd.parse(argc, argv);

	// *** read mesh
	INFO("Reading mesh \"%s\" ... ", mesh_arg.getValue().c_str());
	MeshLib::Mesh * subsurface_mesh(FileIO::readMeshFromFile(mesh_arg.getValue()));
	INFO("done.");
	INFO("Extracting top surface of mesh \"%s\" ... ",
		mesh_arg.getValue().c_str());
	const MathLib::Vector3 dir(0,0,-1);
	double const angle(90);
	std::unique_ptr<MeshLib::Mesh> surface_mesh(
	    MeshLib::MeshSurfaceExtraction::getMeshSurface(*subsurface_mesh, dir,
	                                                   angle));
	INFO("done.");
	delete subsurface_mesh;
	subsurface_mesh = nullptr;

	// *** read geometry
	GeoLib::GEOObjects geometries;
	FileIO::readGeometryFromFile(geometry_fname.getValue(), geometries);

	std::string geo_name;
	{
		std::vector<std::string> geo_names;
		geometries.getGeometryNames(geo_names);
		geo_name = geo_names[0];
	}

	// *** check if the data is usable
	// *** get vector of polylines
	std::vector<GeoLib::Polyline*> const* plys(geometries.getPolylineVec(geo_name));
	if (!plys) {
		ERR("Could not get vector of polylines out of geometry \"%s\".",
			geo_name.c_str());
		return -1;
	}

	MeshGeoToolsLib::SearchLength search_length_strategy;
	if (search_length_arg.isSet()) {
		search_length_strategy =
			MeshGeoToolsLib::SearchLength(search_length_arg.getValue());
	}

	GeoLib::GEOObjects geometry_sets;
	MeshGeoToolsLib::MeshNodeSearcher mesh_searcher(*surface_mesh,
		search_length_strategy);
	for(std::size_t k(0); k<plys->size(); k++) {
		std::vector<std::size_t> ids
			(mesh_searcher.getMeshNodeIDsAlongPolyline(*((*plys)[k])));
		if (ids.empty())
			continue;
		std::string geo_name("Polyline-"+std::to_string(k));
		convertMeshNodesToGeometry(surface_mesh->getNodes(), ids, geo_name,
			geometry_sets);
	}

	// merge all together
//.........这里部分代码省略.........

bool XmlGspInterface::write()
{
	GeoLib::GEOObjects* geoObjects = _project.getGEOObjects();
	QFileInfo fi(QString::fromStdString(_filename));
	std::string path((fi.absolutePath()).toStdString() + "/");

	_out << "<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>\n"; // xml definition
	_out << "<?xml-stylesheet type=\"text/xsl\" href=\"OpenGeoSysProject.xsl\"?>\n\n"; // stylefile definition

	QDomDocument doc("OGS-PROJECT-DOM");
	QDomElement root = doc.createElement("OpenGeoSysProject");
	root.setAttribute( "xmlns:ogs", "http://www.opengeosys.org" );
	root.setAttribute( "xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance" );
	root.setAttribute( "xsi:noNamespaceSchemaLocation",
	                   "http://www.opengeosys.org/images/xsd/OpenGeoSysProject.xsd" );

	doc.appendChild(root);

	// GML
	std::vector<std::string> geoNames;
	geoObjects->getGeometryNames(geoNames);
	for (std::vector<std::string>::const_iterator it(geoNames.begin()); it != geoNames.end(); ++it)
	{
		// write GLI file
		XmlGmlInterface gml(*geoObjects);
		std::string name(*it);
		gml.setNameForExport(name);
		if (gml.writeToFile(std::string(path + name + ".gml")))
		{
			// write entry in project file
			QDomElement geoTag = doc.createElement("geo");
			root.appendChild(geoTag);
			QDomElement fileNameTag = doc.createElement("file");
			geoTag.appendChild(fileNameTag);
			QDomText fileNameText =
			        doc.createTextNode(QString::fromStdString(name + ".gml"));
			fileNameTag.appendChild(fileNameText);
		}
	}

	// MSH
	const std::vector<MeshLib::Mesh*> msh_vec = _project.getMeshObjects();
	for (std::vector<MeshLib::Mesh*>::const_iterator it(msh_vec.begin()); it != msh_vec.end(); ++it)
	{
		// write mesh file
		Legacy::MeshIO meshIO;
		meshIO.setMesh(*it);
		std::string fileName(path + (*it)->getName());
		meshIO.writeToFile(fileName);

		// write entry in project file
		QDomElement mshTag = doc.createElement("msh");
		root.appendChild(mshTag);
		QDomElement fileNameTag = doc.createElement("file");
		mshTag.appendChild(fileNameTag);
		QDomText fileNameText = doc.createTextNode(QString::fromStdString((*it)->getName()));
		fileNameTag.appendChild(fileNameText);
	}

	// STN
	std::vector<std::string> stnNames;
	geoObjects->getStationVectorNames(stnNames);
	for (std::vector<std::string>::const_iterator it(stnNames.begin()); it != stnNames.end(); ++it)
	{
		// write STN file
		XmlStnInterface stn(*geoObjects);
		std::string name(*it);
		stn.setNameForExport(name);

		if (stn.writeToFile(path + name + ".stn"))
		{
			// write entry in project file
			QDomElement geoTag = doc.createElement("stn");
			root.appendChild(geoTag);
			QDomElement fileNameTag = doc.createElement("file");
			geoTag.appendChild(fileNameTag);
			QDomText fileNameText =
			        doc.createTextNode(QString::fromStdString(name + ".stn"));
			fileNameTag.appendChild(fileNameText);
		}
		else
			ERR("XmlGspInterface::writeFile(): Error writing stn-file \"%s\".", name.c_str());
	}

	std::string xml = doc.toString().toStdString();
	_out << xml;
	return true;
}

ProcessVariable::ProcessVariable(BaseLib::ConfigTree const& config,
                                 MeshLib::Mesh& mesh,
                                 GeoLib::GEOObjects const& geometries)
    : _name(config.getConfParam<std::string>("name")),
      _mesh(mesh),
      _n_components(config.getConfParam<int>("components"))
{
	DBUG("Constructing process variable %s", this->_name.c_str());

	// Initial condition
	if (auto ic_config = config.getConfSubtreeOptional("initial_condition"))
	{
		auto const type = ic_config->peekConfParam<std::string>("type");
		if (type == "Uniform")
		{
			_initial_condition =
			    createUniformInitialCondition(*ic_config, _n_components);
		}
		else if (type == "MeshProperty")
		{
			_initial_condition =
			    createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components);
		}
		else
		{
			ERR("Unknown type of the initial condition.");
		}
	}
	else
	{
		INFO("No initial condition found.");
	}

	// Boundary conditions
	if (auto bcs_config = config.getConfSubtreeOptional("boundary_conditions"))
	{
		for (auto bc_config
			 : bcs_config->getConfSubtreeList("boundary_condition"))
		{
			auto const geometrical_set_name =
					bc_config.getConfParam<std::string>("geometrical_set");
			auto const geometry_name =
					bc_config.getConfParam<std::string>("geometry");

			GeoLib::GeoObject const* const geometry =
			    geometries.getGeoObject(geometrical_set_name, geometry_name);
			DBUG(
			    "Found geometry type \"%s\"",
			    GeoLib::convertGeoTypeToString(geometry->getGeoType()).c_str());

			// Construct type dependent boundary condition
			auto const type = bc_config.peekConfParam<std::string>("type");

			if (type == "UniformDirichlet")
			{
				_dirichlet_bc_configs.emplace_back(
				    new UniformDirichletBoundaryCondition(geometry, bc_config));
			}
			else if (type == "UniformNeumann")
			{
				_neumann_bc_configs.emplace_back(
				    new NeumannBcConfig(geometry, bc_config));
			}
			else
			{
				ERR("Unknown type \'%s\' of the boundary condition.",
				    type.c_str());
			}
		}
	} else {
		INFO("No boundary conditions found.");
	}

	// Source Terms
	config.ignoreConfParam("source_terms");
}

	void checkPolylineProperties()
	{
		const GeoLib::PolylineVec *line_vec = geo_objects.getPolylineVecObj(geo_name);
		auto const readerLines = geo_objects.getPolylineVec(geo_name);
		EXPECT_EQ(5u, readerLines->size());

		auto checkPolylineProperty = [this](GeoLib::PolylineVec const* line_vec,
			std::size_t ply_id, std::vector<std::size_t> pnt_ids,
			std::string const& name)
		{
			auto const lines = geo_objects.getPolylineVec(geo_name);
			GeoLib::Polyline* line = (*lines)[ply_id];
			EXPECT_EQ(pnt_ids.size(), line->getNumberOfPoints());
			for (std::size_t k(0); k<pnt_ids.size(); ++k)
				EXPECT_EQ(pnt_ids[k], line->getPointID(k));
			std::string line_name;
			line_vec->getNameOfElementByID(ply_id, line_name);
			EXPECT_EQ(name, line_name);
		};

		checkPolylineProperty(line_vec, 0, {0, 1, 2}, "left");
		checkPolylineProperty(line_vec, 1, {3, 4, 5}, "center");
		checkPolylineProperty(line_vec, 2, {0, 3}, "line1");
		checkPolylineProperty(line_vec, 3, {3, 6}, "line2");
		checkPolylineProperty(line_vec, 4, {6, 7, 8}, "right");
	}