C++ meshlib::Mesh类代码示例

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

std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
    NumLib::LocalToGlobalIndexMap const& dof_table, std::size_t bulk_mesh_id,
    int const variable_id, int const component_id,
    unsigned const integration_order, unsigned const shapefunction_order,
    unsigned const global_dim)
{
    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
    config.checkConfigParameter("type", "Python");

    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Python__bc_object}
    auto const bc_object = config.getConfigParameter<std::string>("bc_object");
    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Python__flush_stdout}
    auto const flush_stdout = config.getConfigParameter("flush_stdout", false);

    // Evaluate Python code in scope of main module
    pybind11::object scope =
        pybind11::module::import("__main__").attr("__dict__");

    if (!scope.contains(bc_object))
        OGS_FATAL(
            "Function `%s' is not defined in the python script file, or there "
            "was no python script file specified.",
            bc_object.c_str());

    auto* bc = scope[bc_object.c_str()]
                   .cast<PythonBoundaryConditionPythonSideInterface*>();

    if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
        component_id >= dof_table.getNumberOfVariableComponents(variable_id))
    {
        OGS_FATAL(
            "Variable id or component id too high. Actual values: (%d, %d), "
            "maximum values: (%d, %d).",
            variable_id, component_id, dof_table.getNumberOfVariables(),
            dof_table.getNumberOfVariableComponents(variable_id));
    }

    // In case of partitioned mesh the boundary could be empty, i.e. there is no
    // boundary condition.
#ifdef USE_PETSC
    // This can be extracted to createBoundaryCondition() but then the config
    // parameters are not read and will cause an error.
    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
    // subtree and move the code up in createBoundaryCondition().
    if (boundary_mesh.getDimension() == 0 &&
        boundary_mesh.getNumberOfNodes() == 0 &&
        boundary_mesh.getNumberOfElements() == 0)
    {
        return nullptr;
    }
#endif  // USE_PETSC

    return std::make_unique<PythonBoundaryCondition>(
        PythonBoundaryConditionData{
            bc, dof_table, bulk_mesh_id,
            dof_table.getGlobalComponent(variable_id, component_id),
            boundary_mesh},
        integration_order, shapefunction_order, global_dim, flush_stdout);
}

IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
                                                     std::string const& name)
{
    if (!mesh.getProperties().existsPropertyVector<char>(
            "IntegrationPointMetaData"))
    {
        OGS_FATAL(
            "Integration point data '%s' is present in the vtk field "
            "data but the required 'IntegrationPointMetaData' array "
            "is not available.",
            name.c_str());
    }
    auto const& mesh_property_ip_meta_data =
        *mesh.getProperties().template getPropertyVector<char>(
            "IntegrationPointMetaData");

    if (mesh_property_ip_meta_data.getMeshItemType() !=
        MeshLib::MeshItemType::IntegrationPoint)
    {
        OGS_FATAL("IntegrationPointMetaData array must be field data.");
    }

    // Find the current integration point data entry and extract the
    // meta data.
    auto const ip_meta_data = extractIntegrationPointMetaData(
        json::parse(mesh_property_ip_meta_data.begin(),
                    mesh_property_ip_meta_data.end()),
        name);

    return ip_meta_data;
}

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

// Creates a PropertyVector<double> and maps it into a vtkDataArray-equivalent
TEST(InSituLibMappedPropertyVector, Int)
{
	const size_t mesh_size = 5;
	const double length = 1.0;

	MeshLib::Mesh* mesh = MeshLib::MeshGenerator::generateRegularHexMesh(length, mesh_size);

	ASSERT_TRUE(mesh != nullptr);
	const std::size_t size(mesh_size*mesh_size*mesh_size);

	std::string const prop_name("TestProperty");
	boost::optional<MeshLib::PropertyVector<int> &> properties(
		mesh->getProperties().createNewPropertyVector<int>(prop_name,
			MeshLib::MeshItemType::Cell));
	(*properties).resize(size);
	std::iota((*properties).begin(), (*properties).end(), 1);

	vtkNew<InSituLib::VtkMappedPropertyVectorTemplate<int> > dataArray;
	dataArray->SetPropertyVector(*properties);

	ASSERT_EQ(dataArray->GetNumberOfComponents(), 1);
	ASSERT_EQ(dataArray->GetNumberOfTuples(), size);

	ASSERT_EQ(dataArray->GetValueReference(0), 1);
	double* range = dataArray->GetRange(0);
	ASSERT_EQ(range[0], 1);
	ASSERT_EQ(range[1], 1 + mesh->getNElements() - 1);

	delete mesh;
}

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

std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition(
    BaseLib::ConfigTree const& config,
    MeshLib::Mesh const& mesh,
    int const n_components)
{
	auto field_name = config.getConfParam<std::string>("field_name");
	DBUG("Using field_name %s", field_name.c_str());

	if (!mesh.getProperties().hasPropertyVector(field_name))
	{
		ERR("The required property %s does not exists in the mesh.",
		    field_name.c_str());
		std::abort();
	}
	auto const& property =
	    mesh.getProperties().template getPropertyVector<double>(field_name);
	if (!property)
	{
		ERR("The required property %s is not of the requested type.",
		    field_name.c_str());
		std::abort();
	}

	if (property->getNumberOfComponents() !=
	    static_cast<std::size_t>(n_components))
	{
		ERR("The required property %s has different number of components %d, "
		    "expected %d.",
		    field_name.c_str(), property->getNumberOfComponents(), n_components);
		std::abort();
	}
	return std::unique_ptr<InitialCondition>(
	    new MeshPropertyInitialCondition(*property));
}

std::unique_ptr<ParameterBase> createMeshPropertyParameter(
    BaseLib::ConfigTree const config, MeshLib::Mesh const& mesh)
{
	auto field_name = config.get_optional<std::string>("field_name");
	if (!field_name)
	{
		ERR("Could not find required parameter field_name.");
		std::abort();
	}
	DBUG("Using field_name %s", field_name->c_str());

	if (!mesh.getProperties().hasPropertyVector(*field_name))
	{
		ERR("The required property %s does not exists in the mesh.",
		    field_name->c_str());
		std::abort();
	}
	auto const& property =
	    mesh.getProperties().template getPropertyVector<double>(*field_name);
	if (!property)
	{
		ERR("The required property %s is not of the requested type.",
		    field_name->c_str());
		std::abort();
	}

	return std::unique_ptr<ParameterBase>(
	    new MeshPropertyParameter<double>(*property));
}

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::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
    BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
    NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
    int const component_id,
    const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters)
{
    DBUG("Constructing DirichletBoundaryCondition from config.");
    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
    config.checkConfigParameter("type", "Dirichlet");

    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Dirichlet__parameter}
    auto const param_name = config.getConfigParameter<std::string>("parameter");
    DBUG("Using parameter %s", param_name.c_str());

    auto& param = findParameter<double>(param_name, parameters, 1);

    // In case of partitioned mesh the boundary could be empty, i.e. there is no
    // boundary condition.
#ifdef USE_PETSC
    // This can be extracted to createBoundaryCondition() but then the config
    // parameters are not read and will cause an error.
    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
    // subtree and move the code up in createBoundaryCondition().
    if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
        bc_mesh.getNumberOfElements() == 0)
    {
        return nullptr;
    }
#endif  // USE_PETSC

    return std::make_unique<DirichletBoundaryCondition>(
        param, bc_mesh, dof_table_bulk, variable_id, component_id);
}

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