+namespace
+{
+ const int Q2TAbs_TmpTriangle = SMDSEntity_Last + 1;
+
+ //================================================================================
+ /*!
+ * \brief Temporary face. It's key feature is that it adds itself to an apex node
+ * as inverse element, so that tmp triangles of a piramid can be easily found
+ */
+ //================================================================================
+
+ class STDMESHERS_EXPORT Q2TAdaptor_Triangle : public SMDS_MeshFace
+ {
+ const SMDS_MeshNode* _nodes[3];
+ public:
+ Q2TAdaptor_Triangle(const SMDS_MeshNode* apexNode,
+ const SMDS_MeshNode* node2,
+ const SMDS_MeshNode* node3)
+ {
+ _nodes[0]=0; ChangeApex(apexNode);
+ _nodes[1]=node2;
+ _nodes[2]=node3;
+ }
+ ~Q2TAdaptor_Triangle() { MarkAsRemoved(); }
+ void ChangeApex(const SMDS_MeshNode* node)
+ {
+ MarkAsRemoved();
+ _nodes[0]=node;
+ const_cast<SMDS_MeshNode*>(node)->AddInverseElement(this);
+ }
+ void MarkAsRemoved()
+ {
+ if ( _nodes[0] )
+ const_cast<SMDS_MeshNode*>(_nodes[0])->RemoveInverseElement(this), _nodes[0] = 0;
+ }
+ bool IsRemoved() const { return !_nodes[0]; }
+ virtual int NbNodes() const { return 3; }
+ virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nodes[ind]; }
+ virtual SMDSAbs_EntityType GetEntityType() const
+ {
+ return SMDSAbs_EntityType( Q2TAbs_TmpTriangle );
+ }
+ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
+ {
+ if ( type == SMDSAbs_Node )
+ return SMDS_ElemIteratorPtr( new SMDS_NodeArrayElemIterator( _nodes, _nodes+3 ));
+ throw SALOME_Exception(LOCALIZED("Not implemented"));
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return true if two nodes of triangles are equal
+ */
+ //================================================================================
+
+ bool EqualTriangles(const SMDS_MeshElement* F1,const SMDS_MeshElement* F2)
+ {
+ return
+ ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) ||
+ ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
+ */
+ //================================================================================
+
+ void MergePiramids( const SMDS_MeshElement* PrmI,
+ const SMDS_MeshElement* PrmJ,
+ SMESHDS_Mesh* meshDS,
+ set<const SMDS_MeshNode*> & nodesToMove)
+ {
+ const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
+ int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
+ SMESH_MeshEditor::TNodeXYZ Pj( Nrem );
+
+ // an apex node to make common to all merged pyramids
+ SMDS_MeshNode* CommonNode = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4));
+ if ( CommonNode == Nrem ) return; // already merged
+ int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume );
+ SMESH_MeshEditor::TNodeXYZ Pi( CommonNode );
+ gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);
+ CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
+
+ nodesToMove.insert( CommonNode );
+ nodesToMove.erase ( Nrem );
+
+ // find and remove coincided faces of merged pyramids
+ SMDS_ElemIteratorPtr triItI = CommonNode->GetInverseElementIterator(SMDSAbs_Face);
+ while ( triItI->more() )
+ {
+ const SMDS_MeshElement* FI = triItI->next();
+ const SMDS_MeshElement* FJEqual = 0;
+ SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face);
+ while ( !FJEqual && triItJ->more() )
+ {
+ const SMDS_MeshElement* FJ = triItJ->next();
+ if ( EqualTriangles( FJ, FI ))
+ FJEqual = FJ;
+ }
+ if ( FJEqual )
+ {
+ ((Q2TAdaptor_Triangle*) FI)->MarkAsRemoved();
+ ((Q2TAdaptor_Triangle*) FJEqual)->MarkAsRemoved();
+ }
+ }
+
+ // set the common apex node to pyramids and triangles merged with J
+ SMDS_ElemIteratorPtr itJ = Nrem->GetInverseElementIterator();
+ while ( itJ->more() )
+ {
+ const SMDS_MeshElement* elem = itJ->next();
+ if ( elem->GetType() == SMDSAbs_Volume ) // pyramid
+ {
+ vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() );
+ nodes[4] = CommonNode;
+ meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size());
+ }
+ else if ( elem->GetEntityType() == Q2TAbs_TmpTriangle ) // tmp triangle
+ {
+ ((Q2TAdaptor_Triangle*) elem )->ChangeApex( CommonNode );
+ }
+ }
+ ASSERT( Nrem->NbInverseElements() == 0 );
+ meshDS->RemoveFreeNode( Nrem,
+ meshDS->MeshElements( Nrem->GetPosition()->GetShapeId()),
+ /*fromGroups=*/false);
+ }
+
+ //================================================================================
+ /*!
+ * \brief Return true if two adjacent pyramids are too close one to another
+ * so that a tetrahedron to built between them whoul have too poor quality
+ */
+ //================================================================================
+
+ bool TooCloseAdjacent( const SMDS_MeshElement* PrmI,
+ const SMDS_MeshElement* PrmJ,
+ const bool hasShape)
+ {
+ const SMDS_MeshNode* nApexI = PrmI->GetNode(4);
+ const SMDS_MeshNode* nApexJ = PrmJ->GetNode(4);
+ if ( nApexI == nApexJ ||
+ nApexI->GetPosition()->GetShapeId() != nApexJ->GetPosition()->GetShapeId() )
+ return false;
+
+ // Find two common base nodes and their indices within PrmI and PrmJ
+ const SMDS_MeshNode* baseNodes[2] = { 0,0 };
+ int baseNodesIndI[2], baseNodesIndJ[2];
+ for ( int i = 0; i < 4 ; ++i )
+ {
+ int j = PrmJ->GetNodeIndex( PrmI->GetNode(i));
+ if ( j >= 0 )
+ {
+ int ind = baseNodes[0] ? 1:0;
+ if ( baseNodes[ ind ])
+ return false; // pyramids with a common base face
+ baseNodes [ ind ] = PrmI->GetNode(i);
+ baseNodesIndI[ ind ] = i;
+ baseNodesIndJ[ ind ] = j;
+ }
+ }
+ if ( !baseNodes[1] ) return false; // not adjacent
+
+ // Get normals of triangles sharing baseNodes
+ gp_XYZ apexI = SMESH_MeshEditor::TNodeXYZ( nApexI );
+ gp_XYZ apexJ = SMESH_MeshEditor::TNodeXYZ( nApexJ );
+ gp_XYZ base1 = SMESH_MeshEditor::TNodeXYZ( baseNodes[0]);
+ gp_XYZ base2 = SMESH_MeshEditor::TNodeXYZ( baseNodes[1]);
+ gp_Vec baseVec( base1, base2 );
+ gp_Vec baI( base1, apexI );
+ gp_Vec baJ( base1, apexJ );
+ gp_Vec nI = baseVec.Crossed( baI );
+ gp_Vec nJ = baseVec.Crossed( baJ );
+
+ // Check angle between normals
+ double angle = nI.Angle( nJ );
+ bool tooClose = ( angle < 15 * PI180 );
+
+ // Check if pyramids collide
+ bool isOutI, isOutJ;
+ if ( !tooClose && baI * baJ > 0 )
+ {
+ // find out if nI points outside of PrmI or inside
+ int dInd = baseNodesIndI[1] - baseNodesIndI[0];
+ isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+
+ // find out sign of projection of nJ to baI
+ double proj = baI * nJ;
+
+ tooClose = isOutI ? proj > 0 : proj < 0;
+ }
+
+ // Check if PrmI and PrmJ are in same domain
+ if ( tooClose && !hasShape )
+ {
+ // check order of baseNodes within pyramids, it must be opposite
+ int dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
+ isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+ if ( isOutJ == isOutI )
+ return false; // other domain
+
+ // check absence of a face separating domains between pyramids
+ TIDSortedElemSet emptySet, avoidSet;
+ int i1, i2;
+ while ( const SMDS_MeshElement* f =
+ SMESH_MeshEditor::FindFaceInSet( baseNodes[0], baseNodes[1],
+ emptySet, avoidSet, &i1, &i2 ))
+ {
+ avoidSet.insert( f );
+
+ // face node other than baseNodes
+ int otherNodeInd = 0;
+ while ( otherNodeInd == i1 || otherNodeInd == i2 ) otherNodeInd++;
+ const SMDS_MeshNode* otherFaceNode = f->GetNode( otherNodeInd );
+
+ // check if f is a base face of either of pyramids
+ if ( f->NbCornerNodes() == 4 &&
+ ( PrmI->GetNodeIndex( otherFaceNode ) >= 0 ||
+ PrmJ->GetNodeIndex( otherFaceNode ) >= 0 ))
+ continue; // f is a base quadrangle
+
+ // check projections of face direction (baOFN) to triange normals (nI and nJ)
+ gp_Vec baOFN( base1, SMESH_MeshEditor::TNodeXYZ( otherFaceNode ));
+ ( isOutI ? nJ : nI ).Reverse();
+ if ( nI * baOFN > 0 && nJ * baOFN > 0 )
+ {
+ tooClose = false; // f is between pyramids
+ break;
+ }
+ }
+ }
+
+ return tooClose;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Merges adjacent pyramids
+ */
+ //================================================================================
+
+ void MergeAdjacent(const SMDS_MeshElement* PrmI,
+ SMESH_Mesh& mesh,
+ set<const SMDS_MeshNode*>& nodesToMove)
+ {
+ TIDSortedElemSet adjacentPyrams, mergedPyrams;
+ for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI
+ {
+ const SMDS_MeshNode* n = PrmI->GetNode(k);
+ SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
+ while ( vIt->more() )
+ {
+ const SMDS_MeshElement* PrmJ = vIt->next();
+ if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second )
+ continue;
+ if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, mesh.HasShapeToMesh() ))
+ {
+ MergePiramids( PrmI, PrmJ, mesh.GetMeshDS(), nodesToMove );
+ mergedPyrams.insert( PrmJ );
+ }
+ }
+ }
+ if ( !mergedPyrams.empty() )
+ {
+ TIDSortedElemSet::iterator prm;
+// for (prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm)
+// MergeAdjacent( *prm, mesh, nodesToMove );
+
+ for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
+ MergeAdjacent( *prm, mesh, nodesToMove );
+ }
+ }
+}
+