diff --git a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
index 5486716bc37fca870cac6a98556749910456e8d2..4e0d1870644d98361eabf1da7e15330e307613df 100644
--- a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
@@ -21,114 +21,12 @@
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/IO/writeMeshToFile.h"
-#include "GeoLib/AnalyticalGeometry.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/Polygon.h"
-#include "MathLib/Vector3.h"
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+#include "MeshGeoToolsLib/MeshEditing/ResetMeshElementProperty.h"
#include "MeshLib/Mesh.h"
-#include "MeshLib/Node.h"
-#include "MeshLib/Elements/Element.h"
-
-static std::vector<bool> markNodesOutSideOfPolygon(
- std::vector<MeshLib::Node*> const& nodes, GeoLib::Polygon const& polygon)
-{
- // *** rotate polygon to xy_plane
- MathLib::Vector3 normal;
- GeoLib::Polygon rot_polygon(GeoLib::rotatePolygonToXY(polygon, normal));
-
- // *** rotate mesh nodes to xy-plane
- // 1 copy all mesh nodes to GeoLib::Points
- std::vector<GeoLib::Point*> rotated_nodes;
- for (auto node : nodes)
- rotated_nodes.push_back(new GeoLib::Point(*node, node->getID()));
- // 2 rotate the Points
- MathLib::DenseMatrix<double> rot_mat(3,3);
- GeoLib::computeRotationMatrixToXY(normal, rot_mat);
- GeoLib::rotatePoints(rot_mat, rotated_nodes);
- // 3 set z coord to zero
- std::for_each(rotated_nodes.begin(), rotated_nodes.end(),
- [] (GeoLib::Point* p) { (*p)[2] = 0.0; }
- );
-
- // *** mark rotated nodes
- std::vector<bool> outside(rotated_nodes.size(), true);
- for (std::size_t k(0); k<rotated_nodes.size(); k++) {
- if (rot_polygon.isPntInPolygon(*(rotated_nodes[k]))) {
- outside[k] = false;
- }
- }
-
- for (auto & rotated_node : rotated_nodes)
- delete rotated_node;
-
- std::vector<GeoLib::Point*> & rot_polygon_pnts(
- const_cast<std::vector<GeoLib::Point*> &>(
- rot_polygon.getPointsVec()
- )
- );
- for (auto & rot_polygon_pnt : rot_polygon_pnts)
- delete rot_polygon_pnt;
-
- return outside;
-}
-
-template <typename PT>
-void resetMeshElementProperty(MeshLib::Mesh& mesh,
- GeoLib::Polygon const& polygon,
- std::string const& property_name,
- PT new_property_value,
- int restrict_to_material_id)
-{
- auto* const pv = MeshLib::getOrCreateMeshProperty<PT>(
- mesh, property_name, MeshLib::MeshItemType::Cell, 1);
-
- if (pv->getMeshItemType() != MeshLib::MeshItemType::Cell)
- {
- ERR("Values of the PropertyVector are not assigned to cells.");
- return;
- }
-
- std::vector<bool> outside(
- markNodesOutSideOfPolygon(mesh.getNodes(), polygon));
-
- auto const* material_ids =
- mesh.getProperties().getPropertyVector<int>("MaterialIDs");
-
- auto hasElementRequiredMaterialID = [&](int const element_id) {
- if (restrict_to_material_id == -1)
- return true;
- if (!material_ids)
- {
- OGS_FATAL(
- "Restriction of reseting a property in a polygonal region "
- "requires that a MaterialIDs data array is available in the "
- "mesh.");
- }
- return (*material_ids)[element_id] == restrict_to_material_id;
- };
-
- for (std::size_t j(0); j < mesh.getElements().size(); ++j)
- {
- bool elem_out(true);
- MeshLib::Element const* const elem(mesh.getElements()[j]);
- for (auto k = decltype(elem->getNumberOfNodes()){0};
- k < elem->getNumberOfNodes() && elem_out;
- ++k)
- {
- if (!outside[elem->getNode(k)->getID()])
- {
- elem_out = false;
- }
- }
- if (!elem_out && hasElementRequiredMaterialID(elem->getID()))
- {
- (*pv)[j] = new_property_value;
- }
- }
-}
int main (int argc, char* argv[])
{
@@ -220,21 +118,21 @@ int main (int argc, char* argv[])
std::string const& property_name(property_name_arg.getValue());
if (char_property_arg.isSet()) {
- resetMeshElementProperty(*mesh, polygon, property_name,
- char_property_arg.getValue(),
- restrict_arg.getValue());
+ MeshGeoToolsLib::resetMeshElementProperty(*mesh, polygon, property_name,
+ char_property_arg.getValue(),
+ restrict_arg.getValue());
}
if (int_property_arg.isSet()) {
- resetMeshElementProperty(*mesh, polygon, property_name,
- int_property_arg.getValue(),
- restrict_arg.getValue());
+ MeshGeoToolsLib::resetMeshElementProperty(*mesh, polygon, property_name,
+ int_property_arg.getValue(),
+ restrict_arg.getValue());
}
if (bool_property_arg.isSet()) {
- resetMeshElementProperty(*mesh, polygon, property_name,
- bool_property_arg.getValue(),
- restrict_arg.getValue());
+ MeshGeoToolsLib::resetMeshElementProperty(*mesh, polygon, property_name,
+ bool_property_arg.getValue(),
+ restrict_arg.getValue());
}
std::vector<std::string> property_names(
diff --git a/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
new file mode 100644
index 0000000000000000000000000000000000000000..219a42b5fb9e76d6c4ed949b9411cc266fdd2389
--- /dev/null
+++ b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
@@ -0,0 +1,68 @@
+/*
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <vector>
+
+#include "GeoLib/AnalyticalGeometry.h"
+#include "GeoLib/Polygon.h"
+
+#include "MathLib/Vector3.h"
+#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+
+#include "MeshLib/Node.h"
+
+namespace MeshGeoToolsLib
+{
+std::vector<bool> markNodesOutSideOfPolygon(
+ std::vector<MeshLib::Node*> const& nodes, GeoLib::Polygon const& polygon)
+{
+ // *** rotate polygon to xy_plane
+ MathLib::Vector3 normal;
+ GeoLib::Polygon rot_polygon(GeoLib::rotatePolygonToXY(polygon, normal));
+
+ // *** rotate mesh nodes to xy-plane
+ // 1 copy all mesh nodes to GeoLib::Points
+ std::vector<GeoLib::Point*> rotated_nodes;
+ for (auto node : nodes)
+ rotated_nodes.push_back(new GeoLib::Point(*node, node->getID()));
+ // 2 rotate the Points
+ MathLib::DenseMatrix<double> rot_mat(3,3);
+ GeoLib::computeRotationMatrixToXY(normal, rot_mat);
+ GeoLib::rotatePoints(rot_mat, rotated_nodes);
+ // 3 set z coord to zero
+ std::for_each(rotated_nodes.begin(), rotated_nodes.end(),
+ [] (GeoLib::Point* p) { (*p)[2] = 0.0; }
+ );
+
+ // *** mark rotated nodes
+ std::vector<bool> outside(rotated_nodes.size(), true);
+ for (std::size_t k(0); k<rotated_nodes.size(); k++) {
+ if (rot_polygon.isPntInPolygon(*(rotated_nodes[k]))) {
+ outside[k] = false;
+ }
+ }
+
+ for (auto & rotated_node : rotated_nodes)
+ delete rotated_node;
+
+ std::vector<GeoLib::Point*> & rot_polygon_pnts(
+ const_cast<std::vector<GeoLib::Point*> &>(
+ rot_polygon.getPointsVec()
+ )
+ );
+ for (auto & rot_polygon_pnt : rot_polygon_pnts)
+ delete rot_polygon_pnt;
+
+ return outside;
+}
+
+} // end namespace MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/MeshEditing/ResetMeshElementProperty.h b/MeshGeoToolsLib/MeshEditing/ResetMeshElementProperty.h
new file mode 100644
index 0000000000000000000000000000000000000000..3f9b6e1940ce0f5f134811404c3bf3a49b0db048
--- /dev/null
+++ b/MeshGeoToolsLib/MeshEditing/ResetMeshElementProperty.h
@@ -0,0 +1,87 @@
+/*
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cstdlib>
+#include <vector>
+
+#include <tclap/CmdLine.h>
+
+#include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/readGeometryFromFile.h"
+
+#include "MeshLib/IO/readMeshFromFile.h"
+#include "MeshLib/IO/writeMeshToFile.h"
+
+#include "GeoLib/GEOObjects.h"
+#include "GeoLib/Polygon.h"
+
+#include "MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h"
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+#include "MeshLib/Elements/Element.h"
+
+namespace MeshGeoToolsLib
+{
+template <typename PT>
+void resetMeshElementProperty(MeshLib::Mesh& mesh,
+ GeoLib::Polygon const& polygon,
+ std::string const& property_name,
+ PT new_property_value,
+ int restrict_to_material_id)
+{
+ auto* const pv = MeshLib::getOrCreateMeshProperty<PT>(
+ mesh, property_name, MeshLib::MeshItemType::Cell, 1);
+
+ if (pv->getMeshItemType() != MeshLib::MeshItemType::Cell)
+ {
+ ERR("Values of the PropertyVector are not assigned to cells.");
+ return;
+ }
+
+ auto is_node_outside =
+ [outside = markNodesOutSideOfPolygon(mesh.getNodes(), polygon)](
+ auto const* node_ptr) { return outside[node_ptr->getID()]; };
+
+ auto const* material_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+
+ if (restrict_to_material_id != -1 && !material_ids)
+ {
+ OGS_FATAL(
+ "Restriction of reseting a property in a polygonal region "
+ "requires that a MaterialIDs data array is available in the "
+ "mesh.");
+ }
+
+ auto has_element_required_material_id = [&](int const element_id) {
+ return restrict_to_material_id == -1 ||
+ (*material_ids)[element_id] == restrict_to_material_id;
+ };
+
+ for (std::size_t j(0); j < mesh.getElements().size(); ++j)
+ {
+ MeshLib::Element const* const elem(mesh.getElements()[j]);
+ if (std::all_of(elem->getNodes(),
+ elem->getNodes() + elem->getNumberOfNodes(),
+ is_node_outside))
+ {
+ continue;
+ }
+ if (has_element_required_material_id(elem->getID()))
+ {
+ (*pv)[j] = new_property_value;
+ }
+ }
+}
+
+} // end namespace