Measure computation and face normals in general dimension. (b8b6b912) · Commits · Tomáš Jakubec / GTMesh

src/GTMesh/UnstructuredMesh/MeshElements/MeshElement.h

+3 −1

Original line number	Original line	Diff line number	Diff line
	@@ -164,7 +164,9 @@ public:
	}		}

	MeshElement(IndexType index = INVALID_INDEX(IndexType))		MeshElement(IndexType index = INVALID_INDEX(IndexType))
	:MeshElementBase<IndexType>(index), CellBoundaryConnection<IndexType> () {}		:MeshElementBase<IndexType>(index),
			std::conditional<ElementDim == MeshDim - 1,CellBoundaryConnection<IndexType>, emptyStruct>::type()
			{}

	};		};

src/GTMesh/UnstructuredMesh/MeshFunctions/ComputeMeasures.h

+118 −2

Original line number	Original line	Diff line number	Diff line
	@@ -6,21 +6,98 @@
	#include "../MeshDataContainer/MeshDataContainer.h"		#include "../MeshDataContainer/MeshDataContainer.h"
	#include "../../NumericStaticArray/GramSchmidt.h"		#include "../../NumericStaticArray/GramSchmidt.h"
	#include <array>		#include <array>
			#include "GetCenters.h"


	namespace Impl {		namespace Impl {


	template <unsigned int CurrentDimension, unsigned int MeshDimension, ComputationMethod Method = METHOD_DEFAULT>		template <unsigned int CurrentDimension, unsigned int MeshDimension, ComputationMethod Method = METHOD_DEFAULT>
	struct _ComputeMeasures{		struct _ComputeMeasures{

	template <typename IndexType, typename Real, unsigned int ...Reserve>		template <typename IndexType, typename Real, unsigned int ...Reserve>
	static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>&, const MeshElements<MeshDimension, IndexType, Real, Reserve...>&){		static void compute(MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>&, const MeshElements<MeshDimension, IndexType, Real, Reserve...>&){
	static_assert (MeshDimension <= 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");		static_assert (MeshDimension <= 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
	}		}

			/**
			* @brief compute
			* @param measures
			* @param centers
			* @param mesh
			*/
			template <typename IndexType, typename Real, unsigned int ...Reserve>
			static void compute(
			MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& measures,
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
			// Calculate the measures of the elements of dimension CurrentDimension
			for(const auto& element : mesh.template getElements<CurrentDimension>()){

			Real measure = 0;

			// Prepare the lambda function to be applied in the mesh.apply
			auto calculationLambda = [&](IndexType, IndexType sube){
			auto res = getCenters<CurrentDimension - 1>(centers, mesh, sube);

			std::array<Vector<MeshDimension, Real>, CurrentDimension> v;
			for (size_t i = 1; i < res.size(); i++){
			v[i-1] = res[i] - res[0];
			}
			v.back() = centers[element] - res[0];
			std::array<Real, CurrentDimension> norms;

			gramSchmidt<CurrentDimension, MeshDimension, IndexType, Real>(v, norms);
			measure += norms.back() * measures.template getDataByDim<CurrentDimension - 1>()[sube];
			};

			// Perform the calculation
			MeshApply<CurrentDimension, CurrentDimension - 1>::apply( element.getIndex(),
			mesh,
			calculationLambda );
			measures[element] = (1.0/CurrentDimension) * measure;
			}
			_ComputeMeasures<CurrentDimension + 1, MeshDimension, Method>::compute(measures, centers, mesh);
			}
	};		};


			template < unsigned int MeshDimension, ComputationMethod Method >
			struct _ComputeMeasures<MeshDimension, MeshDimension, Method>{


			template <typename IndexType, typename Real, unsigned int ...Reserve>
			static void compute(
			MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& measures,
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
			// Calculate the measures of the elements of dimension MeshDimension
			for(const auto& element : mesh.template getElements<MeshDimension>()){

			Real measure = 0;

			// Prepare the lambda function to be applied in the mesh.apply
			auto calculationLambda = [&](IndexType, IndexType sube){
			auto res = getCenters<MeshDimension-1>(centers, mesh, sube);
			std::array<Vector<MeshDimension, Real>, MeshDimension> v;
			for (size_t i = 1; i < res.size(); i++){
			v[i-1] = res[i] - res[0];
			}
			v.back() = centers[element] - res[0];
			std::array<Real, MeshDimension> norms;

			gramSchmidt<MeshDimension, MeshDimension, IndexType, Real>(v, norms);
			measure += norms.back() * measures.template getDataByDim<MeshDimension - 1>()[sube];
			};

			// Perform the calculation
			MeshApply<MeshDimension, MeshDimension - 1>::apply( element.getIndex(),
			mesh,
			calculationLambda );
			measures[element] = (1.0/MeshDimension) * measure;
			}
			}
			};

	template <unsigned int MeshDimension, ComputationMethod Method>		template <unsigned int MeshDimension, ComputationMethod Method>
	struct _ComputeMeasures<1, MeshDimension, Method>{		struct _ComputeMeasures<1, MeshDimension, Method>{
	template <typename IndexType, typename Real, unsigned int ...Reserve>		template <typename IndexType, typename Real, unsigned int ...Reserve>
	@@ -36,6 +113,23 @@ struct _ComputeMeasures<1, MeshDimension, Method>{

	_ComputeMeasures<2, MeshDimension>::compute(measures, mesh);		_ComputeMeasures<2, MeshDimension>::compute(measures, mesh);
	}		}

			template <typename IndexType, typename Real, unsigned int ...Reserve>
			static void compute(
			MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& measures,
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){

			auto& edgeLengths = measures.template getDataByDim<1>();

			for (IndexType edgeIndex = 0; edgeIndex < mesh.template getElements<1>().size(); edgeIndex++) {
			auto& edge = mesh.getEdges().at(edgeIndex);
			edgeLengths.at(edgeIndex) = (mesh.getVertices().at(edge.getVertexAIndex()) -
			mesh.getVertices().at(edge.getVertexBIndex())).normEuclid();
			}

			_ComputeMeasures<2, MeshDimension>::compute(measures, centers, mesh);
			}
	};		};


	@@ -190,6 +284,11 @@ struct _ComputeMeasures<3, 3, METHOD_TESSELLATED>{



			/**
			* @brief Calculates measures of all elements in the mesh
			* with dimension higher than 1. This function is able to work
			* with meshes with dimension lower or equal to 3.
			*/
	template <ComputationMethod Method, unsigned int MeshDimension,typename IndexType, typename Real, unsigned int ...Reserve>		template <ComputationMethod Method, unsigned int MeshDimension,typename IndexType, typename Real, unsigned int ...Reserve>
	MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> computeMeasures(const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){		MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> computeMeasures(const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
	MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> measures(mesh);		MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> measures(mesh);
	@@ -199,5 +298,22 @@ MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, Me
	return measures;		return measures;
	}		}

			/**
			* @brief Calculates measures of all elements in the mesh
			* with dimension higher than 1. This function is able to work
			* with meshes with dimension higher than 3.
			*/
			template <ComputationMethod Method, unsigned int MeshDimension,typename IndexType, typename Real, unsigned int ...Reserve>
			MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> computeMeasures(
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh
			){
			MakeMeshDataContainer_t<Real, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>> measures(mesh);

			Impl::_ComputeMeasures<1, MeshDimension, Method>::compute(measures, centers, mesh);

			return measures;
			}


	#endif // COMPUTEMEASURES_H		#endif // COMPUTEMEASURES_H

src/GTMesh/UnstructuredMesh/MeshFunctions/ComputeNormals.h

+47 −4

Original line number	Original line	Diff line number	Diff line
	@@ -7,17 +7,48 @@
	#include "../../NumericStaticArray/GramSchmidt.h"		#include "../../NumericStaticArray/GramSchmidt.h"
	#include "MeshApply.h"		#include "MeshApply.h"
	#include "MeshFunctionsDefine.h"		#include "MeshFunctionsDefine.h"
			#include "GetCenters.h"

	namespace Impl {		namespace Impl {

	template <unsigned int Dimension, ComputationMethod Method = METHOD_DEFAULT>		template <unsigned int MeshDimension, ComputationMethod Method = METHOD_DEFAULT>
	struct _ComputeNormals{		struct _ComputeNormals{
	template <typename IndexType, typename Real, unsigned int ...Reserve>		template <typename IndexType, typename Real, unsigned int ...Reserve>
	static void compute(MeshDataContainer<Vector<Dimension, Real>, Dimension-1>&,		static void compute(MeshDataContainer<Vector<MeshDimension, Real>, MeshDimension-1>&,
	const MeshElements<Dimension, IndexType, Real, Reserve...>&){		const MeshElements<MeshDimension, IndexType, Real, Reserve...>&){
	static_assert (Dimension > 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");		static_assert (MeshDimension > 3,"The measure computation of mesh of dimension higher than 3 is not implemented yet.");
	throw std::runtime_error("The computation of face normal vectors of mesh of dimension higher than 3 is not implemented yet.");		throw std::runtime_error("The computation of face normal vectors of mesh of dimension higher than 3 is not implemented yet.");
	}		}

			template <typename IndexType, typename Real, unsigned int ...Reserve>
			static void compute(MeshDataContainer<Vector<MeshDimension, Real>, MeshDimension-1>& normals,
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){
			for (auto& face : mesh.getFaces()) {

			double vectorSign = 1.0;
			IndexType cellIndex = face.getCellLeftIndex();
			if (isInvalidIndex( cellIndex ) \|\| isBoundaryIndex( cellIndex )) {
			vectorSign = -1.0;
			cellIndex = face.getCellRightIndex();
			}

			const Vertex<MeshDimension,Real>& cellCenter = centers.template getDataByDim<MeshDimension>().at(cellIndex);

			auto res = getCenters<MeshDimension - 1>(centers, mesh, face.getIndex());

			std::array<Vector<MeshDimension, Real>, MeshDimension> gsVecs;
			for (size_t i = 1; i < res.size(); i++){
			gsVecs[i-1] = res[i] - res[0];
			}
			gsVecs.back() = cellCenter - res[0];
			std::array<Real, MeshDimension> norms;
			gramSchmidt<MeshDimension, MeshDimension, IndexType, Real>(gsVecs, norms);

			normals[face] = vectorSign * gsVecs.back();
			}
			}

	};		};


	@@ -142,4 +173,16 @@ MeshDataContainer<Vector<Dimension, Real>, Dimension-1> computeFaceNormals(const
	}		}


			template <ComputationMethod Method, unsigned int MeshDimension,typename IndexType, typename Real, unsigned int ...Reserve>
			MeshDataContainer<Vector<MeshDimension, Real>, MeshDimension-1> computeFaceNormals(
			const MakeMeshDataContainer_t<Vertex<MeshDimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){

			MeshDataContainer<Vector<MeshDimension, Real>, MeshDimension-1> normals(mesh);

			Impl::_ComputeNormals<MeshDimension, Method>::compute(normals, centers, mesh);

			return normals;
			}

	#endif // COMPUTENORMALS_H		#endif // COMPUTENORMALS_H

src/GTMesh/UnstructuredMesh/MeshFunctions/GetCenters.h

0 → 100644

+86 −0

Original line number	Original line	Diff line number	Diff line
			#ifndef GETCENTERS_H
			#define GETCENTERS_H
			#include "MeshFunctionsDefine.h"
			#include "../MeshDataContainer/MeshDataContainer.h"
			#include "../../NumericStaticArray/GramSchmidt.h"
			#include <array>

			namespace Impl {
			template < unsigned int Dimension, unsigned int CurrentDimension, unsigned int MeshDimension >
			struct GetCenters {
			template <typename IndexType, typename Real, unsigned int ...Reserve, unsigned int Dim = CurrentDimension, typename std::enable_if<Dim == 1, bool>::type = true>
			static void
			get( std::array< Vertex< MeshDimension, Real >, Dimension + 1 >& res,
			const MakeMeshDataContainer_t<Vertex< MeshDimension, Real >, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
			const IndexType& elementIndex ){


			(void) centers;
			res[Dimension - 1] = mesh.getVertices()[mesh.getEdges()[elementIndex].getVertexAIndex()];
			res[Dimension] = mesh.getVertices()[mesh.getEdges()[elementIndex].getVertexBIndex()];


			}

			template <typename IndexType, typename Real, unsigned int ...Reserve, unsigned int Dim = CurrentDimension, typename std::enable_if<Dim != MeshDimension && (Dim > 1), bool>::type = true>
			static void
			get( std::array< Vertex< MeshDimension, Real >, Dimension + 1 >& res,
			const MakeMeshDataContainer_t<Vertex< MeshDimension, Real >, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
			const IndexType& elementIndex ){

			res[Dimension - CurrentDimension] = centers.template getDataByDim<CurrentDimension>()[elementIndex];


			IndexType sube = mesh.template getElements<CurrentDimension>()[elementIndex].getSubelements()[0];
			GetCenters<Dimension, CurrentDimension - 1, MeshDimension>::get(res, centers, mesh, sube);

			}

			template <typename IndexType, typename Real, unsigned int ...Reserve, unsigned int Dim = CurrentDimension, typename std::enable_if<Dim == MeshDimension, bool>::type = true>
			static void
			get( std::array< Vertex< MeshDimension, Real >, Dimension + 1 >& res,
			const MakeMeshDataContainer_t<Vertex< MeshDimension, Real >, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
			const IndexType& elementIndex ){

			res[0] = centers.template getDataByDim<CurrentDimension>()[elementIndex];

			IndexType sube = mesh.template getElements<CurrentDimension>()[elementIndex].getBoundaryElementIndex();
			GetCenters<Dimension, CurrentDimension - 1, MeshDimension>::get(res, centers, mesh, sube);

			return res;

			}



			};


			/**
			* @brief This function returns an array of vertices
			* describing the "plane" the element is lying in.
			* The vertices are chosen as centers of subelements
			* (one per dimension) except edges. This choice prevents
			* obtaining linear dependent vertices. The processed element is
			* specified by its index @a elementIndex and dimension @a Dimension.
			* @param centers: centers of the mesh (can be obtained by the function calculateCenters)
			* @param mesh: mesh to be processed
			* @param elementIndex: index of the processed element
			*/
			template <unsigned int Dimension, unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve>
			std::array< Vertex< MeshDimension, Real >, Dimension + 1 > getCenters( const MakeMeshDataContainer_t<Vertex< MeshDimension, Real >, make_custom_integer_sequence_t<unsigned int, 1, MeshDimension>>& centers,
			const MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh,
			IndexType elementIndex ) {
			std::array< Vertex< MeshDimension, Real >, Dimension + 1 > res;
			GetCenters<Dimension, Dimension, MeshDimension>::get(res, centers, mesh, elementIndex);
			return res;
			}

			}



			#endif // GETCENTERS_H

src/GTMesh/UnstructuredMesh/MeshIO/MeshNativeType.h

+3 −1

Original line number	Original line	Diff line number	Diff line
	@@ -2,7 +2,9 @@
	#define MESHNATIVETYPE_H		#define MESHNATIVETYPE_H
	template <unsigned int MeshDimension>		template <unsigned int MeshDimension>
	struct MeshNativeType{		struct MeshNativeType{
			enum ElementType{
			POLYTOPE = 1000
			};
	};		};

	//! Element types of GTMesh in 2D mesh		//! Element types of GTMesh in 2D mesh