--- /dev/null
+// 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<int, int>(),_middleNode(0) {}
+ TLink(int node1, int node2, int middleNode=0): pair<int, int>(),_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>
+{
+ int _baryNode;
+ TTriaFace(){}
+ TTriaFace(int n1, int n2, int n3, int baryNode):vector<int>(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<double> 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<TLink>::iterator TEdgeIt;
+ typedef set<TTriaFace>::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<int> & inputTria = inputCellsByNode[ iNode-1 ];
+ list<int>::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<int> & inputTetra = inputCellsByNode[ iNode-1 ];
+ list<int>::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<double> & 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
