Vectors for Replay Mode

.gitignore

@@ -27,3 +27,14 @@ testMesh.txt
 *.synctex.gz
 compile_commands.json
 build/
+contrib/colored*
+contrib/mapped/
+contrib/precice-run/
+contrib/bunny*
+*.vtk
+*.vtu
+*.pvtu
+0*
+1*
+*events.json
+*.log

CMakeLists.txt

@@ -53,3 +53,4 @@ file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/eval_mesh.py           DESTINATION ${C
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/vtk_calculator.py      DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/mesh_io.py             DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/mesh.py                DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/precice_to_aste.py     DESTINATION ${CMAKE_CURRENT_BINARY_DIR})

src/mesh_io.py

@@ -5,34 +5,36 @@
 
 import os
 import logging
+import numpy as np
 
 
-def read_mesh(filename, tag = None):
+def read_mesh(filename, tag = None, datadim=1):
     """
     Returns Mesh Points, Mesh Cells, Mesh Celltypes, 
     Mesh Pointdata in this order
     """
     if os.path.splitext(filename)[1] == ".txt":
         return read_txt(filename)
     else:
-        return read_vtk(filename, tag)
+        return read_vtk(filename, tag, datadim)
 
 
-def write_mesh(filename, points, cells = None, cell_types = None, values = None):
+def write_mesh(filename, points, cells = None, cell_types = None, values = None, datadim=1):
     if os.path.splitext(filename)[1] == ".txt":
         return write_txt(filename, points, cells, values)
     else:
-        return write_vtk(filename, points, cells, cell_types, values)
+        return write_vtk(filename, points, cells, cell_types, values, datadim=datadim)
 
 
-def read_vtk(filename, tag = None):
+def read_vtk(filename, tag = None, datadim=1):
     import vtk
     vtkmesh = read_dataset(filename)
     points = []
     cells = []
     pointdata = []
     cell_types = []
     points = [vtkmesh.GetPoint(i) for i in range(vtkmesh.GetNumberOfPoints())]
+    print("Read " +filename)
     for i in range(vtkmesh.GetNumberOfCells()):
         cell = vtkmesh.GetCell(i)
         cell_type = cell.GetCellType()
@@ -45,12 +47,22 @@ def read_vtk(filename, tag = None):
         cells.append(entry)
     if not tag:
         # vtk Python utility method. Same as tag=="scalars"
-        fieldData = vtkmesh.GetPointData().GetScalars()
+        if datadim==1:
+            fieldData = vtkmesh.GetPointData().GetScalars()
+        if datadim >1:
+            fieldData = vtkmesh.GetPointData().GetVectors()
     else:
         fieldData = vtkmesh.GetPointData().GetAbstractArray(tag)
     if fieldData:
         for i in range(vtkmesh.GetNumberOfPoints()):
-            pointdata.append(fieldData.GetTuple1(i))
+            if datadim == 1:
+                pointdata.append(fieldData.GetTuple1(i))
+            elif datadim == 2:
+                pointdata.append(fieldData.GetTuple2(i))
+            elif datadim == 3:
+                pointdata.append(fieldData.GetTuple3(i))
+            else :
+                raise Exception("Something weird is going on")
     return points, cells, cell_types, pointdata
 
 
@@ -74,6 +86,8 @@ def read_dataset(filename):
     else:
         raise Exception("Unkown File extension: " + extension)
     reader.SetFileName(filename)
+    reader.ReadAllScalarsOn()
+    reader.ReadAllVectorsOn()
     reader.Update()
     return reader.GetOutput()
 
@@ -108,11 +122,14 @@ def read_txt(filename):
         for line in fh:
             point = ()
             parts = line.split(" ")
-            for i in range(3):
-                point += (float(parts[i]),)
+            parts_floats = [float(i) for i in parts]
+            point = parts_floats[0:3]
             points.append(point)
-            if len(parts) > 3:
-                pointdata.append(float(parts[3]))
+            if len(parts_floats)==4:
+                data = parts_floats[3]
+            elif len(parts_floats)>4:
+                data = parts_floats[3:]
+                pointdata.append(data)
     base, ext = os.path.splitext(filename)
     connFileName = base + ".conn" + ext
     if os.path.exists(connFileName):
@@ -121,17 +138,24 @@ def read_txt(filename):
     return points, cells, cell_types, pointdata
 
 
-def write_vtk(filename, points, cells = None, cell_types = None, pointdata = None, tag = None):
+def write_vtk(filename, points, cells = None, cell_types = None, pointdata = None, tag = None, datadim=1):
     import vtk
     data = vtk.vtkUnstructuredGrid() # is also vtkDataSet
-    scalars = vtk.vtkDoubleArray()
+    DataArray = vtk.vtkDoubleArray()
+    DataArray.SetNumberOfComponents(datadim)
+    pointdata = np.array(pointdata)
     if tag:
-        scalars.SetName(tag)
+        DataArray.SetName(tag)
     vtkpoints = vtk.vtkPoints()
     for i, point in enumerate(points):
         vtkpoints.InsertPoint(i, point)
         if pointdata is not None and len(pointdata) > 0:
-            scalars.InsertTuple1(i, pointdata[i])
+            if len(pointdata.shape) == 1 : # scalar case
+                DataArray.InsertTuple1(i, pointdata[i])
+            elif len(pointdata.shape) > 1 and pointdata.shape[1]==2: # two dimensional data
+                DataArray.InsertTuple2(i, pointdata[i,0], pointdata[i,1])
+            elif len(pointdata.shape) > 1 and pointdata.shape[1]==3: # 3D-data
+                DataArray.InsertTuple3(i, pointdata[i,0], pointdata[i,1], pointdata[i,2])
     data.SetPoints(vtkpoints)
     if cells:
         cellArray = vtk.vtkCellArray()
@@ -145,7 +169,13 @@ def write_vtk(filename, points, cells = None, cell_types = None, pointdata = Non
             cellArray.InsertNextCell(vtkCell)
         data.SetCells(cell_types, cellArray)
     pointData = data.GetPointData()
-    pointData.SetScalars(scalars)
+    if len(pointdata.shape) == 1:
+        print(pointdata.shape)
+        pointData.SetScalars(DataArray)
+        print("Writing scalar data...")
+    else:
+        pointData.SetVectors(DataArray)
+        print("Writing Vector data...")
     write_dataset(filename, data)
 
     writer = vtk.vtkUnstructuredGridWriter()
@@ -172,7 +202,13 @@ def write_txt(filename, points, cells = [], pointdata = None):
         for i, point in enumerate(points):
             entry = (str(point[0]), str(point[1]), str(point[2]))
             if pointdata is not None and len(pointdata) > 0:
-                entry += (str(float(pointdata[i])),)
+                if type(pointdata[0]) is float:
+                    entry += (str(float(pointdata[i])),)
+                elif len(pointdata[0])==3:
+                    entry += (str(float(pointdata[i][0])), 
+                              str(float(pointdata[i][1])), 
+                              str(float(pointdata[i][2])))
+
             fh.write(" ".join(entry) + "\n")
 
     base, ext = os.path.splitext(filename)

src/partition_mesh.py

@@ -17,15 +17,15 @@ def main():
     if not algorithm:
         logging.info("No algorithm given. Defaulting to \"meshfree\"")
         algorithm = "meshfree"
-    rootmesh = read_mesh(mesh_names[0], args.tag)
+    rootmesh = read_mesh(mesh_names[0], args.tag, args.datadim)
     if args.numparts > 1:
         part = partition(rootmesh, args.numparts, algorithm)
     else:
         part = [0] *  len(rootmesh.points)
 
     for mesh_name in mesh_names:
         logging.info("Processing mesh " + mesh_name)
-        mesh = read_mesh(mesh_name, args.tag)
+        mesh = read_mesh(mesh_name, args.tag, args.datadim)
         logging.debug("Checking if meshes are matching...")
         assert mesh.points == rootmesh.points, ("Non-matching meshes detected!")
         meshes = apply_partition(mesh, part, args.numparts)
@@ -66,12 +66,11 @@ def __str__(self):
         return "Mesh with {} Points and {} Cells ({} Cell Types)".format(len(self.points), len(self.cells), len(self.cell_types))
 
 
-def read_mesh(filename, tag):
-    points, cells, cell_types, pointdata = mesh_io.read_mesh(filename, tag)
+def read_mesh(filename, tag, datadim=1):
+    points, cells, cell_types, pointdata = mesh_io.read_mesh(filename, tag, datadim=datadim)
     return Mesh(points, cells, cell_types, pointdata)
 
 
-
 def partition(mesh, numparts, algorithm):
     """
     Partitions a mesh using METIS or kmeans. This does not call METIS directly, 
@@ -244,7 +243,6 @@ def apply_partition(orig_mesh, part, numparts):
     """
     meshes = [Mesh() for _ in range(numparts)]
     mapping = {}  # Maps global index to partition and local index
-    print(orig_mesh)
     for i in range(len(orig_mesh.points)):
         partition = part[i]
         selected = meshes[partition]
@@ -299,6 +297,7 @@ def parse_args():
     parser.add_argument("--log", "-l", dest="logging", default="INFO", 
             choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], 
             help="Set the log level. Default is INFO")
+    parser.add_argument("--datadim", "-d", dest="datadim", default=1, type=int, help="Dimensions of the data. Default is 1 (scalar data).")
     return parser.parse_args()
 
 if __name__ == "__main__":

src/preciceMap.cpp

@@ -17,6 +17,7 @@
 INITIALIZE_EASYLOGGINGPP
 
 namespace fs = boost::filesystem;
+using VectorData = std::vector<std::array<double, 3>>;
 
 
 /// Returns list of meshes matching meshname or meshname.dtN for the given rank
@@ -58,8 +59,10 @@ struct Mesh {
   std::vector<std::array<size_t, 2>> edges;
   std::vector<std::array<size_t, 3>> triangles;
   std::vector<double> data;
+  size_t dataStride = 1;
 };
 
+
 using VertexID = int;
 using EdgeID = int;
 
@@ -94,18 +97,25 @@ void readMainFile(Mesh& mesh, const std::string& filename, bool read_data)
     std::ifstream mainFile{filename};
     std::string line;
     while (std::getline(mainFile, line)){
-        double x, y, z, val;
-        std::istringstream iss(line);
-        iss >> x >> y >> z >> val; // split up by whitespace
-        std::array<double, 3> vertexPos{x, y, z};
-        mesh.positions.push_back(vertexPos);
-        if (read_data) // val is ignored on B.
-            mesh.data.push_back(val);
+        std::vector<std::string> parts;
+        boost::split(parts, line, [](char c){return c==' ';});
+        std::vector<double> linevalues(parts.size());
+        std::transform(parts.begin(), parts.end(), linevalues.begin(), [](const std::string& s) -> double{return std::stod(s);} );
+        std::array<double, 3> position {linevalues[0], linevalues[1], linevalues[2]};
+        mesh.positions.push_back(position);       
+        if (read_data){
+          mesh.dataStride = linevalues.size() - 3;
+          for (size_t i = 3; i < linevalues.size(); i++){
+            mesh.data.push_back(linevalues[i]);
+          }
+        }
     }
     if (not read_data)
-        mesh.data = std::vector<double>(mesh.positions.size(), 0);
+        mesh.data = std::vector<double>(mesh.positions.size()*mesh.dataStride, 0);
 }
 
+
+
 // Reads the connectivity file containing the triangle and edge information
 void readConnFile(Mesh& mesh, const std::string& filename)
 {
@@ -166,6 +176,7 @@ int main(int argc, char* argv[])
   const int dataID = interface.getDataID("Data", meshID);
 
   // reads in mesh, 0 data for participant B
+
   auto mesh = readMesh(meshes[0], participant == "A");
   VLOG(1) << "The mesh contains:\n " << mesh.positions.size() << " Vertices\n " << mesh.data.size() << " Data\n " << mesh.edges.size()  << " Edges\n " << mesh.triangles.size() << " Triangles";
 
@@ -237,16 +248,16 @@ int main(int argc, char* argv[])
   edgeMap.clear();
   VLOG(1) << "Mesh setup completed on Rank " << MPIrank;
 
-  interface.initialize();
+  auto dt = interface.initialize();
 
   if (interface.isActionRequired(precice::constants::actionWriteInitialData())) {
     VLOG(1) << "Write initial data for participant " << participant;
-    interface.writeBlockScalarData(dataID, mesh.data.size(), vertexIDs.data(), mesh.data.data());
-    std::ostringstream oss;
-    for(size_t i = 0; i<std::min((size_t)10, mesh.data.size()); ++i)
-        oss << mesh.data[i] << ' ';
-    VLOG(1) << "Data written: " << oss.str();
-
+    if (mesh.dataStride>1){
+      interface.writeBlockVectorData(dataID, mesh.data.size(), vertexIDs.data(), mesh.data.data());
+    }
+    else{
+      interface.writeBlockScalarData(dataID, mesh.data.size(), vertexIDs.data(), mesh.data.data());
+    }
     interface.markActionFulfilled(precice::constants::actionWriteInitialData());
   }
   interface.initializeData();
@@ -255,27 +266,26 @@ int main(int argc, char* argv[])
   while (interface.isCouplingOngoing() and round < meshes.size()) {
     if (participant == "A" and not mesh.data.empty()) {
       auto filename = meshes[round];
-
+      std::ostringstream oss;
       VLOG(1) << "Read mesh for t=" << round << " from " << filename;
       if (not boost::filesystem::exists(filename))
         throw std::runtime_error("File does not exist: " + filename);
+
       auto roundmesh = readMesh(filename, participant == "A");
-      VLOG(1) << "This roundmesh contains:\n " << roundmesh.positions.size() << " Vertices\n " << roundmesh.data.size() << " Data\n " << roundmesh.edges.size()  << " Edges\n " << roundmesh.triangles.size() << " Triangles";
-      assert(roundmesh.data.size() == vertexIDs.size());
-      interface.writeBlockScalarData(dataID, roundmesh.data.size(), vertexIDs.data(), roundmesh.data.data());
-    std::ostringstream oss;
-    for(size_t i = 0; i<std::min((size_t)10, mesh.data.size()); ++i)
-        oss << roundmesh.data[i] << ' ';
-    VLOG(1) << "Data written: " << oss.str();
+      if (mesh.dataStride>1){
+          if (roundmesh.data.size() != roundmesh.positions.size()*roundmesh.dataStride){
+            throw std::runtime_error{std::string{"Size of data does not equal number of vertices by datadimension."}};
+          }
+          interface.writeBlockVectorData(dataID, roundmesh.positions.size(), vertexIDs.data(), roundmesh.data.data());
+      }
+      else{ //scalar data
+          interface.writeBlockScalarData(dataID, roundmesh.data.size(), vertexIDs.data(), roundmesh.data.data());
+      }
     }
-    interface.advance(1);
+    interface.advance(dt);
 
     if (participant == "B") {
       interface.readBlockScalarData(dataID, mesh.data.size(), vertexIDs.data(), mesh.data.data());
-      std::ostringstream oss;
-      for(size_t i = 0; i<std::min((size_t)10, mesh.data.size()); ++i)
-          oss << mesh.data[i] << ' ';
-      VLOG(1) << "Data read: " << oss.str();
     }
     round++;
   }
@@ -297,14 +307,13 @@ int main(int argc, char* argv[])
     auto const & positions = mesh.positions;
     auto const & data = mesh.data;
     for (size_t i = 0; i < mesh.positions.size(); i++) {
-      ostream << positions[i][0] << " "
-              << positions[i][1] << " "
-              << positions[i][2] << " "
-              << data[i] << '\n';
-    }
+        ostream << positions[i][0] << " "
+                << positions[i][1] << " "
+                << positions[i][2] << " "
+                << data[i] << '\n';
+      }
     ostream.close();
   }
 
   MPI_Finalize();
 }
-