[PCS/LIE] add post-processing tool

ProcessLib/SmallDeformationWithLIE/Common/PostUtils.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "PostUtils.h"
+
+#include <map>
+#include <vector>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
+#include "MeshLib/MeshEditing/DuplicateMeshComponents.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+namespace
+{
+
+template <typename T>
+inline void sort_unique(std::vector<T> &vec)
+{
+    std::sort(vec.begin(), vec.end());
+    vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
+}
+
+} // no named namespace
+
+PostProcessTool::PostProcessTool(
+    MeshLib::Mesh const& org_mesh,
+    std::vector<MeshLib::Node*> const& vec_fracture_nodes,
+    std::vector<MeshLib::Element*> const& vec_fracutre_matrix_elements)
+    :_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
+    for (auto const* org_node : vec_fracture_nodes)
+    {
+        auto duplicated_node = new MeshLib::Node(org_node->getCoords(), new_nodes.size());
+        new_nodes.push_back(duplicated_node);
+        _map_dup_newNodeIDs[org_node->getID()] = duplicated_node->getID();
+    }
+
+    // split elements using the new duplicated nodes
+    auto prop_levelset = org_mesh.getProperties().getPropertyVector<double>("levelset1");
+    for (auto const* org_e : vec_fracutre_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++)
+        {
+            auto itr = _map_dup_newNodeIDs.find(e->getNodeIndex(i));
+            if (itr == _map_dup_newNodeIDs.end())
+                continue;
+
+            auto itr2 = std::find_if(vec_fracture_nodes.begin(), vec_fracture_nodes.end(),
+                                     [&](MeshLib::Node const*node) { return node->getID()==e->getNodeIndex(i);});
+            if (itr2 == vec_fracture_nodes.end())
+                continue;
+            e->setNode(i, new_nodes[itr->second]);
+        }
+    }
+
+    // new mesh
+    _output_mesh.reset(new MeshLib::Mesh(org_mesh.getName(), new_nodes, new_eles));
+    createProperties<int>();
+    createProperties<double>();
+    copyProperties<int>();
+    copyProperties<double>();
+    calculateTotalDisplacement();
+}
+
+template <typename T>
+void PostProcessTool::createProperties()
+{
+    MeshLib::Properties const& src_properteis = _org_mesh.getProperties();
+    for (auto name : src_properteis.getPropertyVectorNames())
+    {
+        auto const* src_prop = src_properteis.getPropertyVector<T>(name);
+        if (!src_prop)
+            continue;
+        auto const n_src_comp = src_prop->getNumberOfComponents();
+        // 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;
+
+        if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
+        {
+            auto new_prop =
+                _output_mesh->getProperties().createNewPropertyVector<T>(
+                    name, MeshLib::MeshItemType::Node, n_dest_comp);
+            new_prop->resize(_output_mesh->getNumberOfNodes() * n_dest_comp);
+        }
+        else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
+        {
+            auto new_prop =
+                _output_mesh->getProperties().createNewPropertyVector<T>(
+                    name, MeshLib::MeshItemType::Cell, n_dest_comp);
+            new_prop->resize(src_prop->size());
+        }
+    }
+}
+
+template <typename T>
+void PostProcessTool::copyProperties()
+{
+    MeshLib::Properties const& src_properteis = _org_mesh.getProperties();
+    for (auto name : src_properteis.getPropertyVectorNames())
+    {
+        auto const* src_prop = src_properteis.getPropertyVector<T>(name);
+        if (!src_prop)
+            continue;
+        auto* dest_prop = _output_mesh->getProperties().getPropertyVector<T>(name);
+        assert(dest_prop);
+
+        if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Node)
+        {
+            auto const n_src_comp = src_prop->getNumberOfComponents();
+            auto const n_dest_comp = dest_prop->getNumberOfComponents();
+            // copy existing
+            for (unsigned i=0; i<_org_mesh.getNumberOfNodes(); i++)
+            {
+                for (unsigned j=0; j<n_src_comp; j++)
+                    (*dest_prop)[i*n_dest_comp+j] = (*src_prop)[i*n_src_comp+j];
+                // set zero for components not existing in the original
+                for (unsigned j=n_src_comp; j<n_dest_comp; j++)
+                    (*dest_prop)[i*n_dest_comp+j] = 0;
+            }
+            // copy duplicated
+            for (auto itr : _map_dup_newNodeIDs)
+            {
+                for (unsigned j=0; j<n_dest_comp; j++)
+                    (*dest_prop)[itr.second*n_dest_comp + j] = (*dest_prop)[itr.first*n_dest_comp + j];
+            }
+        }
+        else if (src_prop->getMeshItemType() == MeshLib::MeshItemType::Cell)
+        {
+            std::copy(src_prop->begin(), src_prop->end(), dest_prop->begin());
+        }
+    }
+}
+
+void PostProcessTool::calculateTotalDisplacement()
+{
+    // nodal value of levelset
+    std::vector<double> nodal_levelset(_output_mesh->getNodes().size(), 0.0);
+    auto const& ele_levelset =  *_output_mesh->getProperties().getPropertyVector<double>("levelset1");
+    for (MeshLib::Element const* e : _output_mesh->getElements())
+    {
+        if (e->getDimension() != _output_mesh->getDimension())
+            continue;
+        const double e_levelset = ele_levelset[e->getID()];
+        if (e_levelset == 0)
+            continue;
+
+        for (unsigned i=0; i<e->getNumberOfNodes(); i++)
+            nodal_levelset[e->getNodeIndex(i)] = e_levelset;
+    }
+
+    // total displacements
+    auto const& u =  *_output_mesh->getProperties().getPropertyVector<double>("displacement");
+    auto const n_u_comp = u.getNumberOfComponents();
+    assert(u.size() == _output_mesh->getNodes().size() * 3);
+    auto const& g =  *_output_mesh->getProperties().getPropertyVector<double>("displacement_jump1");
+    auto& total_u =
+        *_output_mesh->getProperties().createNewPropertyVector<double>(
+            "u", MeshLib::MeshItemType::Node, n_u_comp);
+    total_u.resize(u.size());
+    for (unsigned i=0; i<_output_mesh->getNodes().size(); i++)
+    {
+        for (unsigned j=0; j<n_u_comp; j++)
+            total_u[i*n_u_comp+j] = u[i*n_u_comp+j] + nodal_levelset[i] * g[i*n_u_comp+j];
+    }
+}
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib

ProcessLib/SmallDeformationWithLIE/Common/PostUtils.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_POSTUTILS_H_
+#define PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_POSTUTILS_H_
+
+#include <map>
+#include <memory>
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationWithLIE
+{
+
+/// 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.
+class PostProcessTool
+{
+public:
+    PostProcessTool(
+        MeshLib::Mesh const& org_mesh,
+        std::vector<MeshLib::Node*> const& vec_fracture_nodes,
+        std::vector<MeshLib::Element*> const& vec_fracutre_matrix_elements);
+
+    MeshLib::Mesh const& getOutputMesh() const { return *_output_mesh; }
+
+private:
+    template <typename T>
+    void createProperties();
+    template <typename T>
+    void copyProperties();
+    void calculateTotalDisplacement();
+
+    MeshLib::Mesh const& _org_mesh;
+    std::unique_ptr<MeshLib::Mesh> _output_mesh;
+    std::map<std::size_t, std::size_t> _map_dup_newNodeIDs;
+};
+
+}  // namespace SmallDeformationWithLIE
+}  // namespace ProcessLib
+
+#endif // PROCESSLIB_SMALLDEFORMATION_WITH_LIE_COMMON_POSTUTILS_H_