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

void HTProcess::initializeConcreteProcess(
    NumLib::LocalToGlobalIndexMap const& dof_table,
    MeshLib::Mesh const& mesh,
    unsigned const integration_order)
{
    // For the staggered scheme, both processes are assumed to use the same
    // element order. Therefore the order of shape function can be fetched from
    // any set of the sets of process variables of the coupled processes. Here,
    // we take the one from the first process by setting process_id = 0.
    const int process_id = 0;
    ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];

    if (_use_monolithic_scheme)
    {
        ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
            mesh.getDimension(), mesh.getElements(), dof_table,
            pv.getShapeFunctionOrder(), _local_assemblers,
            mesh.isAxiallySymmetric(), integration_order,
            *_material_properties);
    }
    else
    {
        ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
            mesh.getDimension(), mesh.getElements(), dof_table,
            pv.getShapeFunctionOrder(), _local_assemblers,
            mesh.isAxiallySymmetric(), integration_order, *_material_properties,
            _heat_transport_process_id, _hydraulic_process_id);
    }

    _secondary_variables.addSecondaryVariable(
        "darcy_velocity",
        makeExtrapolator(mesh.getDimension(), getExtrapolator(),
                         _local_assemblers,
                         &HTLocalAssemblerInterface::getIntPtDarcyVelocity));
}

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

	TCLAP::CmdLine cmd("Reordering of mesh nodes to make OGS Data Explorer 5 meshes compatible with OGS6.\n" \
	                   "Method 1 is the re-ordering between DataExplorer 5 and DataExplorer 6 meshes,\n" \
	                   "Method 2 is the re-ordering with and without InSitu-Lib in OGS6.",
	                   ' ', "0.1");
	TCLAP::UnlabeledValueArg<std::string> input_mesh_arg("input_mesh",
	                                                     "the name of the input mesh file",
	                                                     true, "", "oldmesh.msh");
	cmd.add(input_mesh_arg);
	TCLAP::UnlabeledValueArg<std::string> output_mesh_arg("output_mesh",
	                                                      "the name of the output mesh file",
	                                                      true, "", "newmesh.vtu");
	cmd.add(output_mesh_arg);
	TCLAP::ValueArg<int> method_arg("m", "method", "reordering method selection", false,  1, "value");

	cmd.add(method_arg);
	cmd.parse(argc, argv);

	MeshLib::Mesh* mesh (FileIO::readMeshFromFile(input_mesh_arg.getValue().c_str()));

	INFO("Reordering nodes... ");
	if (!method_arg.isSet() || method_arg.getValue() == 1)
		reorderNodes(const_cast<std::vector<MeshLib::Element*>&>(mesh->getElements()));
	else if (method_arg.getValue() == 2)
		reorderNodes2(const_cast<std::vector<MeshLib::Element*>&>(mesh->getElements()));
	else
	{
		ERR ("Unknown re-ordering method. Exit program...");
		return 1;
	}

	FileIO::VtuInterface writer(mesh);
	writer.writeToFile(output_mesh_arg.getValue().c_str());

	INFO("VTU file written.");

	delete formatter;
	delete logogCout;
	LOGOG_SHUTDOWN();

	return 0;
}

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 HeatTransportBHEProcess::initializeConcreteProcess(
    NumLib::LocalToGlobalIndexMap const& dof_table,
    MeshLib::Mesh const& mesh,
    unsigned const integration_order)
{
    // Quick access map to BHE's through element ids.
    std::unordered_map<std::size_t, BHE::BHETypes*> element_to_bhe_map;
    int const n_BHEs = _process_data._vec_BHE_property.size();
    for (int i = 0; i < n_BHEs; i++)
    {
        auto const& bhe_elements = _bheMeshData.BHE_elements[i];
        for (auto const& e : bhe_elements)
        {
            element_to_bhe_map[e->getID()] =
                &_process_data._vec_BHE_property[i];
        }
    }

    assert(mesh.getDimension() == 3);
    ProcessLib::HeatTransportBHE::createLocalAssemblers<
        HeatTransportBHELocalAssemblerSoil, HeatTransportBHELocalAssemblerBHE>(
        mesh.getElements(), dof_table, _local_assemblers, element_to_bhe_map,
        mesh.isAxiallySymmetric(), integration_order, _process_data);

    // Create BHE boundary conditions for each of the BHEs
    createBHEBoundaryConditionTopBottom(_bheMeshData.BHE_nodes);
}

std::size_t ElementValueModification::setByElementType(MeshLib::Mesh &mesh, MeshElemType ele_type, int const new_value)
{
    MeshLib::PropertyVector<int>* property_value_vector = nullptr;
    try
    {
        property_value_vector = mesh.getProperties().getPropertyVector<int>(
            "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
    }
    catch (std::runtime_error const& e)
    {
        ERR("%s", e.what());
        return 0;
    }

    std::vector<MeshLib::Element*> const& elements(mesh.getElements());
    std::size_t cnt(0);
    for (std::size_t k(0); k < elements.size(); k++)
    {
        if (elements[k]->getGeomType() != ele_type)
        {
            continue;
        }
        (*property_value_vector)[k] = new_value;
        cnt++;
    }

    return cnt;
}

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

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

std::size_t ElementValueModification::setByElementType(MeshLib::Mesh &mesh, MeshElemType ele_type, int const new_value)
{
    boost::optional<MeshLib::PropertyVector<int> &>
        optional_property_value_vec(
            mesh.getProperties().getPropertyVector<int>("MaterialIDs")
        );

    if (!optional_property_value_vec) {
        return 0;
    }

    MeshLib::PropertyVector<int> & property_value_vector(
        optional_property_value_vec.get()
    );

    std::vector<MeshLib::Element*> const& elements(mesh.getElements());
    std::size_t cnt(0);
    for (std::size_t k(0); k<elements.size(); k++) {
        if (elements[k]->getGeomType()!=ele_type)
            continue;
        property_value_vector[k] = new_value;
        cnt++;
    }

    return cnt;
}

void RichardsComponentTransportProcess::initializeConcreteProcess(
    NumLib::LocalToGlobalIndexMap const& dof_table,
    MeshLib::Mesh const& mesh,
    unsigned const integration_order)
{
    const int monolithic_process_id = 0;
    ProcessLib::ProcessVariable const& pv =
        getProcessVariables(monolithic_process_id)[0];
    ProcessLib::createLocalAssemblers<LocalAssemblerData>(
        mesh.getDimension(), mesh.getElements(), dof_table,
        pv.getShapeFunctionOrder(), _local_assemblers,
        mesh.isAxiallySymmetric(), integration_order, _process_data);

    _secondary_variables.addSecondaryVariable(
        "darcy_velocity",
        makeExtrapolator(mesh.getDimension(), getExtrapolator(),
                         _local_assemblers,
                         &RichardsComponentTransportLocalAssemblerInterface::
                             getIntPtDarcyVelocity));

    _secondary_variables.addSecondaryVariable(
        "saturation",
        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                         &RichardsComponentTransportLocalAssemblerInterface::
                             getIntPtSaturation));
}

bool MeshLayerMapper::createRasterLayers(
    MeshLib::Mesh const& mesh,
    std::vector<GeoLib::Raster const*> const& rasters,
    double minimum_thickness,
    double noDataReplacementValue)
{
    const std::size_t nLayers(rasters.size());
    if (nLayers < 2 || mesh.getDimension() != 2)
    {
        ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input.");
        return false;
    }

    auto top = std::make_unique<MeshLib::Mesh>(mesh);
    if (!layerMapping(*top, *rasters.back(), noDataReplacementValue))
        return false;

    auto bottom = std::make_unique<MeshLib::Mesh>(mesh);
    if (!layerMapping(*bottom, *rasters[0], 0))
    {
        return false;
    }

    this->_minimum_thickness = minimum_thickness;
    std::size_t const nNodes = mesh.getNumberOfNodes();
    _nodes.reserve(nLayers * nNodes);

    // number of triangles in the original mesh
    std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
        [](MeshLib::Element const* elem)
            { return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE);}));
    _elements.reserve(nElems * (nLayers-1));
    _materials.reserve(nElems *  (nLayers-1));

    // add bottom layer
    std::vector<MeshLib::Node*> const& nodes = bottom->getNodes();
    for (MeshLib::Node* node : nodes)
        _nodes.push_back(new MeshLib::Node(*node));

    // add the other layers
    for (std::size_t i=0; i<nLayers-1; ++i)
        addLayerToMesh(*top, i, *rasters[i+1]);

    return true;
}

const std::pair<unsigned, unsigned> MeshInformation::getValueBounds(const MeshLib::Mesh &mesh)
{
	const std::vector<MeshLib::Element*> &elements (mesh.getElements());
	const auto minmax = std::minmax_element(elements.cbegin(), elements.cend(),
        [](MeshLib::Element const*const a, MeshLib::Element const*const b)
            {
                return a->getValue() < b->getValue();
        });
	return std::make_pair<unsigned, unsigned>((*minmax.first)->getValue(), (*minmax.second)->getValue());
}

VolumetricSourceTerm::VolumetricSourceTerm(
    MeshLib::Mesh const& source_term_mesh,
    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
    unsigned const integration_order, unsigned const shapefunction_order,
    Parameter<double> const& volumetric_source_term)
    : SourceTerm(std::move(source_term_dof_table)),
      _volumetric_source_term(volumetric_source_term)
{
    ProcessLib::createLocalAssemblers<VolumetricSourceTermLocalAssembler>(
        source_term_mesh.getDimension(), source_term_mesh.getElements(),
        *_source_term_dof_table, shapefunction_order, _local_assemblers,
        source_term_mesh.isAxiallySymmetric(), integration_order,
        _volumetric_source_term);
}

SurfaceFlux::SurfaceFlux(
    MeshLib::Mesh& boundary_mesh,
    std::size_t bulk_property_number_of_components,
    unsigned const integration_order)
{
    DBUG("Create local balance assemblers.");
    // Populate the vector of local assemblers.
    _local_assemblers.resize(boundary_mesh.getElements().size());

    // needed to create dof table
    auto mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
        boundary_mesh, boundary_mesh.getNodes());

    // Collect the mesh subsets in a vector.
    std::vector<MeshLib::MeshSubset> all_mesh_subsets;
    std::generate_n(std::back_inserter(all_mesh_subsets),
                    bulk_property_number_of_components,
                    [&]() { return *mesh_subset_all_nodes; });

    // needed for creation of local assemblers
    auto dof_table = std::make_unique<NumLib::LocalToGlobalIndexMap const>(
        std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION);

    auto const bulk_element_ids =
        boundary_mesh.getProperties().template getPropertyVector<std::size_t>(
            "bulk_element_ids", MeshLib::MeshItemType::Cell, 1);
    auto const bulk_face_ids =
        boundary_mesh.getProperties().template getPropertyVector<std::size_t>(
            "bulk_face_ids", MeshLib::MeshItemType::Cell, 1);

    ProcessLib::createLocalAssemblers<SurfaceFluxLocalAssembler>(
        boundary_mesh.getDimension() + 1,  // or bulk_mesh.getDimension()?
        boundary_mesh.getElements(), *dof_table, 1, _local_assemblers,
        boundary_mesh.isAxiallySymmetric(), integration_order,
        *bulk_element_ids, *bulk_face_ids);
}

void
detectHoles(MeshLib::Mesh const& mesh,
	std::vector<std::size_t> erase_elems,
	std::size_t const expected_n_holes)
{
	std::vector<MeshLib::Node*> nodes = MeshLib::copyNodeVector(mesh.getNodes());
	std::vector<MeshLib::Element*> elems = MeshLib::copyElementVector(mesh.getElements(),nodes);
	for (auto pos : erase_elems)
	{
		delete elems[pos];
		elems.erase(elems.begin()+pos);
	}
	MeshLib::Mesh mesh2("mesh2", nodes, elems);
	ASSERT_EQ(expected_n_holes, MeshLib::MeshValidation::detectHoles(mesh2));
};