+ //================================================================================
+ /*!
+ * \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 Return true if two adjacent pyramids are too close one to another
+ * so that a tetrahedron to built between them would 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->getshapeId() != nApexJ->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_TNodeXYZ( nApexI );
+ gp_XYZ apexJ = SMESH_TNodeXYZ( nApexJ );
+ gp_XYZ base1 = SMESH_TNodeXYZ( baseNodes[0]);
+ gp_XYZ base2 = SMESH_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. * M_PI / 180. );
+
+ // Check if pyramids collide
+ if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 ))
+ {
+ // find out if nI points outside of PrmI or inside
+ int dInd = baseNodesIndI[1] - baseNodesIndI[0];
+ bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+
+ // find out sign of projection of baI to nJ
+ 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;
+ dInd = baseNodesIndI[1] - baseNodesIndI[0];
+ bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+ dInd = baseNodesIndJ[1] - baseNodesIndJ[0];
+ bool isOutJ = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0;
+ if ( isOutJ == isOutI )
+ return false; // other domain
+
+ // direct both normals outside pyramid
+ ( isOutI ? nJ : nI ).Reverse();
+
+ // check absence of a face separating domains between pyramids
+ TIDSortedElemSet emptySet, avoidSet;
+ int i1, i2;
+ while ( const SMDS_MeshElement* f =
+ SMESH_MeshAlgos::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 );
+
+ if ( otherFaceNode == nApexI || otherFaceNode == nApexJ )
+ continue; // f is a temporary triangle
+
+ // 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 triangle normals (nI and nJ)
+ gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode ));
+ if ( nI * baOFN > 0 && nJ * baOFN > 0 &&
+ baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212
+ {
+ tooClose = false; // f is between pyramids
+ break;
+ }
+ }
+ }
+
+ return tooClose;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Move medium nodes of merged quadratic pyramids
+ */
+ //================================================================================
+
+ void UpdateQuadraticPyramids(const set<const SMDS_MeshNode*>& commonApex,
+ SMESHDS_Mesh* meshDS)
+ {
+ typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
+ TStdElemIterator itEnd;
+
+ // shift of node index to get medium nodes between the 4 base nodes and the apex
+ const int base2MediumShift = 9;
+
+ set<const SMDS_MeshNode*>::const_iterator nIt = commonApex.begin();
+ for ( ; nIt != commonApex.end(); ++nIt )
+ {
+ SMESH_TNodeXYZ apex( *nIt );
+
+ vector< const SMDS_MeshElement* > pyrams // pyramids sharing the apex node
+ ( TStdElemIterator( apex._node->GetInverseElementIterator( SMDSAbs_Volume )), itEnd );
+
+ // Select medium nodes to keep and medium nodes to remove
+
+ typedef map < const SMDS_MeshNode*, const SMDS_MeshNode*, TIDCompare > TN2NMap;
+ TN2NMap base2medium; // to keep
+ vector< const SMDS_MeshNode* > nodesToRemove;
+
+ for ( unsigned i = 0; i < pyrams.size(); ++i )
+ for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
+ {
+ SMESH_TNodeXYZ base = pyrams[i]->GetNode( baseIndex );
+ const SMDS_MeshNode* medium = pyrams[i]->GetNode( baseIndex + base2MediumShift );
+ TN2NMap::iterator b2m = base2medium.insert( make_pair( base._node, medium )).first;
+ if ( b2m->second != medium )
+ {
+ nodesToRemove.push_back( medium );
+ }
+ else
+ {
+ // move the kept medium node
+ gp_XYZ newXYZ = 0.5 * ( apex + base );
+ meshDS->MoveNode( medium, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
+ }
+ }
+
+ // Within pyramids, replace nodes to remove by nodes to keep
+
+ for ( unsigned i = 0; i < pyrams.size(); ++i )
+ {
+ vector< const SMDS_MeshNode* > nodes( pyrams[i]->begin_nodes(),
+ pyrams[i]->end_nodes() );
+ for ( int baseIndex = 0; baseIndex < PYRAM_APEX; ++baseIndex )
+ {
+ const SMDS_MeshNode* base = pyrams[i]->GetNode( baseIndex );
+ nodes[ baseIndex + base2MediumShift ] = base2medium[ base ];
+ }
+ meshDS->ChangeElementNodes( pyrams[i], &nodes[0], nodes.size());
+ }
+
+ // Remove the replaced nodes
+
+ if ( !nodesToRemove.empty() )
+ {
+ SMESHDS_SubMesh * sm = meshDS->MeshElements( nodesToRemove[0]->getshapeId() );
+ for ( unsigned i = 0; i < nodesToRemove.size(); ++i )
+ meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
+ }
+ }
+ }
+
+ //================================================================================
+ /*!
+ * \brief Store an error about overlapping faces
+ */
+ //================================================================================
+
+ bool overlapError( SMESH_Mesh& mesh,
+ const SMDS_MeshElement* face1,
+ const SMDS_MeshElement* face2,
+ const TopoDS_Shape& shape = TopoDS_Shape())
+ {
+ if ( !face1 || !face2 ) return false;
+
+ SMESH_Comment msg;
+ msg << "face " << face1->GetID() << " overlaps face " << face2->GetID();
+
+ SMESH_subMesh * sm = 0;
+ if ( shape.IsNull() )
+ {
+ sm = mesh.GetSubMesh( mesh.GetShapeToMesh() );
+ }
+ else if ( shape.ShapeType() >= TopAbs_SOLID )
+ {
+ sm = mesh.GetSubMesh( shape );
+ }
+ else
+ {
+ TopoDS_Iterator it ( shape );
+ if ( it.More() )
+ sm = mesh.GetSubMesh( it.Value() );
+ }
+ if ( sm )
+ {
+ SMESH_ComputeErrorPtr& err = sm->GetComputeError();
+ if ( !err || err->IsOK() )
+ {
+ SMESH_BadInputElements* badElems =
+ new SMESH_BadInputElements( mesh.GetMeshDS(),COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() );
+ badElems->add( face1 );
+ badElems->add( face2 );
+ err.reset( badElems );
+ }
+ }
+
+ return false; // == "algo fails"
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them
+ */
+//================================================================================
+
+void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* PrmI,
+ const SMDS_MeshElement* PrmJ,
+ set<const SMDS_MeshNode*> & nodesToMove)
+{
+ // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " "
+ // << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl;
+ const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
+ //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
+ SMESH_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_TNodeXYZ Pi( CommonNode );
+ gp_XYZ Pnew = /*( nbI*Pi + nbJ*Pj ) / (nbI+nbJ);*/ 0.5 * ( Pi + Pj );
+ CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() );
+
+ nodesToMove.insert( CommonNode );
+ nodesToMove.erase ( Nrem );
+
+ typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
+ TStdElemIterator itEnd;
+
+ // find and remove coincided faces of merged pyramids
+ vector< const SMDS_MeshElement* > inverseElems
+ // copy inverse elements to avoid iteration on changing container
+ ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd);
+ for ( size_t i = 0; i < inverseElems.size(); ++i )
+ {
+ const SMDS_MeshElement* FI = inverseElems[i];
+ 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 )
+ {
+ removeTmpElement( FI );
+ removeTmpElement( FJEqual );
+ myRemovedTrias.insert( FI );
+ myRemovedTrias.insert( FJEqual );
+ }
+ }
+
+ // set the common apex node to pyramids and triangles merged with J
+ vector< const SMDS_MeshNode* > nodes;
+ inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd );
+ for ( size_t i = 0; i < inverseElems.size(); ++i )
+ {
+ const SMDS_MeshElement* elem = inverseElems[i];
+ nodes.assign( elem->begin_nodes(), elem->end_nodes() );
+ nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
+ GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
+ }
+ ASSERT( Nrem->NbInverseElements() == 0 );
+ GetMeshDS()->RemoveFreeNode( Nrem,
+ GetMeshDS()->MeshElements( Nrem->getshapeId()),
+ /*fromGroups=*/false);
+}
+
+//================================================================================
+/*!
+ * \brief Merges adjacent pyramids
+ */
+//================================================================================
+
+void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI,
+ set<const SMDS_MeshNode*>& nodesToMove,
+ const bool isRecursion)
+{
+ TIDSortedElemSet adjacentPyrams;
+ bool mergedPyrams = false;
+ 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 == PrmI || PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second )
+ continue;
+ if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
+ {
+ MergePiramids( PrmI, PrmJ, nodesToMove );
+ mergedPyrams = true;
+ // container of inverse elements can change
+ // vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); -- iterator re-implemented
+ }
+ }
+ }
+ if ( mergedPyrams && !isRecursion )
+ {
+ TIDSortedElemSet::iterator prm;
+ for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
+ MergeAdjacent( *prm, nodesToMove, true );
+ }