C++ Mesh::getNNodes方法代码示例

std::vector<GeoLib::PointWithID*> MshEditor::getSurfaceNodes(const MeshLib::Mesh &mesh, const double *dir)
{
	INFO ("Extracting surface nodes...");
	const std::vector<MeshLib::Element*> all_elements (mesh.getElements());
	const std::vector<MeshLib::Node*> all_nodes (mesh.getNodes());

	std::vector<MeshLib::Element*> sfc_elements;
	get2DSurfaceElements(all_elements, sfc_elements, dir, mesh.getDimension());

	std::vector<MeshLib::Node*> sfc_nodes;
	std::vector<unsigned> node_id_map(mesh.getNNodes());
	get2DSurfaceNodes(all_nodes, sfc_nodes, sfc_elements, node_id_map);

	const unsigned nElements (sfc_elements.size());
	for (unsigned i=0; i<nElements; ++i)
		delete sfc_elements[i];

	const size_t nNodes (sfc_nodes.size());
	std::vector<GeoLib::PointWithID*> surface_pnts(nNodes);
	for (unsigned i=0; i<nNodes; ++i)
	{
		surface_pnts[i] = new GeoLib::PointWithID(sfc_nodes[i]->getCoords(), sfc_nodes[i]->getID());
		delete sfc_nodes[i];
	}
	return surface_pnts;
}

void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster)
{
	const std::size_t nNodes (dem_mesh.getNNodes());
	const std::vector<MeshLib::Node*> &nodes (dem_mesh.getNodes());
	const std::size_t node_id_offset (_nodes.size());
	const std::size_t last_layer_node_offset (node_id_offset-nNodes);

	for (std::size_t i=0; i<nNodes; ++i)
		_nodes.push_back(getNewLayerNode(*nodes[i], *_nodes[last_layer_node_offset + i], raster, _nodes.size()));

	const std::vector<MeshLib::Element*> &layer_elements (dem_mesh.getElements());
	for (MeshLib::Element* elem : layer_elements)
	{
		if (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE)
		{
			std::array<MeshLib::Node*,3> tri_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)],
			                                            _nodes[node_id_offset+elem->getNodeIndex(1)],
			                                            _nodes[node_id_offset+elem->getNodeIndex(2)] }};
			_elements.push_back(new MeshLib::Tri(tri_nodes, layer_id));
		}
		else if (elem->getGeomType() == MeshLib::MeshElemType::QUAD)
		{
			std::array<MeshLib::Node*,4> quad_nodes = {{ _nodes[node_id_offset+elem->getNodeIndex(0)],
			                                             _nodes[node_id_offset+elem->getNodeIndex(1)],
			                                             _nodes[node_id_offset+elem->getNodeIndex(2)],
			                                             _nodes[node_id_offset+elem->getNodeIndex(3)] }};
			_elements.push_back(new MeshLib::Quad(quad_nodes, layer_id));
		}
	}
}

int main (int argc, char* argv[])
{
	LOGOG_INITIALIZE();
	logog::Cout* logog_cout (new logog::Cout);
	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
	logog_cout->SetFormatter(*custom_format);

	TCLAP::CmdLine cmd("Converts VTK mesh into OGS mesh.", ' ', "0.1");
	TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
	                                     "the name of the file containing the input mesh", true,
	                                     "", "file name of input mesh");
	cmd.add(mesh_in);
	TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
	                                      "the name of the file the mesh will be written to", true,
	                                      "", "file name of output mesh");
	cmd.add(mesh_out);
	cmd.parse(argc, argv);

	MeshLib::Mesh* mesh (FileIO::BoostVtuInterface::readVTUFile(mesh_in.getValue()));
	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());

	FileIO::Legacy::MeshIO meshIO;
	meshIO.setMesh(mesh);
	meshIO.writeToFile(mesh_out.getValue());

	delete custom_format;
	delete logog_cout;
	LOGOG_SHUTDOWN();

	return 0;
}

void MeshSurfaceExtraction::getSurfaceAreaForNodes(const MeshLib::Mesh &mesh, std::vector<double> &node_area_vec)
{
	if (mesh.getDimension() == 2)
	{
		double total_area (0);

		// for each node, a vector containing all the element idget every element
		const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
		const size_t nNodes ( mesh.getNNodes() );
		node_area_vec.reserve(nNodes);
		for (size_t n=0; n<nNodes; ++n)
		{
			double node_area (0);

			std::vector<MeshLib::Element*> conn_elems = nodes[n]->getElements();
			const size_t nConnElems (conn_elems.size());

			for (size_t i=0; i<nConnElems; ++i)
			{
				const MeshLib::Element* elem (conn_elems[i]);
				const unsigned nElemParts = (elem->getGeomType() == MeshElemType::TRIANGLE) ? 3 : 4;
				const double area = conn_elems[i]->getContent() / nElemParts;
				node_area += area;
				total_area += area;
			}

			node_area_vec.push_back(node_area);
		}

		INFO ("Total surface Area: %f", total_area);
	}
	else
		ERR ("Error in MeshSurfaceExtraction::getSurfaceAreaForNodes() - Given mesh is no surface mesh (dimension != 2).");
}

void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers)
{
	const unsigned nLayerBoundaries (nLayers-1);
	const std::size_t nNodes (layer.getNNodes());
	const std::vector<MeshLib::Element*> &layer_elements (layer.getElements());
	for (MeshLib::Element* elem : layer_elements)
	{
		const std::size_t nElemNodes (elem->getNBaseNodes());
		for (unsigned i=0; i<nElemNodes; ++i)
			if (elem->getNeighbor(i) == nullptr)
				for (unsigned j=0; j<nLayerBoundaries; ++j)
				{
					const std::size_t offset (j*nNodes);
					MeshLib::Node* n0 = _nodes[offset + elem->getNodeIndex(i)];
					MeshLib::Node* n1 = _nodes[offset + elem->getNodeIndex((i+1)%nElemNodes)];
					MeshLib::Node* n2 = _nodes[offset + nNodes + elem->getNodeIndex((i+1)%nElemNodes)];
					MeshLib::Node* n3 = _nodes[offset + nNodes + elem->getNodeIndex(i)];

					if (MathLib::Vector3(*n1, *n2).getLength() > std::numeric_limits<double>::epsilon())
					{
						const std::array<MeshLib::Node*,3> tri_nodes = {{ n0, n2, n1 }};
						_elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+1+j));
					}
					if (MathLib::Vector3(*n0, *n3).getLength() > std::numeric_limits<double>::epsilon())
					{
						const std::array<MeshLib::Node*,3> tri_nodes = {{ n0, n3, n2 }};
						_elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+1+j));
					}
				}
	}
}

	void testZCoords2D(MeshLib::Mesh const& input, MeshLib::Mesh const& output, double height)
	{
		std::size_t const nNodes (input.getNNodes());
		for (std::size_t i=0; i<nNodes; ++i)
		{
			ASSERT_EQ((*input.getNode(i))[2], (*output.getNode(i))[2]);
			ASSERT_EQ((*input.getNode(i))[2] + height, (*output.getNode(nNodes+i))[2]);
		}
	}

std::vector<GeoLib::Point*> MeshSurfaceExtraction::getSurfaceNodes(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, double angle)
{
	INFO ("Extracting surface nodes...");
	std::vector<MeshLib::Element*> sfc_elements;
	get2DSurfaceElements(mesh.getElements(), sfc_elements, dir, angle, mesh.getDimension());

	std::vector<MeshLib::Node*> sfc_nodes;
	std::vector<std::size_t> node_id_map(mesh.getNNodes());
	get2DSurfaceNodes(sfc_nodes, mesh.getNNodes(), sfc_elements, node_id_map);

	for (auto e : sfc_elements)
		delete e;

	const std::size_t nNodes (sfc_nodes.size());
	std::vector<GeoLib::Point*> surface_pnts(nNodes);
	for (std::size_t i=0; i<nNodes; ++i)
	{
		surface_pnts[i] = new GeoLib::Point(*(sfc_nodes[i]), sfc_nodes[i]->getID());
		delete sfc_nodes[i];
	}
	return surface_pnts;
}

SparsityPattern
computeSparsityPattern(LocalToGlobalIndexMap const& dof_table,
        MeshLib::Mesh const& mesh
        )
{
    MeshLib::NodeAdjacencyTable node_adjacency_table;
    node_adjacency_table.createTable(mesh.getNodes());

    // A mapping   mesh node id -> global indices
    // It acts as a cache for dof table queries.
    std::vector<std::vector<GlobalIndexType> > global_idcs;

    global_idcs.reserve(mesh.getNNodes());
    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
    {
        MeshLib::Location l(mesh.getID(), MeshLib::MeshItemType::Node, n);
        global_idcs.push_back(dof_table.getGlobalIndices(l));
    }

    SparsityPattern sparsity_pattern(dof_table.dofSize());

    // Map adjacent mesh nodes to "adjacent global indices".
    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
    {
        auto const& node_ids = node_adjacency_table.getAdjacentNodes(n);
        for (auto an : node_ids) {
            auto const& row_ids = global_idcs[an];
            auto const num_components = row_ids.size();
            for (auto r : row_ids) {
                // Each component leads to an entry in the row.
                // For the sparsity pattern only the number of entries are needed.
                sparsity_pattern[r] += num_components;
            }
        }
    }

    return sparsity_pattern;
}

TEST(MeshLib, Duplicate)
{
	MeshLib::Mesh* mesh (MeshLib::MeshGenerator::generateRegularQuadMesh(10, 5, 1));

	std::vector<MeshLib::Node*> new_nodes (MeshLib::copyNodeVector(mesh->getNodes()));
	std::vector<MeshLib::Element*> new_elements (MeshLib::copyElementVector(mesh->getElements(), new_nodes));

	MeshLib::Mesh new_mesh ("new", new_nodes, new_elements);

	ASSERT_EQ (mesh->getNElements(), new_mesh.getNElements());
	ASSERT_EQ (mesh->getNNodes(), new_mesh.getNNodes());

	std::vector<std::size_t> del_idx(1,1);
	MeshLib::removeMeshNodes(*mesh, del_idx);

	ASSERT_EQ (mesh->getNElements(), new_mesh.getNElements()-2);
	ASSERT_EQ (mesh->getNNodes(), new_mesh.getNNodes()-2);

	ASSERT_DOUBLE_EQ (4.0, MathLib::sqrDist(*mesh->getNode(0), *new_mesh.getNode(0)));
	ASSERT_DOUBLE_EQ (0.0, MathLib::sqrDist(*mesh->getNode(0), *new_mesh.getNode(2)));

	ASSERT_DOUBLE_EQ (4.0, MathLib::sqrDist(*mesh->getElement(0)->getNode(0), *new_mesh.getElement(0)->getNode(0)));
	ASSERT_DOUBLE_EQ (0.0, MathLib::sqrDist(*mesh->getElement(0)->getNode(0), *new_mesh.getElement(2)->getNode(0)));
}

TEST(MeshLib, ElementStatus)
{
	const unsigned width (100);
	const unsigned elements_per_side (20);
	MeshLib::Mesh* mesh (MeshLib::MeshGenerator::generateRegularQuadMesh(width, elements_per_side));
	MeshLib::ElementStatus status(mesh);
	const std::vector<MeshLib::Element*> elements (mesh->getElements());

	for (unsigned i=0; i<elements_per_side; ++i)
	{
		for (unsigned j=0; j<elements_per_side; ++j)
			elements[i*elements_per_side + j]->setValue(i);
	}

	// all elements active
	ASSERT_EQ (elements.size(), status.getNActiveElements());
	// all nodes active
	ASSERT_EQ (mesh->getNNodes(), status.getNActiveNodes());

	// set material 1 to false
	status.setMaterialStatus(1, false);
	ASSERT_EQ (elements.size()-elements_per_side, status.getNActiveElements());

	// set material 1 to false (again)
	status.setMaterialStatus(1, false);
	ASSERT_EQ (elements.size()-elements_per_side, status.getNActiveElements());

	// set material 0 to false
	status.setMaterialStatus(0, false);
	ASSERT_EQ (elements.size()-(2*elements_per_side), status.getNActiveElements());

	// active elements
	std::vector<std::size_t> active_elements (status.getActiveElements());
	ASSERT_EQ (active_elements.size(), status.getNActiveElements());

	// active nodes
	std::vector<std::size_t> active_nodes (status.getActiveNodes());
	ASSERT_EQ (active_nodes.size(), status.getNActiveNodes());

	// set element 1 to false (yet again)
	status.setElementStatus(1, false);
	status.getElementStatus(1);
	ASSERT_EQ (elements.size()-(2*elements_per_side), status.getNActiveElements());
	ASSERT_EQ (mesh->getNNodes()-(2*(elements_per_side+1)), status.getNActiveNodes());

	// set element 1 to true
	status.setElementStatus(1, true);
	ASSERT_EQ (elements.size()-(2*elements_per_side)+1, status.getNActiveElements());
	ASSERT_EQ (mesh->getNNodes()-(2*(elements_per_side+1))+4, status.getNActiveNodes());
	ASSERT_EQ(status.getElementStatus(1), true);

	std::vector<std::size_t> active_elements_at_node (status.getActiveElementsAtNode(2));
	ASSERT_EQ(1u, active_elements_at_node.size());
	active_elements_at_node = status.getActiveElementsAtNode(22);
	ASSERT_EQ(1u, active_elements_at_node.size());
	active_elements_at_node = status.getActiveElementsAtNode(44);
	ASSERT_EQ(2u, active_elements_at_node.size());
	active_elements_at_node = status.getActiveElementsAtNode(102);
	ASSERT_EQ(4u, active_elements_at_node.size());

	status.setAll(true);
	ASSERT_EQ(elements.size(), status.getNActiveElements());
	ASSERT_EQ(mesh->getNNodes(), status.getNActiveNodes());

	status.setAll(false);
	ASSERT_EQ(0u, status.getNActiveElements());
	ASSERT_EQ(0u, status.getNActiveNodes());
}

int main (int argc, char* argv[])
{
	LOGOG_INITIALIZE();
	logog::Cout* logog_cout (new logog::Cout);
	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
	logog_cout->SetFormatter(*custom_format);

	TCLAP::CmdLine cmd("Add EMI data as a scalar cell array to a 2d mesh.", ' ', "0.1");

	// I/O params
	TCLAP::ValueArg<std::string> poly_out("o", "polydata-output-file",
	                                      "the name of the file the data will be written to", true,
	                                      "", "file name of polydata file");
	cmd.add(poly_out);
	TCLAP::ValueArg<std::string> csv_in("i", "csv-input-file",
	                                    "csv-file containing EMI data", true,
	                                    "", "name of the csv input file");
	cmd.add(csv_in);
	TCLAP::ValueArg<std::string> dem_in("s", "DEM-file",
	                                    "Surface DEM for mapping ERT data", false,
	                                    "", "file name of the Surface DEM");
	cmd.add(dem_in);
	cmd.parse(argc, argv);

	MeshLib::Mesh* mesh (nullptr);
	if (dem_in.isSet())
	{
		mesh = FileIO::VtuInterface::readVTUFile(dem_in.getValue());
		if (mesh == nullptr)
		{
			ERR ("Error reading mesh file.");
			return -2;
		}

		if (mesh->getDimension() != 2)
		{
			ERR ("This utility can handle only 2d meshes at this point.");
			delete mesh;
			return -3;
		}
		INFO("Surface mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());
	}

	GeoLib::GEOObjects geo_objects;
	FileIO::XmlGmlInterface xml(geo_objects);
	//std::vector<GeoLib::Polyline*> *lines = new std::vector<GeoLib::Polyline*>;
	std::array<char, 2> dipol = {{ 'H', 'V' }};
	std::array<char,3> const regions = {{'A', 'B', 'C'}};
	for (std::size_t j=0; j<dipol.size(); ++j)
	{
		std::vector<GeoLib::Point*> *points   = new std::vector<GeoLib::Point*>;
		for (std::size_t i=0; i<regions.size(); ++i)
		{
			//std::size_t const start_idx (points->size());
			getPointsFromFile(*points, csv_in.getValue(), dipol[j], regions[i]);
			//std::size_t const end_idx (points->size());
			//GeoLib::Polyline* line = new GeoLib::Polyline(*points);
			//for (std::size_t j=start_idx; j<end_idx; ++j)
			//	line->addPoint(j);
			//lines->push_back(line);
		}
		std::string geo_name (std::string("EMI Data ").append(1,dipol[j]));
		geo_objects.addPointVec(points, geo_name);
		//geo_objects.addPolylineVec(lines, geo_name);

		if (mesh != nullptr)
		{
			GeoMapper mapper(geo_objects, geo_name);
			mapper.mapOnMesh(mesh);
		}
		
		xml.setNameForExport(geo_name);
		std::string const output_name = poly_out.getValue() + "_" + dipol[j] + ".gml";
		xml.writeToFile(output_name);
		
		std::vector<double> emi;
		for (std::size_t i=0; i<regions.size(); ++i)
			getMeasurements(emi, csv_in.getValue(), dipol[j], regions[i]);
		writeMeasurementsToFile(emi, poly_out.getValue(), dipol[j]);
		std::for_each(points->begin(), points->end(), std::default_delete<GeoLib::Point>());
		delete points;
	}
	
	delete mesh;
	delete custom_format;
	delete logog_cout;
	LOGOG_SHUTDOWN();

	return 0;
}

int main (int argc, char* argv[])
{
	LOGOG_INITIALIZE();
	logog::Cout* logog_cout (new logog::Cout);
	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
	logog_cout->SetFormatter(*custom_format);

	TCLAP::CmdLine cmd("Edit material IDs of mesh elements.", ' ', "0.1");
	TCLAP::SwitchArg replaceArg("r", "replace", "replace material IDs", false);
	TCLAP::SwitchArg condenseArg("c", "condense", "condense material IDs", false);
	TCLAP::SwitchArg specifyArg("s", "specify", "specify material IDs by element types (-e)", false);
	std::vector<TCLAP::Arg*> vec_xors;
	vec_xors.push_back(&replaceArg);
	vec_xors.push_back(&condenseArg);
	vec_xors.push_back(&specifyArg);
	cmd.xorAdd(vec_xors);
	TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
	                                     "the name of the file containing the input mesh", true,
	                                     "", "file name");
	cmd.add(mesh_in);
	TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
	                                      "the name of the file the mesh will be written to", true,
	                                      "", "file name");
	cmd.add(mesh_out);
	TCLAP::MultiArg<unsigned> matIDArg("m", "current-material-id",
	                                      "current material id to be replaced", false, "number");
	cmd.add(matIDArg);
	TCLAP::ValueArg<unsigned> newIDArg("n", "new-material-id",
	                                      "new material id", false, 0, "number");
	cmd.add(newIDArg);
	std::vector<std::string> eleList(MeshLib::getMeshElemTypeStringsShort());
	TCLAP::ValuesConstraint<std::string> allowedVals(eleList);
	TCLAP::ValueArg<std::string> eleTypeArg("e", "element-type",
	                                      "element type", false, "", &allowedVals);
	cmd.add(eleTypeArg);

	cmd.parse(argc, argv);

	if (!replaceArg.isSet() && !condenseArg.isSet() && !specifyArg.isSet()) {
		INFO("Please select editing mode: -r or -c or -s");
		return 0;
	} else if (replaceArg.isSet() && condenseArg.isSet()) {
		INFO("Please select only one editing mode: -r or -c or -s");
		return 0;
	} else if (replaceArg.isSet()) {
		if (!matIDArg.isSet() || !newIDArg.isSet()) {
			INFO("current and new material IDs must be provided for replacement");
			return 0;
		}
	} else if (specifyArg.isSet()) {
		if (!eleTypeArg.isSet() || !newIDArg.isSet()) {
			INFO("element type and new material IDs must be provided to specify elements");
			return 0;
		}
	}

	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());

	if (condenseArg.isSet()) {
		INFO("Condensing material ID...");
		MeshLib::ElementValueModification::condense(*mesh);
	} else if (replaceArg.isSet()) {
		INFO("Replacing material ID...");
		const auto vecOldID = matIDArg.getValue();
		const unsigned newID = newIDArg.getValue();
		for (auto oldID : vecOldID) {
			INFO("%d -> %d", oldID, newID);
			MeshLib::ElementValueModification::replace(*mesh, oldID, newID, true);
		}
	} else if (specifyArg.isSet()) {
		INFO("Specifying material ID...");
		const std::string eleTypeName(eleTypeArg.getValue());
		const MeshLib::MeshElemType eleType = MeshLib::String2MeshElemType(eleTypeName);
		const unsigned newID = newIDArg.getValue();
		unsigned cnt = MeshLib::ElementValueModification::setByElementType(*mesh, eleType, newID);
		INFO("updated %d elements", cnt);
	}

	// write into a file
	FileIO::writeMeshToFile(*mesh, mesh_out.getValue());

	delete custom_format;
	delete logog_cout;
	LOGOG_SHUTDOWN();

	return 0;
}

int main (int argc, char* argv[])
{
	LOGOG_INITIALIZE();
	logog::Cout* logog_cout (new logog::Cout);
	BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
	logog_cout->SetFormatter(*custom_format);

	TCLAP::CmdLine cmd("Remove mesh elements.", ' ', "0.1");
	TCLAP::ValueArg<std::string> mesh_in("i", "mesh-input-file",
	                                     "the name of the file containing the input mesh", true,
	                                     "", "file name of input mesh");
	cmd.add(mesh_in);
	TCLAP::ValueArg<std::string> mesh_out("o", "mesh-output-file",
	                                      "the name of the file the mesh will be written to", true,
	                                      "", "file name of output mesh");
	cmd.add(mesh_out);
	TCLAP::SwitchArg zveArg("z", "zero-volume", "remove zero volume elements", false);
	cmd.add(zveArg);
	TCLAP::MultiArg<std::string> eleTypeArg("t", "element-type",
	                                      "element type to be removed", false, "element type");
	cmd.add(eleTypeArg);
	TCLAP::MultiArg<unsigned> matIDArg("m", "material-id",
	                                      "material id", false, "material id");
	cmd.add(matIDArg);
	cmd.parse(argc, argv);

	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(mesh_in.getValue()));
	INFO("Mesh read: %d nodes, %d elements.", mesh->getNNodes(), mesh->getNElements());

	// search elements IDs to be removed
	std::vector<std::size_t> vec_elementIDs_removed;
	if (zveArg.isSet()) {
		std::vector<std::size_t> vec_matched = searchByZeroContent(mesh->getElements());
		updateUnion(vec_matched, vec_elementIDs_removed);
		INFO("%d zero volume elements found.", vec_matched.size());
	}
	if (eleTypeArg.isSet()) {
		std::vector<std::string> eleTypeNames = eleTypeArg.getValue();
		for (auto typeName : eleTypeNames) {
			MeshElemType type = String2MeshElemType(typeName);
			if (type == MeshElemType::INVALID) continue;
			std::vector<std::size_t> vec_matched = searchByElementType(mesh->getElements(), type);
			updateUnion(vec_matched, vec_elementIDs_removed);
			INFO("%d %s elements found.", vec_matched.size(), typeName.c_str());
		}
	}
	if (matIDArg.isSet()) {
		std::vector<unsigned> vec_matID = matIDArg.getValue();
		for (auto matID : vec_matID) {
			std::vector<std::size_t> vec_matched = searchByMaterialID(mesh->getElements(), matID);
			updateUnion(vec_matched, vec_elementIDs_removed);
			INFO("%d elements with material ID %d found.", vec_matched.size(), matID);
		}
	}

	// remove the elements
	INFO("Removing total %d elements...", vec_elementIDs_removed.size());
	std::vector<MeshLib::Element*> tmp_eles = excludeElements(mesh->getElements(), vec_elementIDs_removed);
	INFO("%d elements remained.", tmp_eles.size());
	std::vector<MeshLib::Node*> new_nodes;
	std::vector<MeshLib::Element*> new_eles;
	copyNodesElements(mesh->getNodes(), tmp_eles, new_nodes, new_eles);

	// create a new mesh object. Unsued nodes are removed while construction
	MeshLib::Mesh* new_mesh(new MeshLib::Mesh(mesh->getName(), new_nodes, new_eles));

	// write into a file
	FileIO::Legacy::MeshIO meshIO;
	meshIO.setMesh(new_mesh);
	meshIO.writeToFile(mesh_out.getValue());

	delete custom_format;
	delete logog_cout;
	LOGOG_SHUTDOWN();

	return 0;
}