From 7c66430d0e8a2bca47086d0af25b0bb90310455e Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 10 Nov 2008 07:39:59 +0000 Subject: [PATCH] dual mesh --- src/INTERP_KERNEL/DualMESH.cxx | 381 +++++++++++++++++++++++++++++++++ src/INTERP_KERNEL/DualMESH.hxx | 30 +++ src/INTERP_KERNEL/Makefile.am | 5 +- 3 files changed, 414 insertions(+), 2 deletions(-) create mode 100644 src/INTERP_KERNEL/DualMESH.cxx create mode 100644 src/INTERP_KERNEL/DualMESH.hxx diff --git a/src/INTERP_KERNEL/DualMESH.cxx b/src/INTERP_KERNEL/DualMESH.cxx new file mode 100644 index 000000000..87d52f574 --- /dev/null +++ b/src/INTERP_KERNEL/DualMESH.cxx @@ -0,0 +1,381 @@ +// File : DualMESH.cxx +// Created : Thu Oct 30 11:35:33 2008 +// Author : Edward AGAPOV (eap) + +#include "DualMESH.hxx" +#include "MEDMEM_Exception.hxx" + +using namespace MEDMEM; +using namespace MED_EN; +using namespace std; + + +namespace INTERP_KERNEL +{ +// ================================================================================ +/*! + * \brief Link of two nodes, stores middle node id, + * is used to find links of same nodes + */ +struct TLink: public pair< int, int > +{ + int _middleNode; + bool _reverse; + TLink():pair(),_middleNode(0) {} + TLink(int node1, int node2, int middleNode=0): pair(),_middleNode(middleNode) { + if ( node1 < node2 ) { first = node1; second = node2; _reverse = false; } + else { first = node2; second = node1; _reverse = true; } + } + int firstNode() const { return _reverse ? second : first; } + int secondNode() const { return _reverse ? first : second; } + int middleNode() const { return _middleNode; } + int otherNode(int n) const { return n==first ? second : n==second ? first : 0; } + //bool operator < (const TLink& l) const { return pair::operator<(l); } +}; + +// ================================================================================ +/*! + * \brief Triangle face, stores id of bary centre node, + * is used to find triangles of same nodes + */ +struct TTriaFace: public vector +{ + int _baryNode; + TTriaFace(){} + TTriaFace(int n1, int n2, int n3, int baryNode):vector(3),_baryNode(baryNode) + { at(0)=n1; at(1)=n2, at(2)=n3; std::sort(begin(), end()); } +}; + +// ================================================================================ +/*! + * \brief Creates a dual mesh from the given one + * + * The dual mesh consists of dual cells corresponding to nodes of the input mesh. + * The dual cell is a holygon in 2D space and a polyhedron in 3D space. The dual + * cell is bound by edges (2D) or triangles (3D) connecting barycentres of cells + * around the input node with middles of edges ending at the input node. + */ +DualMESH::DualMESH( MESH & mesh ):MESH() +{ + _name = "Dual mesh of " + mesh.getName(); + _spaceDimension = mesh.getSpaceDimension(); + _meshDimension = mesh.getMeshDimension(); + + const medEntityMesh entity = MED_CELL; + const medGeometryElement geom = (_spaceDimension == 2 ) ? MED_TRIA3 : MED_TETRA4; + const int nbNodesPerCell = geom % 100; + + // check type of elements in the input mesh + if ( mesh.getNumberOfTypesWithPoly(entity) != 1 || + mesh.getTypesWithPoly(entity)[0] != geom ) + throw MEDEXCEPTION("DualMESH() is only for tetrahedrons and triangles in 2D space!"); + + int nbInputCells = mesh.getNumberOfElements(entity, geom); + int nbDualCells = mesh.getNumberOfNodes(); + int nbCellBary = nbInputCells; + + // Compute bary centres of input cells and reverse nodal connectivity + + vector baryCenters( nbCellBary * _spaceDimension, 0.0 ); + vector< list< int > > inputCellsByNode ( mesh.getNumberOfNodes() ); + + const double * inputCoord = mesh.getCoordinates(MED_FULL_INTERLACE); + const int * inputConn = mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL, entity, geom); + + int iCell, iNode; + const int * cellConn = inputConn; + double * cellBary = & baryCenters[0]; + for ( iCell = 0; iCell < nbInputCells; ++iCell, cellBary += _spaceDimension ) { + for ( iNode = 0; iNode < nbNodesPerCell; ++iNode, ++cellConn ) { + int node = *cellConn - 1; + inputCellsByNode[ node ].push_back( iCell ); + const double * nodeCoord = inputCoord + node * _spaceDimension; + for ( int j = 0; j < _spaceDimension; ++j ) + cellBary[ j ] += nodeCoord[ j ]; + } + for ( int j = 0; j < _spaceDimension; ++j ) + cellBary[ j ] /= nbNodesPerCell; + } + + // Make dual nodal connectivity. + // nodes of dual mesh are: + // 1. nodes at bary centres of input cells + // 2. nodes at middles of input edges + // 3. nodes at bary centres of side faces of tetrahedrons (3D only) + // 4. boundary nodes + + _connectivity = new CONNECTIVITY( /*numberOfTypes = */0 ); + _connectivity->setEntityDimension( _spaceDimension ); + + list< pair< int, int > > bndNodes; // pairs of input and dual mesh boundary nodes + int nbBndNodes = 0; // counter of bndNodes + set< TLink > edges; // links between nodes storing id of middle node + set< TTriaFace > triangles; // triangles storing bary centre node id (3D only) + typedef set::iterator TEdgeIt; + typedef set::iterator TFaceIt; + + // -------------------------------------------------------------------------------- + if ( geom == MED_TRIA3 ) + { + list < int > dualConn; + vector< int > dualIndex; + dualIndex.reserve( nbDualCells + 1 ); + dualIndex.push_back( 1 ); + + for ( iNode = 1; iNode <= nbDualCells; ++iNode ) + { + // to make connectivity of poly around an input node, use map of middle node + // of edge to TLink between the middle node and a neighbor middle node. + // Each TLink represents two edges of a polygon: middle1-bary and bary-middle2 + map< int, TLink > middle2BaryLink; + typedef map< int, TLink >::iterator TMid2LinkIt; + // loop on triangles around iNode + const list & inputTria = inputCellsByNode[ iNode-1 ]; + list::const_iterator tri = inputTria.begin(), triEnd = inputTria.end(); + for ( ; tri != triEnd; ++tri ) + { + const int* triaConn = inputConn + *tri * nbNodesPerCell; // triangle connectivity + // set two other nodes of a triangle in a right order to + // correspond to connectivity convention + int node[2]; + if ( triaConn[0] == iNode ) { + node[0] = triaConn[1]; node[1] = triaConn[2]; + } + else if ( triaConn[1] == iNode ) { + node[0] = triaConn[2]; node[1] = triaConn[0]; + } + else { + node[0] = triaConn[0]; node[1] = triaConn[1]; + } + // get middle nodes on edges iNode-node[0] and iNode-node[1] + int middle[2]; + for ( int i = 0; i < 2; ++i ) { + int nextMiddle = nbCellBary + edges.size() + nbBndNodes + 1; // next free node id + TLink edge( iNode, node[i], nextMiddle ); + TEdgeIt edgeIt = edges.insert( edge ).first; + middle[i] = edgeIt->middleNode(); + } + // fill the map with link between middle[0] and middle[1] storing bary center of a triangle + int baryNode = *tri+1; + TLink link( middle[0], middle[1], baryNode ); + middle2BaryLink[ middle[0] ] = link; + } + + // Form polygon connectivity + + list< int > polyConn; + polyConn.push_back( middle2BaryLink.begin()->first ); + while (1) { + int lastNode = polyConn.back(); + TMid2LinkIt mid2LinkIt = middle2BaryLink.find( lastNode ); + if ( mid2LinkIt == middle2BaryLink.end() ) + break; + const TLink& baryLink = mid2LinkIt->second; + polyConn.push_back( baryLink.middleNode() ); + polyConn.push_back( baryLink.otherNode( lastNode )); + middle2BaryLink.erase( mid2LinkIt ); + } + if ( polyConn.front() != polyConn.back() ) { + // iNode is on mesh boundary. + int nextBnd = nbCellBary + edges.size() + nbBndNodes + 1; // next free node id + bndNodes.push_back( make_pair( iNode, nextBnd )); + nbBndNodes++; + // Prepend segments remaining in middle2BaryLink + TMid2LinkIt mid2LinkIt, mid2LinkEnd; + while ( !middle2BaryLink.empty() ) { + int firstNode = polyConn.front(); + mid2LinkIt = middle2BaryLink.begin(); mid2LinkEnd = middle2BaryLink.end(); + for ( ; mid2LinkIt != mid2LinkEnd; ++mid2LinkIt ) { + const TLink& baryLink = mid2LinkIt->second; + if ( int n = baryLink.otherNode( firstNode )) { + polyConn.push_front( baryLink.middleNode() ); + polyConn.push_front( n ); + middle2BaryLink.erase( mid2LinkIt ); + break; + } + } + } + polyConn.push_front( nextBnd ); + } + else { + polyConn.pop_front(); + } + // add polyConn to global connectivity + dualIndex.push_back( dualIndex.back() + polyConn.size() ); + dualConn.splice( dualConn.end(), polyConn ); + + } // loop on input nodes + + vector< int > conn ( dualConn.begin(), dualConn.end() ); + _connectivity->setPolygonsConnectivity(MED_NODAL, entity, + &conn[0], &dualIndex[0], + conn.size(), nbDualCells); + + } // -------------------------------------------------------------------------------- + + else // here geom == MED_TETRA4 + { + list < int > dualConn, faceIndex; + vector< int > dualIndex; + dualIndex.reserve( nbDualCells + 1 ); + dualIndex.push_back( 1 ); + faceIndex.push_back( 1 ); + + for ( iNode = 1; iNode <= nbDualCells; ++iNode ) + { + // use map of middle-midle links to nb of links for detecting boundary faces + map< TLink, int > middleLinkCounters; + typedef map< TLink, int >::iterator TLinkNbIt; + int nbFaces = 0; + // loop on tetras around iNode + const list & inputTetra = inputCellsByNode[ iNode-1 ]; + list::const_iterator tet = inputTetra.begin(), triEnd = inputTetra.end(); + for ( ; tet != triEnd; ++tet ) + { + const int* tetConn = inputConn + *tet * nbNodesPerCell; // tetrahedron connectivity + + // set tree nodes of a tetrahedron other than iNode in the + // order to form a triangle with a normal pointing inside polyhedron + int node[4]; + if ( tetConn[0] == iNode ) { + node[0] = tetConn[1]; node[1] = tetConn[2]; node[2] = tetConn[3]; + } + else if ( tetConn[1] == iNode ) { + node[0] = tetConn[2]; node[1] = tetConn[0]; node[2] = tetConn[3]; + } + else if ( tetConn[2] == iNode ) { + node[0] = tetConn[3]; node[1] = tetConn[0]; node[2] = tetConn[1]; + } + else { + node[0] = tetConn[0]; node[1] = tetConn[2]; node[2] = tetConn[1]; + } + node[3] = node[0]; + + // get 3 middle nodes on edges iNode-node[i] + // and 3 bary centres of tetrahedron faces sharing iNode. + int middle[4], triaBary[3], i; + for ( i = 0; i < 3; ++i ) { + int nextMiddle = nbCellBary + edges.size() + triangles.size() + nbBndNodes + 1;//free id + TLink edge( iNode, node[i], nextMiddle ); + // store or find edge that will be used later to compute coordinates of middle node + TEdgeIt edgeIt = edges.insert( edge ).first; + middle[i] = edgeIt->middleNode(); + // face bary centre + int nextBary = nbCellBary + edges.size() + triangles.size() + nbBndNodes + 1; // free id + TTriaFace tria( iNode, node[i], node[i+1], nextBary ); + // store or find triangle that will be used later to compute coordinates of its BC + TFaceIt faceIt = triangles.insert( tria ).first; + triaBary[i] = faceIt->_baryNode; + } + middle[3] = middle[0]; + + // Add triangle faces to polyhedron connectivity and fill middleLinkCounters + int tetBary = *tet+1; + for ( i = 0; i < 3; ++i ) + { + dualConn.push_back( middle[i] ); + dualConn.push_back( tetBary ); + dualConn.push_back( triaBary[i] ); + faceIndex.push_back( faceIndex.back() + 3 ); + nbFaces++; + + dualConn.push_back( middle[i+1] ); + dualConn.push_back( triaBary[i] ); + dualConn.push_back( tetBary ); + faceIndex.push_back( faceIndex.back() + 3 ); + nbFaces++; + + TLink link ( middle[i], middle[i+1], triaBary[i] ); + TLinkNbIt linkNb = middleLinkCounters.insert( make_pair( link, 0 )).first; + linkNb->second++; + } + } // loop on tetras around iNode + + // Add faces on mesh boundary + + TLinkNbIt linkNb = middleLinkCounters.begin(), linkNbEnd = middleLinkCounters.end(); + int nextBnd = 0; + for ( ; linkNb != linkNbEnd; ++linkNb ) { + if ( linkNb->second < 2 ) + { + const TLink& link = linkNb->first; + int middle1 = link.firstNode(); + int middle2 = link.secondNode(); + int triBary = link.middleNode(); + if ( !nextBnd ) { + nextBnd = nbCellBary + edges.size() + triangles.size() + nbBndNodes + 1; // free id + bndNodes.push_back( make_pair( iNode, nextBnd )); + nbBndNodes++; + } + dualConn.push_back( middle1 ); + dualConn.push_back( triBary ); + dualConn.push_back( middle2 ); + dualConn.push_back( nextBnd ); + + faceIndex.push_back( faceIndex.back() + 4 ); + nbFaces++; + } + } + dualIndex.push_back( dualIndex.back() + nbFaces ); + + } // loop on input nodes + + vector< int > polyConnectivity( dualConn.begin(), dualConn.end() ); + vector< int > polyFaceIndex( faceIndex.begin(), faceIndex.end() ); + + _connectivity->setPolyhedronConnectivity(MED_NODAL, + & polyConnectivity[0], & dualIndex[0], + polyConnectivity.size(), dualIndex.size()-1, + & polyFaceIndex[0], polyFaceIndex.size()-1); + } + // -------------------------------------------------------------------------------- + + // Create coordinates + + _numberOfNodes = nbCellBary + edges.size() + triangles.size() + bndNodes.size(); + _coordinate = new COORDINATE( _spaceDimension, _numberOfNodes, MED_FULL_INTERLACE ); + + _coordinate->setCoordinatesNames( mesh.getCoordinatesNames() ); + _coordinate->setCoordinatesUnits( mesh.getCoordinatesUnits() ); + _coordinate->setCoordinatesSystem( mesh.getCoordinatesSystem() ); + + vector & coord = baryCenters; + coord.resize( _spaceDimension * _numberOfNodes ); + + // set coords of nodes at middles of input edges + TEdgeIt edge = edges.begin(), edgeEnd = edges.end(); + for ( ; edge != edgeEnd; ++edge ) + { + int n1 = ( edge->firstNode() - 1 ) * _spaceDimension; + int n2 = ( edge->secondNode() - 1 ) * _spaceDimension; + int mi = ( edge->middleNode() - 1 ) * _spaceDimension; + for ( int j = 0; j < _spaceDimension; ++j ) + coord[ mi+j ] = 0.5 * ( inputCoord[ n1+j ] + inputCoord[ n2+j ]); + } + + // set coords of nodes at bary centres of side faces of tetras + TFaceIt tri = triangles.begin(), triEnd = triangles.end(); + for ( ; tri != triEnd; ++tri ) + { + const TTriaFace & triangle = *tri; + int n0 = ( triangle[0] - 1 ) * _spaceDimension; + int n1 = ( triangle[1] - 1 ) * _spaceDimension; + int n2 = ( triangle[2] - 1 ) * _spaceDimension; + int bc = ( triangle._baryNode - 1 ) * _spaceDimension; + for ( int j = 0; j < _spaceDimension; ++j ) + coord[ bc+j ] = ( inputCoord[ n0+j ] + inputCoord[ n1+j ] + inputCoord[ n2+j ]) / 3.; + } + + // set coords of boundary nodes + list< pair< int, int > >::iterator nodeIt = bndNodes.begin(), nodeEnd = bndNodes.end(); + for ( ; nodeIt != nodeEnd; ++nodeIt ) + { + int inputNode = ( nodeIt->first - 1 ) * _spaceDimension; + int dualNode = ( nodeIt->second - 1 ) * _spaceDimension; + for ( int j = 0; j < _spaceDimension; ++j ) + coord[ dualNode+j ] = inputCoord[ inputNode+j ]; + } + _coordinate->setCoordinates( MED_FULL_INTERLACE, & coord[0] ); +} + +} // namespace ParaMEDMEM diff --git a/src/INTERP_KERNEL/DualMESH.hxx b/src/INTERP_KERNEL/DualMESH.hxx new file mode 100644 index 000000000..0c681dc8a --- /dev/null +++ b/src/INTERP_KERNEL/DualMESH.hxx @@ -0,0 +1,30 @@ +// File : DualMESH.hxx +// Created : Thu Oct 30 11:19:22 2008 +// Author : Edward AGAPOV (eap) + + +#ifndef DualMESH_HeaderFile +#define DualMESH_HeaderFile + +#include "MEDMEM_Mesh.hxx" + +namespace INTERP_KERNEL +{ + +/*! + * \brief Creates a dual mesh from the given one + * + * The dual mesh consists of dual cells corresponding to nodes of the input mesh. + * The dual cell is a holygon in 2D space and a polyhedron in 3D space. The dual + * cell is bound by edges (2D) or triangles (3D) connecting barycentres of cells + * around the input node with middles of edges ending at the input node. + */ +class DualMESH: public MEDMEM::MESH +{ +public: + DualMESH( MEDMEM::MESH & mesh ); +}; + +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 8f766a650..ee0302135 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -56,7 +56,7 @@ ConvexIntersector.txx IntersectorHexa.txx PointLocatorAlgos Geometric2DIntersector.txx IntersectorTetra.txx PolygonAlgorithms.txx\ Interpolation2D.txx MEDNormalizedUnstructuredMesh.txx TriangulationIntersector.txx\ Interpolation3DSurf.txx MeshElement.txx VTKNormalizedUnstructuredMesh.txx\ -Interpolation3D.txx MeshRegion.txx +Interpolation3D.txx MeshRegion.txx DualMESH.hxx # Libraries targets @@ -68,7 +68,8 @@ dist_libinterpkernel_la_SOURCES = \ TetraAffineTransform.cxx\ CellModel.cxx \ Remapper.cxx \ - PointLocator.cxx + PointLocator.cxx \ + DualMESH.cxx libinterpkernel_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) ../MEDWrapper/V2_1/Core/libmed_V2_1.la $(STDLIB) ../MEDMEM/libmedmem.la -lutil -lm -- 2.39.2