diff --git a/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.cpp b/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7f043800a4bd5ce994616c41f6df39563b84615a --- /dev/null +++ b/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.cpp @@ -0,0 +1,205 @@ +/** + * \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 diff --git a/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.h b/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.h new file mode 100644 index 0000000000000000000000000000000000000000..60c2f4edc14507d6d918af691d882a8f9ec2a7b2 --- /dev/null +++ b/ProcessLib/SmallDeformationWithLIE/Common/PostUtils.h @@ -0,0 +1,52 @@ +/** + * \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_