patchkernel: rework evaluation of characteristic size of three-dimensional elements

src/IO/STL.cpp

@@ -23,6 +23,7 @@
 \*---------------------------------------------------------------------------*/
 
 #include <cassert>
+#include <cstdint>
 #include <cstring>
 
 #include "bitpit_common.hpp"

src/patchkernel/element.cpp

@@ -1550,8 +1550,77 @@ double Element::evalSize(const std::array<double, 3> *coordinates) const
 {
 	switch (m_type) {
 
-	case ElementType::POLYGON:
 	case ElementType::POLYHEDRON:
+	{
+		// The characteristics size of a polyhedron is evaluated as the
+		// volume divided by the area of the largest face.
+		int nFaces = getFaceCount();
+		BITPIT_CREATE_WORKSPACE(faceVertexCoords, std::array<double BITPIT_COMMA 3>, getVertexCount(), ReferenceElementInfo::MAX_ELEM_VERTICES);
+
+		double volume = evalVolume(coordinates);
+
+		double faceMaxArea = 0;
+		for (int face = 0; face < nFaces; ++face) {
+			ElementType faceType = getFaceType(face);
+			bool faceHasReferenceInfo = ReferenceElementInfo::hasInfo(faceType);
+
+			ConstProxyVector<int> faceLocalVertexIds = getFaceLocalVertexIds(face);
+			for (int n = 0; n < getFaceVertexCount(face); ++n) {
+				faceVertexCoords[n] = coordinates[faceLocalVertexIds[n]];
+			}
+
+			double faceArea;
+			if (faceHasReferenceInfo) {
+				const Reference2DElementInfo &faceInfo = static_cast<const Reference2DElementInfo &>(ReferenceElementInfo::getInfo(faceType));
+				faceArea = faceInfo.evalArea(faceVertexCoords);
+			} else {
+				ConstProxyVector<long> faceConnect = getFaceConnect(face);
+				std::size_t faceConnectSize = faceConnect.size();
+
+				Element faceElement(0, faceType, faceConnect.size());
+				long *faceElementConnect = faceElement.getConnect();
+				for (std::size_t i = 0; i < faceConnectSize; ++i) {
+					faceElementConnect[i] = faceConnect[i];
+				}
+
+				faceArea = faceElement.evalArea(faceVertexCoords);
+			}
+			faceMaxArea = std::max(faceArea, faceMaxArea);
+		}
+
+		double size = volume / faceMaxArea;
+
+		return size;
+	}
+
+	case ElementType::POLYGON:
+	{
+		// The characteristics size of a polygon is evaluated as the
+		// area divided by the length of the longest face.
+		int nFaces = getFaceCount();
+		BITPIT_CREATE_WORKSPACE(faceVertexCoords, std::array<double BITPIT_COMMA 3>, getVertexCount(), ReferenceElementInfo::MAX_ELEM_VERTICES);
+
+		double area = evalArea(coordinates);
+
+		double faceMaxLength = 0;
+		for (int face = 0; face < nFaces; ++face) {
+			ElementType faceType = getFaceType(face);
+			assert(ReferenceElementInfo::hasInfo(faceType));
+
+			ConstProxyVector<int> faceLocalVertexIds = getFaceLocalVertexIds(face);
+			for (int n = 0; n < getFaceVertexCount(face); ++n) {
+				faceVertexCoords[n] = coordinates[faceLocalVertexIds[n]];
+			}
+
+			const Reference1DElementInfo &faceInfo = static_cast<const Reference1DElementInfo &>(ReferenceElementInfo::getInfo(faceType));
+			faceMaxLength = std::max(faceInfo.evalLength(faceVertexCoords), faceMaxLength);
+		}
+
+		double size = area / faceMaxLength;
+
+		return size;
+	}
+
 	case ElementType::UNDEFINED:
 	{
 		return 0.;