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

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

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

GlobalSparsityPattern computeSparsityPatternNonPETSc(
    NumLib::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.getNumberOfNodes());
    for (std::size_t n = 0; n < mesh.getNumberOfNodes(); ++n)
    {
        MeshLib::Location l(mesh.getID(), MeshLib::MeshItemType::Node, n);
        global_idcs.push_back(dof_table.getGlobalIndices(l));
    }

    GlobalSparsityPattern sparsity_pattern(dof_table.dofSizeWithGhosts());

    // Map adjacent mesh nodes to "adjacent global indices".
    for (std::size_t n = 0; n < mesh.getNumberOfNodes(); ++n)
    {
        unsigned n_connected_dof = 0;
        for (auto an : node_adjacency_table.getAdjacentNodes(n))
            n_connected_dof += global_idcs[an].size();
        for (auto global_index : global_idcs[n])
            sparsity_pattern[global_index] = n_connected_dof;
    }

    return sparsity_pattern;
}

bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const& raster, double noDataReplacementValue = 0.0)
{
    if (new_mesh.getDimension() != 2)
    {
        ERR("MshLayerMapper::layerMapping() - requires 2D mesh");
        return false;
    }

    GeoLib::RasterHeader const& header (raster.getHeader());
    const double x0(header.origin[0]);
    const double y0(header.origin[1]);
    const double delta(header.cell_size);

    const std::pair<double, double> xDim(x0, x0 + header.n_cols * delta); // extension in x-dimension
    const std::pair<double, double> yDim(y0, y0 + header.n_rows * delta); // extension in y-dimension

    const std::size_t nNodes (new_mesh.getNumberOfNodes());
    const std::vector<MeshLib::Node*> &nodes = new_mesh.getNodes();
    for (unsigned i = 0; i < nNodes; ++i)
    {
        if (!raster.isPntOnRaster(*nodes[i]))
        {
            // use either default value or elevation from layer above
            nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], noDataReplacementValue);
            continue;
        }

        double elevation (raster.interpolateValueAtPoint(*nodes[i]));
        if (std::abs(elevation - header.no_data) < std::numeric_limits<double>::epsilon())
            elevation = noDataReplacementValue;
        nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], elevation);
    }

    return true;
}

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

bool MeshLayerMapper::mapToStaticValue(MeshLib::Mesh &mesh, double value)
{
    if (mesh.getDimension() != 2)
    {
        ERR("MshLayerMapper::mapToStaticValue() - requires 2D mesh");
        return false;
    }

    std::vector<MeshLib::Node*> const& nodes (mesh.getNodes());
    for (MeshLib::Node* node : nodes)
        node->updateCoordinates((*node)[0], (*node)[1], value);
    return true;
}

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

DirichletBoundaryConditionWithinTimeInterval::
    DirichletBoundaryConditionWithinTimeInterval(
        std::unique_ptr<BaseLib::TimeInterval> time_interval,
        Parameter<double> const& parameter, MeshLib::Mesh const& bc_mesh,
        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
        int const variable_id, int const component_id)
    : _parameter(parameter),
      _bc_mesh(bc_mesh),
      _nodes_in_bc_mesh(bc_mesh.getNodes()),
      _variable_id(variable_id),
      _component_id(component_id),
      _time_interval(std::move(time_interval))
{
    config(dof_table_bulk);
}

bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh,
                                       const std::vector<GeoLib::Raster const*> &rasters,
                                       double minimum_thickness,
                                       double noDataReplacementValue)
{
	if (mesh.getDimension() != 2)
		return false;

	_elevation_epsilon = calcEpsilon(*rasters[0], *rasters.back());
	if (_elevation_epsilon <= 0)
		return false;

	// remove line elements, only tri + quad remain
	MeshLib::ElementSearch ex(mesh);
	ex.searchByElementType(MeshLib::MeshElemType::LINE);
	MeshLib::Mesh* top (removeElements(mesh, ex.getSearchedElementIDs(), "MeshLayer"));
	if (top==nullptr)
		top = new MeshLib::Mesh(mesh);

	if (!MeshLib::MeshLayerMapper::layerMapping(*top, *rasters.back(), noDataReplacementValue))
		return false;

	MeshLib::Mesh* bottom (new MeshLib::Mesh(*top));
	if (!MeshLib::MeshLayerMapper::layerMapping(*bottom, *rasters[0], 0))
	{
		delete top;
		return false;
	}

	this->_minimum_thickness = minimum_thickness;
	_nodes = MeshLib::copyNodeVector(bottom->getNodes());
	_elements = MeshLib::copyElementVector(bottom->getElements(), _nodes);
	delete bottom;

	// map each layer and attach to subsurface mesh
	const std::size_t nRasters (rasters.size());
	for (std::size_t i=1; i<nRasters; ++i)
		this->addLayerToMesh(*top, i, *rasters[i]);

	// close boundaries between layers
	this->addLayerBoundaries(*top, nRasters);
	this->removeCongruentElements(nRasters, top->getNElements());
	delete top;
	return true;
}

/// Returns a vector of values where each value is associated with a
/// particular node. Since a node is part of elements, it is possible to
/// assign an area per element to this node. Each value of the return vector
/// is the sum of the assigned area (per element) multiplied by the property
/// value of the element.
/// @param mesh a surface mesh containing a property \c prop_name assigned
/// to cells
/// @param prop_name name of the cell based property within the \c mesh
/// @return vector of integration values associated to the surface mesh nodes
std::vector<double> getSurfaceIntegratedValuesForNodes(
    const MeshLib::Mesh& mesh, std::string const& prop_name)
{
    if (mesh.getDimension() != 2)
    {
        ERR("Error in "
            "MeshSurfaceExtraction::getSurfaceIntegratedValuesForNodes() - "
            "Given mesh is no surface mesh (dimension != 2).");
        return std::vector<double>();
    }

    if (!mesh.getProperties().existsPropertyVector<double>(prop_name))
    {
        ERR("Need element property, but the property '%s' is not "
            "available.",
            prop_name.c_str());
        return std::vector<double>();
    }
    auto const* const elem_pv = mesh.getProperties().getPropertyVector<double>(
        prop_name, MeshLib::MeshItemType::Cell, 1);

    std::vector<double> integrated_node_area_vec;
    double total_area(0);

    for (auto const* node : mesh.getNodes())
    {
        double node_area(0);
        double integrated_node_area(0);
        for (auto const& connected_elem : node->getElements())
        {
            double const area = connected_elem->getContent() /
                                connected_elem->getNumberOfBaseNodes();
            node_area += area;
            integrated_node_area += area * (*elem_pv)[connected_elem->getID()];
            total_area += area;
        }

        integrated_node_area_vec.push_back(integrated_node_area);
    }

    INFO ("Total surface area: %g", total_area);

    return integrated_node_area_vec;
}

    ExtrapolationTestProcess(MeshLib::Mesh const& mesh,
                             unsigned const integration_order)
        : _integration_order(integration_order),
          _mesh_subset_all_nodes(mesh, mesh.getNodes())
    {
        std::vector<MeshLib::MeshSubset> all_mesh_subsets{
            _mesh_subset_all_nodes};

        _dof_table = std::make_unique<NumLib::LocalToGlobalIndexMap>(
            std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_COMPONENT);

        // Passing _dof_table works, because this process has only one variable
        // and the variable has exactly one component.
        _extrapolator =
            std::make_unique<ExtrapolatorImplementation>(*_dof_table);

        // createAssemblers(mesh);
        ProcessLib::createLocalAssemblers<LocalAssemblerData>(
            mesh.getDimension(), mesh.getElements(), *_dof_table, 1,
            _local_assemblers, mesh.isAxiallySymmetric(), _integration_order);
    }

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

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