From 84e1ccb0b4d2d1db824405cce5fabf4117d42d27 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Thu, 3 Aug 2017 14:31:19 +0200
Subject: [PATCH] added mapping to from boundary to bulk
---
.../map-boundary-mesh-to-bulk-mesh/README.md | 18 ++++
.../map-boundary-to-bulk.py | 87 +++++++++++++++++++
2 files changed, 105 insertions(+)
create mode 100644 pre/BCs/map-boundary-mesh-to-bulk-mesh/README.md
create mode 100755 pre/BCs/map-boundary-mesh-to-bulk-mesh/map-boundary-to-bulk.py
diff --git a/pre/BCs/map-boundary-mesh-to-bulk-mesh/README.md b/pre/BCs/map-boundary-mesh-to-bulk-mesh/README.md
new file mode 100644
index 0000000..c793ddf
--- /dev/null
+++ b/pre/BCs/map-boundary-mesh-to-bulk-mesh/README.md
@@ -0,0 +1,18 @@
+# Map the nodes of a boundary mesh to the corresponding nodes of the given bulk mesh
+
+## Requirements
+
+* Both meshes have to be conforming, i.e., the coordinates of the boundary mesh
+ must exactly coincide with the coordinates of some bulk mesh nodes
+* Put more explicitly: The coordinates in bulk and boundary must exactly match,
+ i.e., bit-wise, because this script does not apply any tolerances.
+
+
+## What is output?
+
+* A nodal field `bulk_mesh_node_ids` is added to the given `--boundary` mesh.
+ This nodal field maps the nodes of the boundary mesh to the corresponding
+ nodes of the `--bulk` mesh.
+* The generated `--output` mesh can be used with OGS to specify (inhomogeneous)
+ Boundary conditions along arbitrary surfaces.
+
diff --git a/pre/BCs/map-boundary-mesh-to-bulk-mesh/map-boundary-to-bulk.py b/pre/BCs/map-boundary-mesh-to-bulk-mesh/map-boundary-to-bulk.py
new file mode 100755
index 0000000..cb3bbfd
--- /dev/null
+++ b/pre/BCs/map-boundary-mesh-to-bulk-mesh/map-boundary-to-bulk.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env pvpython
+
+import inspect
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("--bulk", type=str, required=True)
+parser.add_argument("--boundary", type=str, required=True)
+parser.add_argument("--output", type=str, required=True)
+args = parser.parse_args()
+
+#### import the simple module from the paraview
+from paraview.simple import *
+#### disable automatic camera reset on 'Show'
+paraview.simple._DisableFirstRenderCameraReset()
+
+# create a new 'XML Unstructured Grid Reader'
+reader_bulk = XMLUnstructuredGridReader(FileName=[args.bulk])
+reader_boundary = XMLUnstructuredGridReader(FileName=[args.boundary])
+
+
+def do_mapping():
+ bulk = self.GetInputDataObject(0, 1)
+ boundary = self.GetInputDataObject(0, 0)
+
+ points_bulk = bulk.GetPoints()
+ points_boundary = boundary.GetPoints()
+
+ N_boundary = points_boundary.GetNumberOfPoints()
+ N_bulk = points_bulk.GetNumberOfPoints()
+
+ # point_ids = vtk.vtkIdTypeArray() ## TODO why does that not work?
+ # point_ids = vtk.vtkLongArray()
+ point_ids = vtk.vtkUnsignedLongArray()
+ point_ids.SetName("bulk_mesh_node_ids")
+ point_ids.SetNumberOfComponents(1)
+ point_ids.SetNumberOfTuples(N_boundary)
+
+ i_bulk = 0
+ map_coords_to_bulk_point_ids = {}
+
+ for i_boundary in range(N_boundary):
+ coords_boundary = points_boundary.GetPoint(i_boundary)
+
+ try:
+ point_id_bulk = map_coords_to_bulk_point_ids[coords_boundary]
+ except KeyError:
+ # If point has not been mapped yet, search for it
+ # Note: This approach is not very efficient, and in the worst case
+ # it copies all bulk mesh nodes.
+ for i_bulk in range(i_bulk, N_bulk):
+ coords_bulk = points_bulk.GetPoint(i_bulk)
+
+ if coords_bulk == coords_boundary:
+ point_id_bulk = i_bulk
+ # Breaking here makes sure that duplicate nodes are not
+ # found and an error message will be generated
+ break
+
+ map_coords_to_bulk_point_ids[coords_bulk] = i_bulk
+ else:
+ raise ValueError("The boundary coordinates " + str(coords_boundary)
+ + " have not been found in the bulk mesh. Reasons might be"
+ + " (i) those coordinates do not match any in the bulk mesh, or"
+ + " (ii) those coordinates are duplicate in the boundary mesh")
+
+ point_ids.SetComponent(i_boundary, 0, point_id_bulk)
+
+
+ self.GetOutputDataObject(0).GetPointData().AddArray(point_ids)
+
+ # TODO map cells
+
+
+boundary_bulk_mapping = ProgrammableFilter(Input=[reader_boundary, reader_bulk])
+boundary_bulk_mapping.Script = inspect.getsource(do_mapping) + "\n\ndo_mapping()"
+boundary_bulk_mapping.RequestInformationScript = ''
+boundary_bulk_mapping.RequestUpdateExtentScript = ''
+boundary_bulk_mapping.PythonPath = ''
+boundary_bulk_mapping.CopyArrays = 1
+
+
+# save data
+SaveData(args.output, proxy=boundary_bulk_mapping, DataMode='Binary',
+ EncodeAppendedData=1,
+ CompressorType='ZLib')
+
--
GitLab