add build option to disable the viewer and use the system vtk

CMakeLists.txt

@@ -2,7 +2,37 @@ cmake_minimum_required(VERSION 3.15)
 
 project(pfe)
 
+set(VTK_COMPONENTS
+  VTK::CommonCore
+  VTK::IOCore
+  VTK::FiltersCore
+  VTK::CommonColor
+  VTK::CommonDataModel
+  VTK::IOLegacy
+  VTK::IOXML)
+set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.")
+if(ENABLE_VIEWER)
+  list(APPEND VTK_COMPONENTS
+	VTK::RenderingCore
+	VTK::ViewsCore
+	VTK::GUISupportQt
+	VTK::RenderingQt
+	VTK::ViewsQt
+	VTK::RenderingVolume
+	VTK::RenderingVolumeOpenGL2)
+endif()
+
+set(USE_SYSTEM_VTK NO CACHE BOOL
+  "Use the version of vtk installed in the system instead of downloading and compiling it ourselves.")
+if(USE_SYSTEM_VTK)
+  list(TRANSFORM VTK_COMPONENTS REPLACE "VTK::" ""
+	OUTPUT_VARIABLE VTK_PACKAGE_COMPONENTS)
+  message("VTK_COMPONENTS: ${VTK_COMPONENTS}")
+  message("VTK_PACKAGE_COMPONENTS: ${VTK_PACKAGE_COMPONENTS}")
+  find_package(VTK COMPONENTS ${VTK_PACKAGE_COMPONENTS})
+else()
   add_subdirectory(external/vtk)
+endif()
 
 add_executable(pfe)
 
@@ -15,20 +45,12 @@ target_sources(pfe PRIVATE
   src/aspect_ratio_filter.cc
   src/aspect_ratio_filter.h
   src/dihedral_angles_filter.cc
-  src/dihedral_angles_filter.h)
+  src/dihedral_angles_filter.h
+  src/external_points_filter.cc
+  src/external_points_filter.h)
 
-target_link_libraries(pfe PRIVATE
+target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS})
-  VTK::CommonCore
+
-  VTK::ViewsCore
+if(ENABLE_VIEWER)
-  VTK::RenderingCore
+  target_compile_definitions(pfe PRIVATE ENABLE_VIEWER)
-  VTK::IOCore
+endif()
-  VTK::FiltersCore
-  VTK::CommonColor
-  VTK::GUISupportQt
-  VTK::RenderingQt
-  VTK::ViewsQt
-  VTK::RenderingVolume
-  VTK::CommonDataModel
-  VTK::IOLegacy
-  VTK::IOXML
-  VTK::RenderingVolumeOpenGL2)

data/TetMesh2.vtk

@@ -18,6 +18,3 @@ LOOKUP_TABLE default
 1
 1
 1
-FIELD law 1
-cell_law 1 1 int
-0

src/external_points_filter.cc

@@ -0,0 +1,93 @@
+#include "external_points_filter.h"
+
+#include <vtkIdList.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCellIterator.h>
+#include <vtkTriangle.h>
+#include <vtkStaticCellLinks.h>
+#include <vtkIntArray.h>
+
+
+
+/* Remove from ids_a the elements that are not also in ids_b. */
+static void keep_only_dupes(std::vector<vtkIdType> &ids_a,
+                            vtkIdType n_ids_b, const vtkIdType *ids_b) {
+	for (auto it = ids_a.begin(); it != ids_a.end();) {
+		bool duplicate = false;
+		for (vtkIdType i = 0; i < n_ids_b; i++) {
+			if (*it == ids_b[i]) {
+				duplicate = true;
+				break;
+			}
+		}
+		if (!duplicate) {
+			it = ids_a.erase(it);
+		} else {
+			++it;
+		}
+	}
+}
+
+
+vtkStandardNewMacro(ExternalPointsFilter);
+
+
+vtkTypeBool ExternalPointsFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+
+	/* Copy point and cell data unmodified. */
+	output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+	output->GetCellData()->PassData(input->GetCellData());
+
+	/* Create an array to store the IDs of external points. */
+	vtkNew<vtkIdTypeArray> externalPoints;
+	externalPoints->SetName("external_points");
+	externalPoints->SetNumberOfComponents(1);
+	vtkNew<vtkIntArray> isExternal;
+	isExternal->SetName("bl");
+	isExternal->SetNumberOfComponents(1);
+	isExternal->SetNumberOfTuples(input->GetNumberOfPoints());
+	isExternal->FillValue(0);
+
+	/* Generate cell links, these tell us which cell a point belongs to. */
+	vtkNew<vtkStaticCellLinks> links;
+	links->BuildLinks(input);
+
+	// TODO: try to use vtkDataSet::GetCellNeighbors instead
+	auto *it = input->NewCellIterator();
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		vtkIdList *pointIds = it->GetPointIds();
+		for (vtkIdType faceId = 0; faceId < 4; faceId++) {
+			vtkIdType firstPoint = pointIds->GetId(faceId);
+			vtkIdType n_cells = links->GetNcells(firstPoint);
+			vtkIdType *cells = links->GetCells(firstPoint);
+			std::vector<vtkIdType> intersection(cells, cells + n_cells);
+			for (vtkIdType i = 1; i < 3; i++) {
+				vtkIdType pointId = pointIds->GetId((faceId + i) % 4);
+				keep_only_dupes(intersection,
+				                links->GetNcells(pointId), links->GetCells(pointId));
+			}
+			if (intersection.size() == 1) { // The face only belongs to one cell.
+				for (vtkIdType i = 0; i < 3; i++) {
+					externalPoints->InsertNextValue(pointIds->GetId((faceId + i) % 4));
+					isExternal->SetValue(pointIds->GetId((faceId + i) % 4), 1);
+				}
+			}
+		}
+	}
+
+	output->GetPointData()->SetScalars(isExternal);
+	output->GetFieldData()->AddArray(externalPoints);
+	return true;
+}

src/external_points_filter.h

@@ -0,0 +1,22 @@
+#ifndef EXTERNAL_POINTS_FILTER_H
+#define EXTERNAL_POINTS_FILTER_H
+
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+
+class ExternalPointsFilter : public vtkUnstructuredGridAlgorithm {
+public:
+	static ExternalPointsFilter *New();
+	vtkTypeMacro(ExternalPointsFilter, vtkUnstructuredGridAlgorithm);
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+
+protected:
+	~ExternalPointsFilter() override = default;
+};
+
+
+#endif

src/main.cc

@@ -1,11 +1,16 @@
 #include "angles_filter.h"
 #include "aspect_ratio_filter.h"
 #include "dihedral_angles_filter.h"
+#include "external_points_filter.h"
 
 #include <vtkCellData.h>
 #include <vtkUnstructuredGrid.h>
-#include <vtkCamera.h>
+#include <vtkUnstructuredGridReader.h>
+#include <vtkUnstructuredGridWriter.h>
+
+#ifdef USE_VIEWER
 #include <vtkNamedColors.h>
+#include <vtkCamera.h>
 #include <vtkVolumeProperty.h>
 #include <vtkRenderWindow.h>
 #include <vtkRenderWindowInteractor.h>
@@ -13,33 +18,10 @@
 #include <vtkVolumeMapper.h>
 #include <vtkVolume.h>
 #include <vtkOpenGLProjectedTetrahedraMapper.h>
-#include <vtkUnstructuredGridReader.h>
 #include <vtkPiecewiseFunction.h>
+#endif
 
 #include <array>
-#include <vector>
-#include <algorithm>
-
-
-// template <typename T>
-// T average(const std::vector<T> &data) {
-// 	T avg = 0;
-// 	for (const T &t : data) {
-// 		avg += t;
-// 	}
-// 	return avg / data.size();
-// }
-
-
-// template <typename T>
-// T standard_deviation(const std::vector<T> &data) {
-// 	T avg = average(data);
-// 	T stddev = 0;
-// 	for (const T &t : data) {
-// 		stddev += (t - avg) * (t - avg);
-// 	}
-// 	return std::sqrt(stddev / data.size());
-// }
 
 
 int main(int argc, char **argv) {
@@ -58,49 +40,65 @@ int main(int argc, char **argv) {
 	anglesFilter->SetInputConnection(reader->GetOutputPort());
 
 	/* Display angles in console. */
-	anglesFilter->Update();
+	// anglesFilter->Update();
-	vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
+	// vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
-	auto *angles = grid->GetCellData()->GetArray("angles");
+	// auto *angles = grid->GetCellData()->GetArray("angles");
-	for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
+	// for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": ";
+	// 	std::cerr << "cell " << i << ": ";
-		for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
+	// 	for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
-			std::cout << angles->GetTuple(i)[j] << ", ";
+	// 		std::cerr << angles->GetTuple(i)[j] << ", ";
-		}
+	// 	}
-		std::cout << "\b\b \n";
+	// 	std::cerr << "\b\b \n";
-	}
+	// }
 
 	/* Aspect ratio */
 	vtkNew<AspectRatioFilter> aspectRatioFilter;
 	aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort());
 
 	/* Display aspect ratios in console. */
-	aspectRatioFilter->Update();
+	// aspectRatioFilter->Update();
-	grid = aspectRatioFilter->GetOutput();
+	// grid = aspectRatioFilter->GetOutput();
-	auto *aspectRatios = grid->GetCellData()->GetArray("aspect_ratio");
+	// auto *aspectRatios = grid->GetCellData()->GetArray("aspect_ratio");
-	for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
+	// for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": "
+	// 	std::cerr << "cell " << i << ": "
-				  << aspectRatios->GetTuple1(i) << "\n";
+	// 			  << aspectRatios->GetTuple1(i) << "\n";
-	}
+	// }
 
 	/* Dihedral angles filter */
 	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
 	dihedralAnglesFilter->SetInputConnection(aspectRatioFilter->GetOutputPort());
 
 	/* Display dihedral angles in console. */
-	dihedralAnglesFilter->Update();
+	// dihedralAnglesFilter->Update();
-	grid = dihedralAnglesFilter->GetOutput();
+	// grid = dihedralAnglesFilter->GetOutput();
-	auto dihedralAngles = grid->GetCellData()->GetArray("dihedral_angles");
+	// auto dihedralAngles = grid->GetCellData()->GetArray("dihedral_angles");
-	for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++) {
+	// for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++) {
-		std::cout << "dihedral angles ";
+	// 	std::cerr << "dihedral angles ";
-		for (vtkIdType j = 0; j < 6; j++) {
+	// 	for (vtkIdType j = 0; j < 6; j++) {
-			std::cout << dihedralAngles->GetTuple(i)[j] << ", ";
+	// 		std::cerr << dihedralAngles->GetTuple(i)[j] << ", ";
-		}
+	// 	}
-		std::cout << "\b\b\n";
+	// 	std::cerr << "\b\b\n";
-	}
+	// }
 
+	/* External faces filter. */
+	// vtkNew<vtkmExternalFaces> externalFacesFilter;
+	// externalFacesFilter->DebugOn();
+	// externalFacesFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+
+	/* External points filter. */
+	vtkNew<ExternalPointsFilter> externalPointsFilter;
+	externalPointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+
+	vtkNew<vtkUnstructuredGridWriter> writer;
+	writer->SetInputConnection(externalPointsFilter->GetOutputPort());
+	writer->SetFileTypeToASCII();
+	writer->SetFileName("out.vtk");
+	writer->Write();
+
+#ifdef USE_VIEWER
 	/* Volume rendering properties */
 	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
-	volumeMapper->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+	volumeMapper->SetInputConnection(externalPointsFilter->GetOutputPort());
 	vtkNew<vtkVolume> volume;
 	volume->SetMapper(volumeMapper);
 	vtkNew<vtkPiecewiseFunction> transferFunction;
@@ -126,6 +124,6 @@ int main(int argc, char **argv) {
 	renderWindowInteractor->SetRenderWindow(renderWindow);
 	renderWindow->Render();
 	renderWindowInteractor->Start();
-
+#endif
 	return EXIT_SUCCESS;
 }