Detailed Description

A tool for post-processing results from the LIE approach

The tool creates a new mesh containing duplicated fracture nodes to represent geometric discontinuities in visualization.

Definition at line 24 of file PostUtils.h.

#include <PostUtils.h>

Collaboration diagram for ProcessLib::LIE::PostProcessTool:

[legend]

Public Member Functions
	PostProcessTool (MeshLib::Mesh const &org_mesh, std::vector< int > const &vec_fracture_mat_IDs, std::vector< std::vector< MeshLib::Node * > > const &vec_vec_fracture_nodes, std::vector< std::vector< MeshLib::Element * > > const &vec_vec_fracture_matrix_elements, std::vector< std::pair< std::size_t, std::vector< int > > > const &vec_branch_nodeID_matIDs, std::vector< std::pair< std::size_t, std::vector< int > > > const &vec_junction_nodeID_matIDs)
MeshLib::Mesh const &	getOutputMesh () const

Private Member Functions
template<typename T>
MeshLib::PropertyVector< T > *	createProperty (MeshLib::PropertyVector< T > const &property)
template<typename T>
void	copyPropertyValues (MeshLib::PropertyVector< T > const &source_property, MeshLib::PropertyVector< T > *const destination_property)
void	calculateTotalDisplacement (unsigned const n_fractures, unsigned const n_junctions)

Private Attributes
MeshLib::Mesh const &	_org_mesh
std::unique_ptr< MeshLib::Mesh >	_output_mesh
std::map< std::size_t, std::vector< std::size_t > >	_map_dup_newNodeIDs

Constructor & Destructor Documentation

◆ PostProcessTool()

ProcessLib::LIE::PostProcessTool::PostProcessTool	(	MeshLib::Mesh const &	org_mesh,
		std::vector< int > const &	vec_fracture_mat_IDs,
		std::vector< std::vector< MeshLib::Node * > > const &	vec_vec_fracture_nodes,
		std::vector< std::vector< MeshLib::Element * > > const &	vec_vec_fracture_matrix_elements,
		std::vector< std::pair< std::size_t, std::vector< int > > > const &	vec_branch_nodeID_matIDs,
		std::vector< std::pair< std::size_t, std::vector< int > > > const &	vec_junction_nodeID_matIDs )

Definition at line 40 of file PostUtils.cpp.

    : _org_mesh(org_mesh)
{
    if (!org_mesh.getProperties().hasPropertyVector("displacement") ||
        !org_mesh.getProperties().hasPropertyVector("displacement_jump1") ||
        !org_mesh.getProperties().hasPropertyVector("levelset1"))
    {
        OGS_FATAL("The given mesh does not have relevant properties");
    }
 
    // clone nodes and elements
    std::vector<MeshLib::Node*> new_nodes(
        MeshLib::copyNodeVector(org_mesh.getNodes()));
    std::vector<MeshLib::Element*> new_eles(
        MeshLib::copyElementVector(org_mesh.getElements(), new_nodes));
 
    // duplicate fracture nodes (two dup. nodes created at a branch)
    for (auto const& vec_fracture_nodes : vec_vec_fracture_nodes)
    {
        for (auto const* org_node : vec_fracture_nodes)
        {
            auto duplicated_node =
                new MeshLib::Node(org_node->data(), new_nodes.size());
            new_nodes.push_back(duplicated_node);
            _map_dup_newNodeIDs[org_node->getID()].push_back(
                duplicated_node->getID());
        }
    }
    // at a junction, generate one more duplicated node (total 4 nodes)
    for (auto const& entry : vec_junction_nodeID_matIDs)
    {
        auto* org_node = org_mesh.getNode(entry.first);
        auto duplicated_node =
            new MeshLib::Node(org_node->data(), new_nodes.size());
        new_nodes.push_back(duplicated_node);
        _map_dup_newNodeIDs[org_node->getID()].push_back(
            duplicated_node->getID());
    }
 
    // split elements using the new duplicated nodes
    for (unsigned fracture_id = 0;
         fracture_id < vec_vec_fracture_matrix_elements.size();
         fracture_id++)
    {
        auto const frac_matid = vec_fracture_mat_IDs[fracture_id];
        auto const& vec_fracture_matrix_elements =
            vec_vec_fracture_matrix_elements[fracture_id];
        auto const& vec_fracture_nodes = vec_vec_fracture_nodes[fracture_id];
        auto prop_levelset = org_mesh.getProperties().getPropertyVector<double>(
            "levelset" + std::to_string(fracture_id + 1));
        for (auto const* org_e : vec_fracture_matrix_elements)
        {
            // only matrix elements
            if (org_e->getDimension() != org_mesh.getDimension())
            {
                continue;
            }
 
            auto const eid = org_e->getID();
            // keep original if the element has levelset<=0
            if ((*prop_levelset)[eid] <= 0)
            {
                continue;
            }
 
            // replace fracture nodes with duplicated ones
            MeshLib::Element* e = new_eles[eid];
            for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
            {
                const auto node_id = getNodeIndex(*e, i);
                if (!includesNodeID(vec_fracture_nodes, node_id))
                {
                    continue;
                }
 
                // list of duplicated node IDs
                auto itr = _map_dup_newNodeIDs.find(node_id);
                if (itr == _map_dup_newNodeIDs.end())
                {
                    continue;
                }
                const auto& dup_newNodeIDs = itr->second;
 
                // choose new node id
                unsigned new_node_id = 0;
                if (dup_newNodeIDs.size() == 1)
                {
                    // non-intersected nodes
                    new_node_id = dup_newNodeIDs[0];
                }
                else if (dup_newNodeIDs.size() == 2)
                {
                    // branch nodes
                    const auto& br_matids = getMaterialIdsForNode(
                        vec_branch_nodeID_matIDs, node_id);
                    auto const frac2_matid = BaseLib::findFirstNotEqualElement(
                        br_matids, frac_matid);
                    assert(frac2_matid);
                    auto const frac2_id =
                        BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
                    assert(frac2_id != std::numeric_limits<std::size_t>::max());
                    auto prop_levelset2 =
                        org_mesh.getProperties().getPropertyVector<double>(
                            "levelset" + std::to_string(frac2_id + 1));
                    std::size_t pos = 0;
                    if ((*prop_levelset2)[eid] <= 0)
                    {
                        // index of this fracture
                        pos = BaseLib::findIndex(br_matids, frac_matid);
                    }
                    else
                    {
                        // index of the other fracture
                        pos = BaseLib::findIndex(br_matids, *frac2_matid);
                    }
                    assert(pos != std::numeric_limits<std::size_t>::max());
                    new_node_id = dup_newNodeIDs[pos];
                }
                else
                {
                    // junction nodes
                    const auto& jct_matids = getMaterialIdsForNode(
                        vec_junction_nodeID_matIDs, node_id);
                    auto const frac2_matid = BaseLib::findFirstNotEqualElement(
                        jct_matids, frac_matid);
                    assert(frac2_matid);
                    auto const frac2_id =
                        BaseLib::findIndex(vec_fracture_mat_IDs, *frac2_matid);
                    assert(frac2_id != std::numeric_limits<std::size_t>::max());
                    auto prop_levelset2 =
                        org_mesh.getProperties().getPropertyVector<double>(
                            "levelset" + std::to_string(frac2_id + 1));
 
                    //
                    if ((*prop_levelset2)[eid] <= 0)
                    {
                        // index of this frac
                        auto const pos =
                            BaseLib::findIndex(jct_matids, frac_matid);
                        assert(pos != std::numeric_limits<std::size_t>::max());
                        new_node_id = dup_newNodeIDs[pos];
                    }
                    else
                    {
                        // set the last duplicated node
                        new_node_id = dup_newNodeIDs.back();
                    }
                }
 
                // replace node
                e->setNode(i, new_nodes[new_node_id]);
            }
        }
    }
 
    // new mesh
    _output_mesh = std::make_unique<MeshLib::Mesh>(org_mesh.getName(),
                                                   new_nodes, new_eles);
 
    for (auto [name, property] : _org_mesh.getProperties())
    {
        if (auto const* p =
                dynamic_cast<MeshLib::PropertyVector<double>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<float>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<int>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<unsigned>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<long>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<long long>*>(
                         property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<unsigned long>*>(
                         property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<unsigned long long>*>(
                         property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<std::size_t>*>(
                         property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else if (auto const* p =
                     dynamic_cast<MeshLib::PropertyVector<char>*>(property))
        {
            copyPropertyValues(*p, createProperty(*p));
        }
        else
        {
            OGS_FATAL(
                "Mesh property '{:s}' of unhandled data type '{:s}'. Please "
                "check the data type of the mesh properties. The available "
                "data types are:"
                "\n\t double,"
                "\n\t float,"
                "\n\t int,"
                "\n\t unsigned,"
                "\n\t long,"
                "\n\t long long,"
                "\n\t unsigned long,"
                "\n\t unsigned long long,"
                "\n\t char.",
                property->getPropertyName(),
                BaseLib::typeToString(*property));
        }
    }
    calculateTotalDisplacement(vec_vec_fracture_nodes.size(),
                               vec_junction_nodeID_matIDs.size());
}

References _map_dup_newNodeIDs, _org_mesh, _output_mesh, calculateTotalDisplacement(), MeshLib::copyElementVector(), MeshLib::copyNodeVector(), copyPropertyValues(), createProperty(), BaseLib::findFirstNotEqualElement(), BaseLib::findIndex(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getName(), MeshLib::Mesh::getNode(), MeshLib::Mesh::getNodes(), MeshLib::Element::getNumberOfNodes(), MeshLib::Mesh::getProperties(), MeshLib::Properties::getPropertyVector(), MeshLib::Properties::hasPropertyVector(), MeshLib::Node, OGS_FATAL, MeshLib::Element::setNode(), and BaseLib::typeToString().

Member Function Documentation

◆ calculateTotalDisplacement()

void ProcessLib::LIE::PostProcessTool::calculateTotalDisplacement	(	unsigned const	n_fractures,
		unsigned const	n_junctions )

private

Definition at line 377 of file PostUtils.cpp.

{
    auto const& u = *_output_mesh->getProperties().getPropertyVector<double>(
        "displacement");
    auto const n_u_comp = u.getNumberOfGlobalComponents();
    assert(u.size() == _output_mesh->getNodes().size() * 3);
    auto& total_u =
        *_output_mesh->getProperties().createNewPropertyVector<double>(
            "u", MeshLib::MeshItemType::Node, n_u_comp);
    total_u.resize(u.size());
    ranges::copy(u, total_u.begin());
 
    for (unsigned enrich_id = 0; enrich_id < n_fractures + n_junctions;
         enrich_id++)
    {
        // nodal value of levelset
        std::vector<double> nodal_levelset(
            _output_mesh->getNodes().size(),
            std::numeric_limits<double>::quiet_NaN());
 
        auto const& ele_levelset =
            *_output_mesh->getProperties().getPropertyVector<double>(
                "levelset" + std::to_string(enrich_id + 1));
        for (MeshLib::Element const* e : _output_mesh->getElements())
        {
            if (e->getDimension() != _output_mesh->getDimension())
            {
                continue;
            }
            const double e_levelset = ele_levelset[e->getID()];
 
            for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
            {
                nodal_levelset[getNodeIndex(*e, i)] = e_levelset;
            }
        }
 
        // update total displacements
        auto const& g =
            *_output_mesh->getProperties().getPropertyVector<double>(
                "displacement_jump" + std::to_string(enrich_id + 1));
        for (unsigned i = 0; i < _output_mesh->getNodes().size(); i++)
        {
            for (int j = 0; j < n_u_comp; j++)
            {
                total_u[i * n_u_comp + j] +=
                    nodal_levelset[i] * g[i * n_u_comp + j];
            }
        }
    }
}

References _output_mesh, and MeshLib::Node.

Referenced by PostProcessTool().

◆ copyPropertyValues()

template<typename T>

void ProcessLib::LIE::PostProcessTool::copyPropertyValues	(	MeshLib::PropertyVector< T > const &	source_property,
		MeshLib::PropertyVector< T > *const	destination_property )

private

Definition at line 324 of file PostUtils.cpp.

{
    if (destination_property == nullptr)
    {
        // skip the copy, because the destination wasn't created.
        return;
    }
 
    auto const item_type = source_property.getMeshItemType();
    if (item_type == MeshLib::MeshItemType::Node)
    {
        auto const n_src_comp = source_property.getNumberOfGlobalComponents();
        auto const n_dest_comp =
            destination_property->getNumberOfGlobalComponents();
        // copy existing
        for (unsigned i = 0; i < _org_mesh.getNumberOfNodes(); i++)
        {
            auto last = std::copy_n(&source_property[i * n_src_comp],
                                    n_src_comp,
                                    &(*destination_property)[i * n_dest_comp]);
            // set zero for components not existing in the original
            std::fill_n(last, n_dest_comp - n_src_comp, T{0});
        }
        // copy duplicated
        for (auto itr : _map_dup_newNodeIDs)
        {
            for (unsigned k = 0; k < itr.second.size(); k++)
            {
                std::copy_n(
                    &(*destination_property)[itr.first * n_dest_comp],
                    n_dest_comp,
                    &(*destination_property)[itr.second[k] * n_dest_comp]);
            }
        }
    }
    else if (item_type == MeshLib::MeshItemType::Cell)
    {
        assert(source_property.size() == destination_property->size());
        std::copy(source_property.begin(), source_property.end(),
                  destination_property->begin());
    }
    else
    {
        OGS_FATAL(
            "Property '{:s}' values cannot be copied because its mesh item "
            "type '{:s}' is not supported. Unexpected error, because the "
            "destination property was created.",
            source_property.getPropertyName(), toString(item_type));
    }
}

References _map_dup_newNodeIDs, _org_mesh, MeshLib::PropertyVector< PROP_VAL_TYPE >::begin(), MeshLib::Cell, MeshLib::PropertyVector< PROP_VAL_TYPE >::end(), MeshLib::PropertyVectorBase::getMeshItemType(), MeshLib::PropertyVectorBase::getNumberOfGlobalComponents(), MeshLib::PropertyVectorBase::getPropertyName(), MeshLib::Node, OGS_FATAL, and MeshLib::PropertyVector< PROP_VAL_TYPE >::size().

Referenced by PostProcessTool().

◆ createProperty()

template<typename T>

MeshLib::PropertyVector< T > * ProcessLib::LIE::PostProcessTool::createProperty ( MeshLib::PropertyVector< T > const & property )

private

Definition at line 288 of file PostUtils.cpp.

{
    auto const item_type = property.getMeshItemType();
    auto const n_src_comp = property.getNumberOfGlobalComponents();
    // convert 2D vector to 3D. Otherwise Paraview Calculator filter does
    // not recognize it as a vector
    auto const n_dest_comp = (n_src_comp == 2) ? 3 : n_src_comp;
 
    auto new_property = MeshLib::getOrCreateMeshProperty<T>(
        *_output_mesh, property.getPropertyName(), item_type, n_dest_comp);
 
    if (item_type == MeshLib::MeshItemType::Node)
    {
        assert(new_property->size() ==
               _output_mesh->getNumberOfNodes() * n_dest_comp);
    }
    else if (item_type == MeshLib::MeshItemType::Cell)
    {
        assert(new_property->size() ==
               _output_mesh->getNumberOfElements() * n_dest_comp);
    }
    else
    {
        WARN(
            "Property '{:s}' cannot be created because its mesh item type "
            "'{:s}' is not supported.",
            property.getPropertyName(), toString(item_type));
        _output_mesh->getProperties().removePropertyVector(
            std::string(new_property->getPropertyName()));
        return nullptr;
    }
    return new_property;
}