Loading Unstructured_mesh/MeshDataContainer.h +41 −12 Original line number Diff line number Diff line Loading @@ -2,7 +2,22 @@ #define MESHDATACONTAINER_H #include "MeshElement.h" #include "../debug/debug.h" #include "../debug/Debug.h" template<typename DataType, unsigned int Position, unsigned int MappedDimenion> struct DataContainer : public std::vector<DataType> { using type = DataType; constexpr unsigned int getPosition() { return Position; } constexpr unsigned int getMappedDimension() { return MappedDimenion; } }; /** Loading @@ -29,16 +44,22 @@ private: return DimensionPos<dim, 0, std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>::res(); } template<unsigned int pos> static constexpr unsigned int dimensionAt(){ return std::get<pos>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...}); } public: template<typename _DataType, unsigned int _Dim> struct _DataContainer : _DataContainer<_DataType,_Dim - 1> { std::vector<_DataType> _data; DataContainer<_DataType, _Dim, dimensionAt<_Dim>()> _data; }; template<typename _DataType> struct _DataContainer<_DataType, 0>{ std::vector<_DataType> _data; struct _DataContainer<_DataType, 0> : public std::vector<_DataType>{ DataContainer<_DataType, 0, dimensionAt<0U>()> _data; }; private: template<unsigned int pos, typename dummy = void> struct Alocator{ MeshDataContainer<DataType, Dimensions...>& parent; Loading Loading @@ -95,13 +116,13 @@ private: public: template<unsigned int dim> std::vector<DataType>& getDataByDim(){ DataContainer<DataType, dimensionIndex<dim>(), dim>& getDataByDim(){ return data._DataContainer<DataType, dimensionIndex<dim>()>::_data; } template<unsigned int pos> std::vector<DataType>& getDataByPos(){ DataContainer<DataType, pos, dimensionAt<pos>()>& getDataByPos(){ return data._DataContainer<DataType,pos>::_data; } Loading Loading @@ -184,6 +205,11 @@ public: static constexpr unsigned int dimensionIndex(){ return DimensionPos<dim, 0, std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>::res(); } template<unsigned int pos> static constexpr unsigned int dimensionAt(){ return std::get<pos>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...}); } public: template<unsigned int pos> Loading @@ -191,12 +217,14 @@ public: template<unsigned int Pos, typename Dummy = void> struct _DataContainer : _DataContainer<Pos - 1, Dummy>{ std::vector<DataType<Pos>> _data; DataContainer<DataType<Pos>, Pos, dimensionAt<Pos>()> _data; //std::vector<DataType<Pos>> _data; }; template<typename Dummy> struct _DataContainer<0, Dummy>{ std::vector<DataType<0>> _data; DataContainer<DataType<0>, 0, dimensionAt<0>()> _data; //std::vector<DataType<0>> _data; }; template<unsigned int pos, typename _DataType, typename... _DataTypes> Loading Loading @@ -260,7 +288,8 @@ public: * @return */ template<unsigned int dim> std::vector<std::tuple_element_t<dimensionIndex<dim>(), std::tuple<DataTypes...>>>& getDataByDim(){ DataContainer<std::tuple_element_t<dimensionIndex<dim>(), std::tuple<DataTypes...>>, dimensionIndex<dim>(), dim>& getDataByDim(){ return data._DataContainer<dimensionIndex<dim>()>::_data; } Loading @@ -270,7 +299,7 @@ public: * @return */ template<unsigned int pos> std::vector<DataType<pos>>& getDataByPos(){ DataContainer<DataType<pos>, pos, dimensionAt<pos>()>& getDataByPos(){ return data._DataContainer<pos>::_data; } Loading Unstructured_mesh/MeshElement.h +2 −2 Original line number Diff line number Diff line Loading @@ -102,7 +102,7 @@ public: return this->begin() + getNumberOfSubElements(); } typename std::array<Subelement<IndexType>, Reserve>::const_iterator cend(){ typename std::array<Subelement<IndexType>, Reserve>::const_iterator cend() const { return this->cbegin() + getNumberOfSubElements(); } }; Loading Unstructured_mesh/MeshFunctions.h +99 −19 Original line number Diff line number Diff line Loading @@ -447,9 +447,10 @@ MeshDataContainer<Real, Dimension-1> ComputeCellsDistance(MeshElements<Dimension namespace temp1 { template <unsigned int CurrentDimension, unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool End, bool Ascend> template <unsigned int CurrentDimension, unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool End, bool Descend> struct MeshRun { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> Loading @@ -461,14 +462,14 @@ struct MeshRun { auto i = mesh.template getElements<CurrentDimension>().at(index); for (auto sube: mesh.template getElement<CurrentDimension>(i.getIndex()).getSubelements()) MeshRun< CurrentDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == CurrentDimension - 1, Ascend>::run(mesh, origElementIndex, sube.index, fun); MeshRun< CurrentDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == CurrentDimension - 1, Descend>::run(mesh, origElementIndex, sube.index, fun); } }; template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, false, Ascend> { template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, false, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, Loading @@ -479,7 +480,7 @@ struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, fa auto& cell = mesh.getCells().at(index); IndexType tmpFace = cell.getBoundaryElementIndex(); do { MeshRun<MeshDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == MeshDimension - 1, Ascend>::run(mesh, origElementIndex, tmpFace, fun); MeshRun<MeshDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == MeshDimension - 1, Descend>::run(mesh, origElementIndex, tmpFace, fun); tmpFace = mesh.getFaces().at(tmpFace).getNextBElem(cell.getIndex()); } while (tmpFace != cell.getBoundaryElementIndex()); Loading @@ -487,8 +488,8 @@ struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, fa }; template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Ascend> { template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, Loading @@ -497,25 +498,25 @@ struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Ascend> Func fun){ auto& edge = mesh.getEdges().at(index); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Ascend>::run(mesh, origElementIndex, edge.getVertexAIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Ascend>::run(mesh, origElementIndex, edge.getVertexBIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Descend>::run(mesh, origElementIndex, edge.getVertexAIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Descend>::run(mesh, origElementIndex, edge.getVertexBIndex(), fun); } }; template <unsigned int CurrentDimension,unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<CurrentDimension, StartDimension, TargetDimension, MeshDimension, true, Ascend> { template <unsigned int CurrentDimension,unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<CurrentDimension, StartDimension, TargetDimension, MeshDimension, true, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& , IndexType origElementIndex, IndexType index, Func fun){ if(Ascend){ fun(StartDimension, TargetDimension, origElementIndex, index); if(Descend){ fun(origElementIndex, index); }else{ fun(TargetDimension, StartDimension, index, origElementIndex); fun(index, origElementIndex); } } }; Loading Loading @@ -549,7 +550,7 @@ struct MeshConnections { MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh ) { MeshDataContainer<std::set<IndexType>, StartDim> result(mesh); MeshApply<StartDim, TargetDim, MeshDimension>::apply(mesh, [&result](unsigned int, unsigned int, IndexType ori, IndexType element){ MeshApply<StartDim, TargetDim, MeshDimension>::apply(mesh, [&result](IndexType ori, IndexType element){ result.template getDataByPos<0>().at(ori).insert(element); }); Loading Loading @@ -578,7 +579,7 @@ struct MeshColouring { run(mesh, startElement.getIndex(), startElement.getIndex(), [&possibleColours, &attachedColours](unsigned int, unsigned int, IndexType, IndexType element){ [&possibleColours, &attachedColours](IndexType, IndexType element){ DBGTRY(possibleColours &= !attachedColours.template getDataByPos<0>().at(element);) } ); Loading @@ -593,7 +594,7 @@ struct MeshColouring { run(mesh, startElement.getIndex(), startElement.getIndex(), [selectedColour, &attachedColours](unsigned int, unsigned int, IndexType, IndexType element){ [selectedColour, &attachedColours](IndexType, IndexType element){ DBGTRY(attachedColours.template getDataByPos<0>().at(element)[selectedColour] = true;) } ); Loading Loading @@ -652,7 +653,86 @@ static MeshDataContainer<unsigned int, FromDim> colour( } }; template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<2, IndexType, Real, Reserve...>& mesh, typename MeshElements<2, IndexType, Real, Reserve...>::template ElementType<2>& face, typename MeshElements<2, IndexType, Real, Reserve...>::Edge& edge ) { Vertex<2, Real> AminC = mesh.getVertices().at(edge.getVertexAIndex()) - face.getCenter(); Vertex<2, Real> BminC = mesh.getVertices().at(edge.getVertexBIndex()) - face.getCenter(); Real res = AminC[0]*BminC[1]-BminC[0]*AminC[1]; return res > 0; throw std::runtime_error("can not determine orientation of edge " + std::to_string(edge.getIndex()) + " wrt face: " + std::to_string(face.getIndex())); } template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<3, IndexType, Real, Reserve...>& mesh, typename MeshElements<3, IndexType, Real, Reserve...>::template ElementType<2>& face, typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge, Vector<3, Real> faceNormal ) { Vertex<3, Real> AminC = mesh.getVertices().at(edge.getVertexAIndex()) - face.getCenter(); Vertex<3, Real> BminC = mesh.getVertices().at(edge.getVertexBIndex()) - face.getCenter(); Real res = Real(0); for (IndexType i = 0; i < 3; i++){ IndexType ipo = (i+1)%(3); IndexType ipt = (i+2)%(3); res += AminC[i]*BminC[ipo]*faceNormal[ipt]-AminC[ipt]*BminC[ipo]*faceNormal[i]; } return res > 0; } template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, IndexType faceIndex, IndexType edgeIndex) { typename MeshElements<MeshDimension, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(edgeIndex); typename MeshElements<MeshDimension, IndexType, Real, Reserve...>::template ElementType<2>& face = mesh.template getElements<2>().at(faceIndex); return edgeIsLeft(mesh, face, edge); } template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<3, IndexType, Real, Reserve...>& mesh, IndexType faceIndex, IndexType edgeIndex) { typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(edgeIndex); typename MeshElements<3, IndexType, Real, Reserve...>::template ElementType<2>& face = mesh.template getElements<2>().at(faceIndex); auto normals = ComputeFaceNormals(mesh); return edgeIsLeft(mesh, face, edge, normals[face]); } template<typename IndexType, typename Real, unsigned int ...Reserve> MeshDataContainer<std::vector<bool>, 2> edgesOrientation(MeshElements<3, IndexType, Real, Reserve...>& mesh) { MeshDataContainer<std::vector<bool>, 2> orientations(mesh); auto normals = ComputeFaceNormals(mesh); for (auto& face : mesh.getFaces()) { orientations[face].resize(face.getSubelements().getNumberOfSubElements()); for (IndexType i = 0; i < face.getSubelements().getNumberOfSubElements(); i++){ typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(face.getSubelements()[i].index); orientations[face][i] = edgeIsLeft(mesh, face, edge, normals[face]); } } return orientations; } #endif // MESH_FUNCTIONS_H Unstructured_mesh/MeshReader.h +8 −0 Original line number Diff line number Diff line Loading @@ -14,5 +14,13 @@ public: using type = MeshNativeType<2>; }; template <typename IndexType, typename Real> class MeshReader<3, IndexType, Real> { public: using type = MeshNativeType<3>; }; #endif // MESHREADER_H Unstructured_mesh/MeshWriter.h +56 −7 Original line number Diff line number Diff line #ifndef MESHWRITER_H #define MESHWRITER_H #include "MeshNativeType.h" #include "Vertex.h" #include "MeshElement.h" template<unsigned int MeshDimension, typename IndexType, typename Real> class MeshWriter{ }; enum { hmmm protected: /** * @brief The MeshHash struct<HR> * A container to store cumulative data about a mesh * to recognize changes in the mesh. * It uses number of elements and sum of vertices coordinations. */ struct MeshHash { IndexType numberOfElements = 0; Real totalVert = 0; bool operator== (const MeshHash& rhs) const { return numberOfElements == rhs.numberOfElements && totalVert == rhs.totalVert; } bool operator!= (const MeshHash& rhs) const { return !((*this)==rhs); } }; template <typename IndexType, typename Real> class MeshWriter<2, IndexType, Real> { private: template<unsigned int Dim, typename Dummy = void> struct sumOfMeshElements{ template<unsigned int ...Reserve> static IndexType sum(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ return mesh.template getElements<Dim>().size() + sumOfMeshElements<Dim - 1>::sum(mesh); } }; template<typename Dummy> struct sumOfMeshElements<0, Dummy>{ template<unsigned int ...Reserve> static IndexType sum(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ return mesh.template getElements<0>().size(); } }; public: using type = MeshNativeType<2>; using type = MeshNativeType<MeshDimension>; template<unsigned int ...Reserve> static MeshHash computeHash(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ MeshHash res; // A vector of ones that simplifies the sum of coordinates Vertex<MeshDimension, Real> ones; for(unsigned int dim = 0; dim < MeshDimension; dim++){ ones[dim] = 1.0; } for(IndexType i = 0; i < mesh.getVertices().size(); i++) { res.totalVert += mesh.getVertices().at(i) * ones; } res.numberOfElements = sumOfMeshElements<MeshDimension>::sum(mesh); return res; } }; #endif // MESHWRITER_H Loading
Unstructured_mesh/MeshDataContainer.h +41 −12 Original line number Diff line number Diff line Loading @@ -2,7 +2,22 @@ #define MESHDATACONTAINER_H #include "MeshElement.h" #include "../debug/debug.h" #include "../debug/Debug.h" template<typename DataType, unsigned int Position, unsigned int MappedDimenion> struct DataContainer : public std::vector<DataType> { using type = DataType; constexpr unsigned int getPosition() { return Position; } constexpr unsigned int getMappedDimension() { return MappedDimenion; } }; /** Loading @@ -29,16 +44,22 @@ private: return DimensionPos<dim, 0, std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>::res(); } template<unsigned int pos> static constexpr unsigned int dimensionAt(){ return std::get<pos>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...}); } public: template<typename _DataType, unsigned int _Dim> struct _DataContainer : _DataContainer<_DataType,_Dim - 1> { std::vector<_DataType> _data; DataContainer<_DataType, _Dim, dimensionAt<_Dim>()> _data; }; template<typename _DataType> struct _DataContainer<_DataType, 0>{ std::vector<_DataType> _data; struct _DataContainer<_DataType, 0> : public std::vector<_DataType>{ DataContainer<_DataType, 0, dimensionAt<0U>()> _data; }; private: template<unsigned int pos, typename dummy = void> struct Alocator{ MeshDataContainer<DataType, Dimensions...>& parent; Loading Loading @@ -95,13 +116,13 @@ private: public: template<unsigned int dim> std::vector<DataType>& getDataByDim(){ DataContainer<DataType, dimensionIndex<dim>(), dim>& getDataByDim(){ return data._DataContainer<DataType, dimensionIndex<dim>()>::_data; } template<unsigned int pos> std::vector<DataType>& getDataByPos(){ DataContainer<DataType, pos, dimensionAt<pos>()>& getDataByPos(){ return data._DataContainer<DataType,pos>::_data; } Loading Loading @@ -184,6 +205,11 @@ public: static constexpr unsigned int dimensionIndex(){ return DimensionPos<dim, 0, std::get<0>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...})>::res(); } template<unsigned int pos> static constexpr unsigned int dimensionAt(){ return std::get<pos>(std::array<unsigned int, sizeof... (Dimensions)>{Dimensions...}); } public: template<unsigned int pos> Loading @@ -191,12 +217,14 @@ public: template<unsigned int Pos, typename Dummy = void> struct _DataContainer : _DataContainer<Pos - 1, Dummy>{ std::vector<DataType<Pos>> _data; DataContainer<DataType<Pos>, Pos, dimensionAt<Pos>()> _data; //std::vector<DataType<Pos>> _data; }; template<typename Dummy> struct _DataContainer<0, Dummy>{ std::vector<DataType<0>> _data; DataContainer<DataType<0>, 0, dimensionAt<0>()> _data; //std::vector<DataType<0>> _data; }; template<unsigned int pos, typename _DataType, typename... _DataTypes> Loading Loading @@ -260,7 +288,8 @@ public: * @return */ template<unsigned int dim> std::vector<std::tuple_element_t<dimensionIndex<dim>(), std::tuple<DataTypes...>>>& getDataByDim(){ DataContainer<std::tuple_element_t<dimensionIndex<dim>(), std::tuple<DataTypes...>>, dimensionIndex<dim>(), dim>& getDataByDim(){ return data._DataContainer<dimensionIndex<dim>()>::_data; } Loading @@ -270,7 +299,7 @@ public: * @return */ template<unsigned int pos> std::vector<DataType<pos>>& getDataByPos(){ DataContainer<DataType<pos>, pos, dimensionAt<pos>()>& getDataByPos(){ return data._DataContainer<pos>::_data; } Loading
Unstructured_mesh/MeshElement.h +2 −2 Original line number Diff line number Diff line Loading @@ -102,7 +102,7 @@ public: return this->begin() + getNumberOfSubElements(); } typename std::array<Subelement<IndexType>, Reserve>::const_iterator cend(){ typename std::array<Subelement<IndexType>, Reserve>::const_iterator cend() const { return this->cbegin() + getNumberOfSubElements(); } }; Loading
Unstructured_mesh/MeshFunctions.h +99 −19 Original line number Diff line number Diff line Loading @@ -447,9 +447,10 @@ MeshDataContainer<Real, Dimension-1> ComputeCellsDistance(MeshElements<Dimension namespace temp1 { template <unsigned int CurrentDimension, unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool End, bool Ascend> template <unsigned int CurrentDimension, unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool End, bool Descend> struct MeshRun { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> Loading @@ -461,14 +462,14 @@ struct MeshRun { auto i = mesh.template getElements<CurrentDimension>().at(index); for (auto sube: mesh.template getElement<CurrentDimension>(i.getIndex()).getSubelements()) MeshRun< CurrentDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == CurrentDimension - 1, Ascend>::run(mesh, origElementIndex, sube.index, fun); MeshRun< CurrentDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == CurrentDimension - 1, Descend>::run(mesh, origElementIndex, sube.index, fun); } }; template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, false, Ascend> { template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, false, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, Loading @@ -479,7 +480,7 @@ struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, fa auto& cell = mesh.getCells().at(index); IndexType tmpFace = cell.getBoundaryElementIndex(); do { MeshRun<MeshDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == MeshDimension - 1, Ascend>::run(mesh, origElementIndex, tmpFace, fun); MeshRun<MeshDimension - 1, StartDimension, TargetDimension, MeshDimension, TargetDimension == MeshDimension - 1, Descend>::run(mesh, origElementIndex, tmpFace, fun); tmpFace = mesh.getFaces().at(tmpFace).getNextBElem(cell.getIndex()); } while (tmpFace != cell.getBoundaryElementIndex()); Loading @@ -487,8 +488,8 @@ struct MeshRun<MeshDimension, StartDimension, TargetDimension, MeshDimension, fa }; template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Ascend> { template <unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, Loading @@ -497,25 +498,25 @@ struct MeshRun<1, StartDimension, TargetDimension, MeshDimension, false, Ascend> Func fun){ auto& edge = mesh.getEdges().at(index); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Ascend>::run(mesh, origElementIndex, edge.getVertexAIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Ascend>::run(mesh, origElementIndex, edge.getVertexBIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Descend>::run(mesh, origElementIndex, edge.getVertexAIndex(), fun); MeshRun<0, StartDimension, TargetDimension, MeshDimension, TargetDimension == 0, Descend>::run(mesh, origElementIndex, edge.getVertexBIndex(), fun); } }; template <unsigned int CurrentDimension,unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Ascend> struct MeshRun<CurrentDimension, StartDimension, TargetDimension, MeshDimension, true, Ascend> { template <unsigned int CurrentDimension,unsigned int StartDimension, unsigned int TargetDimension, unsigned int MeshDimension, bool Descend> struct MeshRun<CurrentDimension, StartDimension, TargetDimension, MeshDimension, true, Descend> { template<typename Func, typename IndexType, typename Real, unsigned int ...Reserve> static void run(MeshElements<MeshDimension, IndexType, Real, Reserve...>& , IndexType origElementIndex, IndexType index, Func fun){ if(Ascend){ fun(StartDimension, TargetDimension, origElementIndex, index); if(Descend){ fun(origElementIndex, index); }else{ fun(TargetDimension, StartDimension, index, origElementIndex); fun(index, origElementIndex); } } }; Loading Loading @@ -549,7 +550,7 @@ struct MeshConnections { MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh ) { MeshDataContainer<std::set<IndexType>, StartDim> result(mesh); MeshApply<StartDim, TargetDim, MeshDimension>::apply(mesh, [&result](unsigned int, unsigned int, IndexType ori, IndexType element){ MeshApply<StartDim, TargetDim, MeshDimension>::apply(mesh, [&result](IndexType ori, IndexType element){ result.template getDataByPos<0>().at(ori).insert(element); }); Loading Loading @@ -578,7 +579,7 @@ struct MeshColouring { run(mesh, startElement.getIndex(), startElement.getIndex(), [&possibleColours, &attachedColours](unsigned int, unsigned int, IndexType, IndexType element){ [&possibleColours, &attachedColours](IndexType, IndexType element){ DBGTRY(possibleColours &= !attachedColours.template getDataByPos<0>().at(element);) } ); Loading @@ -593,7 +594,7 @@ struct MeshColouring { run(mesh, startElement.getIndex(), startElement.getIndex(), [selectedColour, &attachedColours](unsigned int, unsigned int, IndexType, IndexType element){ [selectedColour, &attachedColours](IndexType, IndexType element){ DBGTRY(attachedColours.template getDataByPos<0>().at(element)[selectedColour] = true;) } ); Loading Loading @@ -652,7 +653,86 @@ static MeshDataContainer<unsigned int, FromDim> colour( } }; template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<2, IndexType, Real, Reserve...>& mesh, typename MeshElements<2, IndexType, Real, Reserve...>::template ElementType<2>& face, typename MeshElements<2, IndexType, Real, Reserve...>::Edge& edge ) { Vertex<2, Real> AminC = mesh.getVertices().at(edge.getVertexAIndex()) - face.getCenter(); Vertex<2, Real> BminC = mesh.getVertices().at(edge.getVertexBIndex()) - face.getCenter(); Real res = AminC[0]*BminC[1]-BminC[0]*AminC[1]; return res > 0; throw std::runtime_error("can not determine orientation of edge " + std::to_string(edge.getIndex()) + " wrt face: " + std::to_string(face.getIndex())); } template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<3, IndexType, Real, Reserve...>& mesh, typename MeshElements<3, IndexType, Real, Reserve...>::template ElementType<2>& face, typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge, Vector<3, Real> faceNormal ) { Vertex<3, Real> AminC = mesh.getVertices().at(edge.getVertexAIndex()) - face.getCenter(); Vertex<3, Real> BminC = mesh.getVertices().at(edge.getVertexBIndex()) - face.getCenter(); Real res = Real(0); for (IndexType i = 0; i < 3; i++){ IndexType ipo = (i+1)%(3); IndexType ipt = (i+2)%(3); res += AminC[i]*BminC[ipo]*faceNormal[ipt]-AminC[ipt]*BminC[ipo]*faceNormal[i]; } return res > 0; } template<unsigned int MeshDimension, typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh, IndexType faceIndex, IndexType edgeIndex) { typename MeshElements<MeshDimension, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(edgeIndex); typename MeshElements<MeshDimension, IndexType, Real, Reserve...>::template ElementType<2>& face = mesh.template getElements<2>().at(faceIndex); return edgeIsLeft(mesh, face, edge); } template<typename IndexType, typename Real, unsigned int ...Reserve> bool edgeIsLeft(MeshElements<3, IndexType, Real, Reserve...>& mesh, IndexType faceIndex, IndexType edgeIndex) { typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(edgeIndex); typename MeshElements<3, IndexType, Real, Reserve...>::template ElementType<2>& face = mesh.template getElements<2>().at(faceIndex); auto normals = ComputeFaceNormals(mesh); return edgeIsLeft(mesh, face, edge, normals[face]); } template<typename IndexType, typename Real, unsigned int ...Reserve> MeshDataContainer<std::vector<bool>, 2> edgesOrientation(MeshElements<3, IndexType, Real, Reserve...>& mesh) { MeshDataContainer<std::vector<bool>, 2> orientations(mesh); auto normals = ComputeFaceNormals(mesh); for (auto& face : mesh.getFaces()) { orientations[face].resize(face.getSubelements().getNumberOfSubElements()); for (IndexType i = 0; i < face.getSubelements().getNumberOfSubElements(); i++){ typename MeshElements<3, IndexType, Real, Reserve...>::Edge& edge = mesh.getEdges().at(face.getSubelements()[i].index); orientations[face][i] = edgeIsLeft(mesh, face, edge, normals[face]); } } return orientations; } #endif // MESH_FUNCTIONS_H
Unstructured_mesh/MeshReader.h +8 −0 Original line number Diff line number Diff line Loading @@ -14,5 +14,13 @@ public: using type = MeshNativeType<2>; }; template <typename IndexType, typename Real> class MeshReader<3, IndexType, Real> { public: using type = MeshNativeType<3>; }; #endif // MESHREADER_H
Unstructured_mesh/MeshWriter.h +56 −7 Original line number Diff line number Diff line #ifndef MESHWRITER_H #define MESHWRITER_H #include "MeshNativeType.h" #include "Vertex.h" #include "MeshElement.h" template<unsigned int MeshDimension, typename IndexType, typename Real> class MeshWriter{ }; enum { hmmm protected: /** * @brief The MeshHash struct<HR> * A container to store cumulative data about a mesh * to recognize changes in the mesh. * It uses number of elements and sum of vertices coordinations. */ struct MeshHash { IndexType numberOfElements = 0; Real totalVert = 0; bool operator== (const MeshHash& rhs) const { return numberOfElements == rhs.numberOfElements && totalVert == rhs.totalVert; } bool operator!= (const MeshHash& rhs) const { return !((*this)==rhs); } }; template <typename IndexType, typename Real> class MeshWriter<2, IndexType, Real> { private: template<unsigned int Dim, typename Dummy = void> struct sumOfMeshElements{ template<unsigned int ...Reserve> static IndexType sum(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ return mesh.template getElements<Dim>().size() + sumOfMeshElements<Dim - 1>::sum(mesh); } }; template<typename Dummy> struct sumOfMeshElements<0, Dummy>{ template<unsigned int ...Reserve> static IndexType sum(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ return mesh.template getElements<0>().size(); } }; public: using type = MeshNativeType<2>; using type = MeshNativeType<MeshDimension>; template<unsigned int ...Reserve> static MeshHash computeHash(MeshElements<MeshDimension, IndexType, Real, Reserve...>& mesh){ MeshHash res; // A vector of ones that simplifies the sum of coordinates Vertex<MeshDimension, Real> ones; for(unsigned int dim = 0; dim < MeshDimension; dim++){ ones[dim] = 1.0; } for(IndexType i = 0; i < mesh.getVertices().size(); i++) { res.totalVert += mesh.getVertices().at(i) * ones; } res.numberOfElements = sumOfMeshElements<MeshDimension>::sum(mesh); return res; } }; #endif // MESHWRITER_H