-// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
//
-// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-// SMESH SMESH : idl implementation based on 'SMESH' unit's classes
+// SMESH SMESH : idl implementation based on 'SMESH' unit's classes
// File : SMESH_MeshEditor.cxx
// Created : Mon Apr 12 16:10:22 2004
// Author : Edward AGAPOV (eap)
-//
+
#define CHRONODEF
#include "SMESH_MeshEditor.hxx"
#include "SMESH_OctreeNode.hxx"
#include "SMESH_subMesh.hxx"
+#include <Basics_OCCTVersion.hxx>
+
#include "utilities.h"
#include <BRepAdaptor_Surface.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
#include <BRep_Tool.hxx>
#include <ElCLib.hxx>
#include <gp_XY.hxx>
#include <gp_XYZ.hxx>
-#include <math.h>
+#include <cmath>
#include <map>
#include <set>
#include <numeric>
#include <limits>
+#include <algorithm>
+#include <sstream>
#define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
if ( aMesh->ShapeToMesh().IsNull() )
return 0;
- if ( theElem->GetType() == SMDSAbs_Node )
- {
- int aShapeID = theElem->getshapeId();
- if (aShapeID <= 0)
- return 0;
- else
- return aShapeID;
- }
-
- TopoDS_Shape aShape; // the shape a node is on
- SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
- while ( nodeIt->more() ) {
- const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
- int aShapeID = node->getshapeId();
- if (aShapeID > 0) {
- SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID );
- if ( sm ) {
- if ( sm->Contains( theElem ))
- return aShapeID;
- if ( aShape.IsNull() )
- aShape = aMesh->IndexToShape( aShapeID );
- }
- else {
- //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID );
+ int aShapeID = theElem->getshapeId();
+ if ( aShapeID < 1 )
+ return 0;
+
+ if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
+ if ( sm->Contains( theElem ))
+ return aShapeID;
+
+ if ( theElem->GetType() == SMDSAbs_Node ) {
+ MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
+ }
+ else {
+ MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
+ }
+
+ TopoDS_Shape aShape; // the shape a node of theElem is on
+ if ( theElem->GetType() != SMDSAbs_Node )
+ {
+ SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
+ while ( nodeIt->more() ) {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
+ if ((aShapeID = node->getshapeId()) > 0) {
+ if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
+ if ( sm->Contains( theElem ))
+ return aShapeID;
+ if ( aShape.IsNull() )
+ aShape = aMesh->IndexToShape( aShapeID );
+ }
}
}
}
// None of nodes is on a proper shape,
// find the shape among ancestors of aShape on which a node is
- if ( aShape.IsNull() ) {
- //MESSAGE ("::FindShape() - NONE node is on shape")
- return 0;
+ if ( !aShape.IsNull() ) {
+ TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
+ for ( ; ancIt.More(); ancIt.Next() ) {
+ SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
+ if ( sm && sm->Contains( theElem ))
+ return aMesh->ShapeToIndex( ancIt.Value() );
+ }
}
- TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
- for ( ; ancIt.More(); ancIt.Next() ) {
- SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
- if ( sm && sm->Contains( theElem ))
- return aMesh->ShapeToIndex( ancIt.Value() );
+ else
+ {
+ const map<int,SMESHDS_SubMesh*>& id2sm = GetMeshDS()->SubMeshes();
+ map<int,SMESHDS_SubMesh*>::const_iterator id_sm = id2sm.begin();
+ for ( ; id_sm != id2sm.end(); ++id_sm )
+ if ( id_sm->second->Contains( theElem ))
+ return id_sm->first;
}
//MESSAGE ("::FindShape() - SHAPE NOT FOUND")
aNodes[5] = nd2;
}
+//=======================================================================
+//function : edgeConnectivity
+//purpose : auxilary
+// return number of the edges connected with the theNode.
+// if theEdges has connections with the other type of the
+// elements, return -1
+//=======================================================================
+static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
+{
+ SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
+ int nb=0;
+ while(elemIt->more()) {
+ elemIt->next();
+ nb++;
+ }
+ return nb;
+}
+
+
//=======================================================================
//function : GetNodesFromTwoTria
//purpose : auxilary
myLastCreatedNodes.Append(newN);
// create a new element
- const SMDS_MeshNode* N[6];
if ( aBadRate1 <= aBadRate2 ) {
- N[0] = aNodes[0];
- N[1] = aNodes[1];
- N[2] = aNodes[2];
- N[3] = aNodes[4];
- N[4] = aNodes[5];
- N[5] = newN;
newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
aNodes[6], aNodes[7], newN );
newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
newN, aNodes[4], aNodes[5] );
}
else {
- N[0] = aNodes[1];
- N[1] = aNodes[2];
- N[2] = aNodes[3];
- N[3] = aNodes[5];
- N[4] = aNodes[6];
- N[5] = newN;
newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
aNodes[7], aNodes[4], newN );
newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
helper.SetIsQuadratic( false );
}
vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
+ helper.SetElementsOnShape( true );
if ( splitMethod._baryNode )
{
// make a node at barycenter
}
// make tetras
- helper.SetElementsOnShape( true );
vector<const SMDS_MeshElement* > tetras( splitMethod._nbTetra ); // splits of a volume
const int* tetConn = splitMethod._connectivity;
for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 )
helper.SetElementsOnShape( false );
vector< const SMDS_MeshElement* > triangles;
+ // find submesh to add new triangles in
+ if ( !fSubMesh || !fSubMesh->Contains( face ))
+ {
+ int shapeID = FindShape( face );
+ fSubMesh = GetMeshDS()->MeshElements( shapeID );
+ }
map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
if ( iF_n != splitMethod._faceBaryNode.end() )
{
if ( !volTool.IsFaceExternal( iF ))
swap( n2, n3 );
triangles.push_back( helper.AddFace( n1,n2,n3 ));
+
+ if ( fSubMesh && n3->getshapeId() < 1 )
+ fSubMesh->AddNode( n3 );
}
}
else
volNodes[ facet->_n3 ]));
}
}
- // find submesh to add new triangles in
- if ( !fSubMesh || !fSubMesh->Contains( face ))
- {
- int shapeID = FindShape( face );
- fSubMesh = GetMeshDS()->MeshElements( shapeID );
- }
for ( int i = 0; i < triangles.size(); ++i )
{
if ( !triangles[i] ) continue;
// create a new element
const SMDS_MeshElement* newElem1 = 0;
const SMDS_MeshElement* newElem2 = 0;
- const SMDS_MeshNode* N[6];
if ( the13Diag ) {
- N[0] = aNodes[0];
- N[1] = aNodes[1];
- N[2] = aNodes[2];
- N[3] = aNodes[4];
- N[4] = aNodes[5];
- N[5] = newN;
newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
aNodes[6], aNodes[7], newN );
newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
newN, aNodes[4], aNodes[5] );
}
else {
- N[0] = aNodes[1];
- N[1] = aNodes[2];
- N[2] = aNodes[3];
- N[3] = aNodes[5];
- N[4] = aNodes[6];
- N[5] = newN;
newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
aNodes[7], aNodes[4], newN );
newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
if ( projector.IsDone() ) {
double u, v, minVal = DBL_MAX;
for ( int i = projector.NbExt(); i > 0; i-- )
+#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
+ if ( projector.SquareDistance( i ) < minVal ) {
+ minVal = projector.SquareDistance( i );
+#else
if ( projector.Value( i ) < minVal ) {
minVal = projector.Value( i );
+#endif
projector.Point( i ).Parameter( u, v );
}
result.SetCoord( u, v );
Handle(Geom_Surface) surface;
SMESHDS_SubMesh* faceSubMesh = 0;
TopoDS_Face face;
- double fToler2 = 0, vPeriod = 0., uPeriod = 0., f,l;
+ double fToler2 = 0, f,l;
double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
bool isUPeriodic = false, isVPeriodic = false;
if ( *fId ) {
fToler2 *= fToler2 * 10.;
isUPeriodic = surface->IsUPeriodic();
if ( isUPeriodic )
- vPeriod = surface->UPeriod();
+ surface->UPeriod();
isVPeriodic = surface->IsVPeriodic();
if ( isVPeriodic )
- uPeriod = surface->VPeriod();
+ surface->VPeriod();
surface->Bounds( u1, u2, v1, v2 );
}
// ---------------------------------------------------------
// fix nodes on mesh boundary
if ( checkBoundaryNodes ) {
- map< NLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
- map< NLink, int >::iterator link_nb;
+ map< SMESH_TLink, int > linkNbMap; // how many times a link encounters in elemsOnFace
+ map< SMESH_TLink, int >::iterator link_nb;
// put all elements links to linkNbMap
list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
const SMDS_MeshElement* elem = (*elemIt);
- int nbn = elem->NbNodes();
- if(elem->IsQuadratic())
- nbn = nbn/2;
+ int nbn = elem->NbCornerNodes();
// loop on elem links: insert them in linkNbMap
- const SMDS_MeshNode* curNode, *prevNode = elem->GetNodeWrap( nbn );
for ( int iN = 0; iN < nbn; ++iN ) {
- curNode = elem->GetNode( iN );
- NLink link;
- if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
- else link = make_pair( prevNode , curNode );
- prevNode = curNode;
- link_nb = linkNbMap.find( link );
- if ( link_nb == linkNbMap.end() )
- linkNbMap.insert( make_pair ( link, 1 ));
- else
- link_nb->second++;
+ const SMDS_MeshNode* n1 = elem->GetNode( iN );
+ const SMDS_MeshNode* n2 = elem->GetNode(( iN+1 ) % nbn);
+ SMESH_TLink link( n1, n2 );
+ link_nb = linkNbMap.insert( make_pair( link, 0 )).first;
+ link_nb->second++;
}
}
// remove nodes that are in links encountered only once from setMovableNodes
for ( link_nb = linkNbMap.begin(); link_nb != linkNbMap.end(); ++link_nb ) {
if ( link_nb->second == 1 ) {
- setMovableNodes.erase( link_nb->first.first );
- setMovableNodes.erase( link_nb->first.second );
+ setMovableNodes.erase( link_nb->first.node1() );
+ setMovableNodes.erase( link_nb->first.node2() );
}
}
}
for ( int iStep = 0; iStep < nbSteps; iStep++ ) {
vTool.Set( *v );
vTool.SetExternalNormal();
+ const int nextShift = vTool.IsForward() ? +1 : -1;
list< int >::iterator ind = freeInd.begin();
list< const SMDS_MeshElement* >::iterator srcEdge = srcEdges.begin();
for ( ; ind != freeInd.end(); ++ind, ++srcEdge ) // loop on free faces
{
const SMDS_MeshNode** nodes = vTool.GetFaceNodes( *ind );
int nbn = vTool.NbFaceNodes( *ind );
- switch ( nbn ) {
- case 3: { ///// triangle
- const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
- if ( !f )
- myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
- else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+ if ( ! (*v)->IsPoly() )
+ switch ( nbn ) {
+ case 3: { ///// triangle
+ const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ]);
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ]) + nextShift ))
{
- myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] ));
- aMesh->RemoveElement(f);
+ const SMDS_MeshNode* newOrder[3] = { nodes[ 1 - nextShift ],
+ nodes[ 1 ],
+ nodes[ 1 + nextShift ] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ] ));
}
- break;
- }
- case 4: { ///// quadrangle
- const SMDS_MeshFace * f = aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
- if ( !f )
- myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
- else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
+ break;
+ }
+ case 4: { ///// quadrangle
+ const SMDS_MeshFace * f =
+ aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ]);
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ]) + nextShift ))
{
- myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ] ));
- aMesh->RemoveElement(f);
- }
- break;
- }
- default:
- if( (*v)->IsQuadratic() ) {
- if(nbn==6) { /////// quadratic triangle
- const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
- nodes[1], nodes[3], nodes[5] );
- if ( !f ) {
- myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
- nodes[1], nodes[3], nodes[5]));
- }
- else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 )) {
- const SMDS_MeshNode** tmpnodes = new const SMDS_MeshNode*[6];
- tmpnodes[0] = nodes[0];
- tmpnodes[1] = nodes[2];
- tmpnodes[2] = nodes[4];
- tmpnodes[3] = nodes[1];
- tmpnodes[4] = nodes[3];
- tmpnodes[5] = nodes[5];
- myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4],
- nodes[1], nodes[3], nodes[5]));
- aMesh->RemoveElement(f);
- }
+ const SMDS_MeshNode* newOrder[4] = { nodes[ 0 ], nodes[ 2-nextShift ],
+ nodes[ 2 ], nodes[ 2+nextShift ] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ]));
}
- else { /////// quadratic quadrangle
- const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
- nodes[1], nodes[3], nodes[5], nodes[7] );
- if ( !f ) {
- myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
- nodes[1], nodes[3], nodes[5], nodes[7]));
- }
- else if ( nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 )) {
- const SMDS_MeshNode** tmpnodes = new const SMDS_MeshNode*[8];
- tmpnodes[0] = nodes[0];
- tmpnodes[1] = nodes[2];
- tmpnodes[2] = nodes[4];
- tmpnodes[3] = nodes[6];
- tmpnodes[4] = nodes[1];
- tmpnodes[5] = nodes[3];
- tmpnodes[6] = nodes[5];
- tmpnodes[7] = nodes[7];
- myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6],
- nodes[1], nodes[3], nodes[5], nodes[7]));
- aMesh->RemoveElement(f);
- }
+ break;
+ }
+ case 6: { /////// quadratic triangle
+ const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4],
+ nodes[1], nodes[3], nodes[5] );
+ if ( !f ||
+ nodes[2] != f->GetNodeWrap( f->GetNodeIndex( nodes[0] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[6] = { nodes[2 - 2*nextShift],
+ nodes[2],
+ nodes[2 + 2*nextShift],
+ nodes[3 - 2*nextShift],
+ nodes[3],
+ nodes[3 + 2*nextShift]};
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace( newOrder[ 0 ],
+ newOrder[ 1 ],
+ newOrder[ 2 ],
+ newOrder[ 3 ],
+ newOrder[ 4 ],
+ newOrder[ 5 ] ));
}
+ break;
}
- else { //////// polygon
- vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
- const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
- if ( !f )
- myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes));
- else if ( nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 1 ))
- {
- // TODO problem ChangeElementNodes : not the same number of nodes, not the same type
- MESSAGE("ChangeElementNodes");
- aMesh->ChangeElementNodes( f, nodes, nbn );
- }
+ default: /////// quadratic quadrangle
+ const SMDS_MeshFace * f = aMesh->FindFace( nodes[0], nodes[2], nodes[4], nodes[6],
+ nodes[1], nodes[3], nodes[5], nodes[7] );
+ if ( !f ||
+ nodes[ 2 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + 2*nextShift ))
+ {
+ const SMDS_MeshNode* newOrder[8] = { nodes[0],
+ nodes[4 - 2*nextShift],
+ nodes[4],
+ nodes[4 + 2*nextShift],
+ nodes[1],
+ nodes[5 - 2*nextShift],
+ nodes[5],
+ nodes[5 + 2*nextShift] };
+ if ( f )
+ aMesh->ChangeElementNodes( f, &newOrder[0], nbn );
+ else
+ myLastCreatedElems.Append(aMesh->AddFace(newOrder[ 0 ], newOrder[ 1 ],
+ newOrder[ 2 ], newOrder[ 3 ],
+ newOrder[ 4 ], newOrder[ 5 ],
+ newOrder[ 6 ], newOrder[ 7 ]));
+ }
+ } // switch ( nbn )
+
+ else { //////// polygon
+
+ vector<const SMDS_MeshNode*> polygon_nodes ( nodes, &nodes[nbn] );
+ const SMDS_MeshFace * f = aMesh->FindFace( polygon_nodes );
+ if ( !f ||
+ nodes[ 1 ] != f->GetNodeWrap( f->GetNodeIndex( nodes[ 0 ] ) + nextShift ))
+ {
+ if ( !vTool.IsForward() )
+ std::reverse( polygon_nodes.begin(), polygon_nodes.end());
+ if ( f )
+ aMesh->ChangeElementNodes( f, &polygon_nodes[0], nbn );
+ else
+ AddElement(polygon_nodes, SMDSAbs_Face, polygon_nodes.size()>4);
}
}
+
while ( srcElements.Length() < myLastCreatedElems.Length() )
srcElements.Append( *srcEdge );
// go to the next volume
iVol = 0;
while ( iVol++ < nbVolumesByStep ) v++;
- }
- }
+
+ } // loop on steps
+ } // loop on volumes of one step
} // sweep free links into faces
// Make a ceiling face with a normal external to a volume
}
//Extrusion_Error err =
MakeEdgePathPoints(aPrms, aTrackEdge, (aN1==theN1), fullList);
- }
- else if( aS.ShapeType() == TopAbs_WIRE ) {
+ } else if( aS.ShapeType() == TopAbs_WIRE ) {
list< SMESH_subMesh* > LSM;
TopTools_SequenceOfShape Edges;
SMESH_subMeshIteratorPtr itSM = theTrack->getDependsOnIterator(false,true);
list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
int startNid = theN1->GetID();
TColStd_MapOfInteger UsedNums;
+
int NbEdges = Edges.Length();
int i = 1;
for(; i<=NbEdges; i++) {
list<SMESH_MeshEditor_PathPoint> fullList;
const TopoDS_Shape& aS = theTrack->GetShapeToMesh();
- // Sub shape for the Pattern must be an Edge or Wire
- if( aS.ShapeType() == TopAbs_EDGE ) {
+
+ if( aS == SMESH_Mesh::PseudoShape() ) {
+ //Mesh without shape
+ const SMDS_MeshNode* currentNode = NULL;
+ const SMDS_MeshNode* prevNode = theN1;
+ std::vector<const SMDS_MeshNode*> aNodesList;
+ aNodesList.push_back(theN1);
+ int nbEdges = 0, conn=0;
+ const SMDS_MeshElement* prevElem = NULL;
+ const SMDS_MeshElement* currentElem = NULL;
+ int totalNbEdges = theTrack->NbEdges();
+ SMDS_ElemIteratorPtr nIt;
+ bool isClosed = false;
+
+ //check start node
+ if( !theTrack->GetMeshDS()->Contains(theN1) ) {
+ return EXTR_BAD_STARTING_NODE;
+ }
+
+ conn = nbEdgeConnectivity(theN1);
+ if(conn > 2)
+ return EXTR_PATH_NOT_EDGE;
+
+ aItE = theN1->GetInverseElementIterator();
+ prevElem = aItE->next();
+ currentElem = prevElem;
+ //Get all nodes
+ if(totalNbEdges == 1 ) {
+ nIt = currentElem->nodesIterator();
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+ } else {
+ nIt = currentElem->nodesIterator();
+ while( nIt->more() ) {
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ if(currentNode == prevNode)
+ currentNode = static_cast<const SMDS_MeshNode*>(nIt->next());
+ aNodesList.push_back(currentNode);
+
+ //case of the closed mesh
+ if(currentNode == theN1) {
+ nbEdges++;
+ isClosed = true;
+ break;
+ }
+
+ conn = nbEdgeConnectivity(currentNode);
+ if(conn > 2) {
+ return EXTR_PATH_NOT_EDGE;
+ }else if( conn == 1 && nbEdges > 0 ) {
+ //End of the path
+ nbEdges++;
+ break;
+ }else {
+ prevNode = currentNode;
+ aItE = currentNode->GetInverseElementIterator();
+ currentElem = aItE->next();
+ if( currentElem == prevElem)
+ currentElem = aItE->next();
+ nIt = currentElem->nodesIterator();
+ prevElem = currentElem;
+ nbEdges++;
+ }
+ }
+ }
+
+ if(nbEdges != totalNbEdges)
+ return EXTR_PATH_NOT_EDGE;
+
+ TopTools_SequenceOfShape Edges;
+ double x1,x2,y1,y2,z1,z2;
+ list< list<SMESH_MeshEditor_PathPoint> > LLPPs;
+ int startNid = theN1->GetID();
+ for(int i = 1; i < aNodesList.size(); i++) {
+ x1 = aNodesList[i-1]->X();x2 = aNodesList[i]->X();
+ y1 = aNodesList[i-1]->Y();y2 = aNodesList[i]->Y();
+ z1 = aNodesList[i-1]->Z();z2 = aNodesList[i]->Z();
+ TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2));
+ list<SMESH_MeshEditor_PathPoint> LPP;
+ aPrms.clear();
+ MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP);
+ LLPPs.push_back(LPP);
+ if( aNodesList[i-1]->GetID() == startNid ) startNid = aNodesList[i]->GetID();
+ else startNid = aNodesList[i-1]->GetID();
+
+ }
+
+ list< list<SMESH_MeshEditor_PathPoint> >::iterator itLLPP = LLPPs.begin();
+ list<SMESH_MeshEditor_PathPoint> firstList = *itLLPP;
+ list<SMESH_MeshEditor_PathPoint>::iterator itPP = firstList.begin();
+ for(; itPP!=firstList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+
+ SMESH_MeshEditor_PathPoint PP1 = fullList.back();
+ SMESH_MeshEditor_PathPoint PP2;
+ fullList.pop_back();
+ itLLPP++;
+ for(; itLLPP!=LLPPs.end(); itLLPP++) {
+ list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
+ itPP = currList.begin();
+ PP2 = currList.front();
+ gp_Dir D1 = PP1.Tangent();
+ gp_Dir D2 = PP2.Tangent();
+ gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
+ (D1.Z()+D2.Z())/2 ) );
+ PP1.SetTangent(Dnew);
+ fullList.push_back(PP1);
+ itPP++;
+ for(; itPP!=currList.end(); itPP++) {
+ fullList.push_back( *itPP );
+ }
+ PP1 = fullList.back();
+ fullList.pop_back();
+ }
+ fullList.push_back(PP1);
+
+ } // Sub shape for the Pattern must be an Edge or Wire
+ else if( aS.ShapeType() == TopAbs_EDGE ) {
aTrackEdge = TopoDS::Edge( aS );
// the Edge must not be degenerated
if ( BRep_Tool::Degenerated( aTrackEdge ) )
list<SMESH_MeshEditor_PathPoint> currList = *itLLPP;
itPP = currList.begin();
SMESH_MeshEditor_PathPoint PP2 = currList.front();
- gp_Pnt P1 = PP1.Pnt();
- //cout<<" PP1: Pnt("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
- gp_Pnt P2 = PP2.Pnt();
gp_Dir D1 = PP1.Tangent();
gp_Dir D2 = PP2.Tangent();
gp_Dir Dnew( gp_Vec( (D1.X()+D2.X())/2, (D1.Y()+D2.Y())/2,
* \param theCopy - if true, create translated copies of theElems
* \param theMakeGroups - if true and theCopy, create translated groups
* \param theTargetMesh - mesh to copy translated elements into
- * \retval SMESH_MeshEditor::PGroupIDs - list of ids of created groups
+ * \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups
*/
//================================================================================
PGroupIDs newGroupIDs;
- if ( theMakeGroups && theCopy ||
- theMakeGroups && theTargetMesh )
+ if ( ( theMakeGroups && theCopy ) ||
+ ( theMakeGroups && theTargetMesh ) )
newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh );
return newGroupIDs;
//=======================================================================
/*!
- * \brief Return SMESH_ElementSearcher
+ * \brief Return SMESH_ElementSearcher acting on a sub-set of elements
*/
//=======================================================================
}
curNodes[ iCur ] = n;
bool isUnique = nodeSet.insert( n ).second;
- if ( isUnique )
+ if ( isUnique ) {
uniqueNodes[ iUnique++ ] = n;
+ if ( nbRepl && iRepl[ nbRepl-1 ] == iCur )
+ --nbRepl; // n do not stick to a node of the elem
+ }
iCur++;
}
}
continue;
- }
+ } // poly element
// Regular elements
// TODO not all the possible cases are solved. Find something more generic?
isOk = false;
SMDS_VolumeTool hexa (elem);
hexa.SetExternalNormal();
- if ( nbUniqueNodes == 4 && nbRepl == 6 ) {
- //////////////////////// ---> tetrahedron
+ if ( nbUniqueNodes == 4 && nbRepl == 4 ) {
+ //////////////////////// HEX ---> 1 tetrahedron
for ( int iFace = 0; iFace < 6; iFace++ ) {
const int *ind = hexa.GetFaceNodesIndices( iFace ); // indices of face nodes
if (curNodes[ind[ 0 ]] == curNodes[ind[ 1 ]] &&
int iOppFace = hexa.GetOppFaceIndex( iFace );
ind = hexa.GetFaceNodesIndices( iOppFace );
int nbStick = 0;
- iUnique = 2; // reverse a tetrahedron bottom
for ( iCur = 0; iCur < 4 && nbStick < 2; iCur++ ) {
if ( curNodes[ind[ iCur ]] == curNodes[ind[ iCur + 1 ]] )
nbStick++;
- else if ( iUnique >= 0 )
- uniqueNodes[ iUnique-- ] = curNodes[ind[ iCur ]];
}
if ( nbStick == 1 ) {
// ... and the opposite one - into a triangle.
}
}
}
+ else if ( nbUniqueNodes == 6 && nbRepl == 2 ) {
+ //////////////////////// HEX ---> 1 prism
+ int nbTria = 0, iTria[3];
+ const int *ind; // indices of face nodes
+ // look for triangular faces
+ for ( int iFace = 0; iFace < 6 && nbTria < 3; iFace++ ) {
+ ind = hexa.GetFaceNodesIndices( iFace );
+ TIDSortedNodeSet faceNodes;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ faceNodes.insert( curNodes[ind[iCur]] );
+ if ( faceNodes.size() == 3 )
+ iTria[ nbTria++ ] = iFace;
+ }
+ // check if triangles are opposite
+ if ( nbTria == 2 && iTria[0] == hexa.GetOppFaceIndex( iTria[1] ))
+ {
+ isOk = true;
+ // set nodes of the bottom triangle
+ ind = hexa.GetFaceNodesIndices( iTria[ 0 ]);
+ vector<int> indB;
+ for ( iCur = 0; iCur < 4; iCur++ )
+ if ( ind[iCur] != iRepl[0] && ind[iCur] != iRepl[1])
+ indB.push_back( ind[iCur] );
+ if ( !hexa.IsForward() )
+ std::swap( indB[0], indB[2] );
+ for ( iCur = 0; iCur < 3; iCur++ )
+ uniqueNodes[ iCur ] = curNodes[indB[iCur]];
+ // set nodes of the top triangle
+ const int *indT = hexa.GetFaceNodesIndices( iTria[ 1 ]);
+ for ( iCur = 0; iCur < 3; ++iCur )
+ for ( int j = 0; j < 4; ++j )
+ if ( hexa.IsLinked( indB[ iCur ], indT[ j ] ))
+ {
+ uniqueNodes[ iCur + 3 ] = curNodes[ indT[ j ]];
+ break;
+ }
+ }
+ break;
+ }
else if (nbUniqueNodes == 5 && nbRepl == 4 ) {
//////////////////// HEXAHEDRON ---> 2 tetrahedrons
for ( int iFace = 0; iFace < 6; iFace++ ) {
}
}
} // if ( nbUniqueNodes == 6 && nbRepl == 4 )
+ else
+ {
+ MESSAGE("MergeNodes() removes hexahedron "<< elem);
+ }
break;
} // HEXAHEDRON
} // if ( nbNodes != nbUniqueNodes ) // some nodes stick
- if ( isOk ) {
- if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume) {
+ if ( isOk ) { // the elem remains valid after sticking nodes
+ if (elem->IsPoly() && elem->GetType() == SMDSAbs_Volume)
+ {
// Change nodes of polyedre
const SMDS_VtkVolume* aPolyedre =
dynamic_cast<const SMDS_VtkVolume*>( elem );
aMesh->ChangePolyhedronNodes( elem, poly_nodes, quantities );
}
}
- else {
- //int elemId = elem->GetID();
- //MESSAGE("Change regular element or polygon " << elemId);
- SMDSAbs_ElementType etyp = elem->GetType();
+ else // replace non-polyhedron elements
+ {
+ const SMDSAbs_ElementType etyp = elem->GetType();
+ const int elemId = elem->GetID();
+ const bool isPoly = (elem->GetEntityType() == SMDSEntity_Polygon);
uniqueNodes.resize(nbUniqueNodes);
- SMDS_MeshElement* newElem = 0;
- if (elem->GetEntityType() == SMDSEntity_Polygon)
- newElem = this->AddElement(uniqueNodes, etyp, true);
- else
- newElem = this->AddElement(uniqueNodes, etyp, false);
- if (newElem)
- {
- myLastCreatedElems.Append(newElem);
- if ( aShapeId )
- aMesh->SetMeshElementOnShape( newElem, aShapeId );
- }
- aMesh->RemoveElement(elem);
+
+ SMESHDS_SubMesh * sm = aShapeId > 0 ? aMesh->MeshElements(aShapeId) : 0;
+
+ aMesh->RemoveFreeElement(elem, sm, /*fromGroups=*/false);
+ SMDS_MeshElement* newElem = this->AddElement(uniqueNodes, etyp, isPoly, elemId);
+ if ( sm && newElem )
+ sm->AddElement( newElem );
+ if ( elem != newElem )
+ ReplaceElemInGroups( elem, newElem, aMesh );
}
}
else {
// Remove invalid regular element or invalid polygon
- //MESSAGE("Remove invalid " << elem->GetID());
rmElemIds.push_back( elem->GetID() );
}
//=======================================================================
/*!
* \brief Convert elements contained in a submesh to quadratic
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
*/
//=======================================================================
if( !elem || elem->IsQuadratic() ) continue;
int id = elem->GetID();
- //MESSAGE("elem " << id);
- id = 0; // get a free number for new elements
int nbNodes = elem->NbNodes();
SMDSAbs_ElementType aType = elem->GetType();
if ( elem->GetEntityType() == SMDSEntity_Polyhedra )
nbNodeInFaces = static_cast<const SMDS_VtkVolume* >( elem )->GetQuantities();
+ GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
+
const SMDS_MeshElement* NewElem = 0;
switch( aType )
ReplaceElemInGroups( elem, NewElem, GetMeshDS());
if( NewElem )
theSm->AddElement( NewElem );
-
- GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false);
}
// if (!GetMeshDS()->isCompacted())
// GetMeshDS()->compactMesh();
}
}
- if ( !theForce3d && !getenv("NO_FixQuadraticElements"))
+ if ( !theForce3d )
{ // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
aHelper.FixQuadraticElements();
}
- if (!GetMeshDS()->isCompacted())
- GetMeshDS()->compactMesh();
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements quadratic
+ * \param theForce3d - if true, the medium nodes will be placed in the middle of link
+ * \param theElements - elements to make quadratic
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d,
+ TIDSortedElemSet& theElements)
+{
+ if ( theElements.empty() ) return;
+
+ // we believe that all theElements are of the same type
+ SMDSAbs_ElementType elemType = (*theElements.begin())->GetType();
+
+ // get all nodes shared by theElements
+ TIDSortedNodeSet allNodes;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() );
+
+ // complete theElements with elements of lower dim whose all nodes are in allNodes
+
+ TIDSortedElemSet quadAdjacentElems [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements
+ TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ];
+ TIDSortedNodeSet::iterator nIt = allNodes.begin();
+ for ( ; nIt != allNodes.end(); ++nIt )
+ {
+ const SMDS_MeshNode* n = *nIt;
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ if ( e->IsQuadratic() )
+ {
+ quadAdjacentElems[ e->GetType() ].insert( e );
+ continue;
+ }
+ if ( e->GetType() >= elemType )
+ {
+ continue; // same type of more complex linear element
+ }
+
+ if ( !checkedAdjacentElems[ e->GetType() ].insert( e ).second )
+ continue; // e is already checked
+
+ // check nodes
+ bool allIn = true;
+ SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
+ while ( nodeIt->more() && allIn )
+ allIn = allNodes.count( cast2Node( nodeIt->next() ));
+ if ( allIn )
+ theElements.insert(e );
+ }
+ }
+
+ SMESH_MesherHelper helper(*myMesh);
+ helper.SetIsQuadratic( true );
+
+ // add links of quadratic adjacent elements to the helper
+
+ if ( !quadAdjacentElems[SMDSAbs_Edge].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Edge].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Face].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Face].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) );
+ }
+ if ( !quadAdjacentElems[SMDSAbs_Volume].empty() )
+ for ( eIt = quadAdjacentElems[SMDSAbs_Volume].begin();
+ eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt )
+ {
+ helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) );
+ }
+
+ // make quadratic elements instead of linear ones
+
+ SMESHDS_Mesh* meshDS = GetMeshDS();
+ SMESHDS_SubMesh* smDS = 0;
+ for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* elem = *eIt;
+ if( elem->IsQuadratic() || elem->NbNodes() < 2 || elem->IsPoly() )
+ continue;
+
+ int id = elem->GetID();
+ SMDSAbs_ElementType type = elem->GetType();
+ vector<const SMDS_MeshNode *> nodes ( elem->begin_nodes(), elem->end_nodes());
+
+ if ( !smDS || !smDS->Contains( elem ))
+ smDS = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false);
+
+ SMDS_MeshElement * newElem = 0;
+ switch( nodes.size() )
+ {
+ case 4: // cases for most multiple element types go first (for optimization)
+ if ( type == SMDSAbs_Volume )
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ else
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d);
+ break;
+ case 8:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d);
+ break;
+ case 3:
+ newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], id, theForce3d);
+ break;
+ case 2:
+ newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d);
+ break;
+ case 5:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], id, theForce3d);
+ break;
+ case 6:
+ newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3],
+ nodes[4], nodes[5], id, theForce3d);
+ break;
+ default:;
+ }
+ ReplaceElemInGroups( elem, newElem, meshDS);
+ if( newElem && smDS )
+ smDS->AddElement( newElem );
+ }
+
+ if ( !theForce3d && !getenv("NO_FixQuadraticElements"))
+ { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+ helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
+ helper.FixQuadraticElements();
+ }
}
//=======================================================================
/*!
* \brief Convert quadratic elements to linear ones and remove quadratic nodes
- * \retval int - nb of checked elements
+ * \return int - nb of checked elements
*/
//=======================================================================
{
int nbElem = 0;
SMESHDS_Mesh* meshDS = GetMeshDS();
- const bool notFromGroups = false;
while( theItr->more() )
{
nbElem++;
if( elem && elem->IsQuadratic())
{
- int id = elem->GetID();
- int nbNodes = elem->NbNodes();
- vector<const SMDS_MeshNode *> nodes, mediumNodes;
- nodes.reserve( nbNodes );
- mediumNodes.reserve( nbNodes );
-
- for(int i = 0; i < nbNodes; i++)
- {
- const SMDS_MeshNode* n = elem->GetNode(i);
-
- if( elem->IsMediumNode( n ) )
- mediumNodes.push_back( n );
- else
- nodes.push_back( n );
- }
- if( nodes.empty() ) continue;
+ int id = elem->GetID();
+ int nbCornerNodes = elem->NbCornerNodes();
SMDSAbs_ElementType aType = elem->GetType();
- //remove old quadratic element
- meshDS->RemoveFreeElement( elem, theSm, notFromGroups );
+ vector<const SMDS_MeshNode *> nodes( elem->begin_nodes(), elem->end_nodes() );
- SMDS_MeshElement * NewElem = AddElement( nodes, aType, false, id );
- ReplaceElemInGroups(elem, NewElem, meshDS);
- if( theSm && NewElem )
- theSm->AddElement( NewElem );
+ //remove a quadratic element
+ if ( !theSm || !theSm->Contains( elem ))
+ theSm = meshDS->MeshElements( elem->getshapeId() );
+ meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false );
// remove medium nodes
- vector<const SMDS_MeshNode*>::iterator nIt = mediumNodes.begin();
- for ( ; nIt != mediumNodes.end(); ++nIt ) {
- const SMDS_MeshNode* n = *nIt;
- if ( n->NbInverseElements() == 0 ) {
- if ( n->getshapeId() != theShapeID )
- meshDS->RemoveFreeNode( n, meshDS->MeshElements
- ( n->getshapeId() ));
- else
- meshDS->RemoveFreeNode( n, theSm );
- }
- }
+ for ( unsigned i = nbCornerNodes; i < nodes.size(); ++i )
+ if ( nodes[i]->NbInverseElements() == 0 )
+ meshDS->RemoveFreeNode( nodes[i], theSm );
+
+ // add a linear element
+ nodes.resize( nbCornerNodes );
+ SMDS_MeshElement * newElem = AddElement( nodes, aType, false, id );
+ ReplaceElemInGroups(elem, newElem, meshDS);
+ if( theSm && newElem )
+ theSm->AddElement( newElem );
}
}
return nbElem;
//function : ConvertFromQuadratic
//purpose :
//=======================================================================
-bool SMESH_MeshEditor::ConvertFromQuadratic()
+
+bool SMESH_MeshEditor::ConvertFromQuadratic()
{
int nbCheckedElems = 0;
if ( myMesh->HasShapeToMesh() )
return true;
}
+namespace
+{
+ //================================================================================
+ /*!
+ * \brief Return true if all medium nodes of the element are in the node set
+ */
+ //================================================================================
+
+ bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet )
+ {
+ for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i )
+ if ( !nodeSet.count( elem->GetNode(i) ))
+ return false;
+ return true;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Makes given elements linear
+ */
+//================================================================================
+
+void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements)
+{
+ if ( theElements.empty() ) return;
+
+ // collect IDs of medium nodes of theElements; some of these nodes will be removed
+ set<int> mediumNodeIDs;
+ TIDSortedElemSet::iterator eIt = theElements.begin();
+ for ( ; eIt != theElements.end(); ++eIt )
+ {
+ const SMDS_MeshElement* e = *eIt;
+ for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i )
+ mediumNodeIDs.insert( e->GetNode(i)->GetID() );
+ }
+
+ // replace given elements by linear ones
+ typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::iterator> TSetIterator;
+ SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() ));
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+
+ // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs
+ // except those elements sharing medium nodes of quadratic element whose medium nodes
+ // are not all in mediumNodeIDs
+
+ // get remaining medium nodes
+ TIDSortedNodeSet mediumNodes;
+ set<int>::iterator nIdsIt = mediumNodeIDs.begin();
+ for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt )
+ if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt ))
+ mediumNodes.insert( mediumNodes.end(), n );
+
+ // find more quadratic elements to convert
+ TIDSortedElemSet moreElemsToConvert;
+ TIDSortedNodeSet::iterator nIt = mediumNodes.begin();
+ for ( ; nIt != mediumNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e = invIt->next();
+ if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes ))
+ {
+ // find a more complex element including e and
+ // whose medium nodes are not in mediumNodes
+ bool complexFound = false;
+ for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type )
+ {
+ SMDS_ElemIteratorPtr invIt2 =
+ (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type ));
+ while ( invIt2->more() )
+ {
+ const SMDS_MeshElement* eComplex = invIt2->next();
+ if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes))
+ {
+ int nbCommonNodes = SMESH_Algo::GetCommonNodes( e, eComplex ).size();
+ if ( nbCommonNodes == e->NbNodes())
+ {
+ complexFound = true;
+ type = SMDSAbs_NbElementTypes; // to quit from the outer loop
+ break;
+ }
+ }
+ }
+ }
+ if ( !complexFound )
+ moreElemsToConvert.insert( e );
+ }
+ }
+ }
+ elemIt = SMDS_ElemIteratorPtr
+ (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() ));
+ removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 );
+}
+
//=======================================================================
//function : SewSideElements
//purpose :
}
// check similarity of elements of the sides
- if (aResult == SEW_OK && ( face[0] && !face[1] ) || ( !face[0] && face[1] )) {
+ if (aResult == SEW_OK && (( face[0] && !face[1] ) || ( !face[0] && face[1] ))) {
MESSAGE("Correspondent face not found on side " << ( face[0] ? 1 : 0 ));
if ( nReplaceMap.size() == 2 ) { // faces on input nodes not found
aResult = ( face[0] ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES );
* \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1
* \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1
* \param nReplaceMap - output map of corresponding nodes
- * \retval bool - is a success or not
+ * \return bool - is a success or not
*/
//================================================================================
_extremum.Perform(aPnt);
if ( _extremum.IsDone() )
for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
+#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
+ _state = ( _extremum.SquareDistance(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+#else
_state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+#endif
}
TopAbs_State State() const
{
return DoubleNodes( theElems, theNodesNot, anAffected );
}
+/*!
+ * \brief compute an oriented angle between two planes defined by four points.
+ * The vector (p0,p1) defines the intersection of the 2 planes (p0,p1,g1) and (p0,p1,g2)
+ * @param p0 base of the rotation axe
+ * @param p1 extremity of the rotation axe
+ * @param g1 belongs to the first plane
+ * @param g2 belongs to the second plane
+ */
+double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2)
+{
+// MESSAGE(" p0: " << p0.X() << " " << p0.Y() << " " << p0.Z());
+// MESSAGE(" p1: " << p1.X() << " " << p1.Y() << " " << p1.Z());
+// MESSAGE(" g1: " << g1.X() << " " << g1.Y() << " " << g1.Z());
+// MESSAGE(" g2: " << g2.X() << " " << g2.Y() << " " << g2.Z());
+ gp_Vec vref(p0, p1);
+ gp_Vec v1(p0, g1);
+ gp_Vec v2(p0, g2);
+ gp_Vec n1 = vref.Crossed(v1);
+ gp_Vec n2 = vref.Crossed(v2);
+ return n2.AngleWithRef(n1, vref);
+}
+
/*!
* \brief Double nodes on shared faces between groups of volumes and create flat elements on demand.
* The list of groups must describe a partition of the mesh volumes.
* The nodes of the internal faces at the boundaries of the groups are doubled.
* In option, the internal faces are replaced by flat elements.
* Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * The flat elements are stored in groups of volumes.
* @param theElems - list of groups of volumes, where a group of volume is a set of
* SMDS_MeshElements sorted by Id.
* @param createJointElems - if TRUE, create the elements
MESSAGE("----------------------------------------------");
SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
- meshDS->BuildDownWardConnectivity(false);
+ meshDS->BuildDownWardConnectivity(true);
CHRONO(50);
SMDS_UnstructuredGrid *grid = meshDS->getGrid();
cellDomains.clear();
nodeDomains.clear();
std::map<int,int> emptyMap;
+ std::set<int> emptySet;
emptyMap.clear();
for (int idom = 0; idom < theElems.size(); idom++)
}
}
- MESSAGE("Number of shared faces " << faceDomains.size());
+ //MESSAGE("Number of shared faces " << faceDomains.size());
std::map<DownIdType, std::map<int, int>, DownIdCompare>::iterator itface;
// --- explore the shared faces domain by domain,
continue;
int vtkType = grid->GetCellType(vtkId);
int downId = grid->CellIdToDownId(vtkId);
+ if (downId < 0)
+ {
+ MESSAGE("doubleNodesOnGroupBoundaries: internal algorithm problem");
+ continue; // not OK at this stage of the algorithm:
+ //no cells created after BuildDownWardConnectivity
+ }
DownIdType aCell(downId, vtkType);
if (celldom.count(vtkId))
continue;
// for each shared face, get the nodes
// for each node, for each domain of the face, create a clone of the node
+ // --- edges at the intersection of 3 or 4 domains, with the order of domains to build
+ // junction elements of type prism or hexa. the key is the pair of nodesId (lower first)
+ // the value is the ordered domain ids. (more than 4 domains not taken into account)
+
+ std::map<std::vector<int>, std::vector<int> > edgesMultiDomains; // nodes of edge --> ordered domains
+ std::map<int, std::vector<int> > mutipleNodes; // nodes muti domains with domain order
+
for (int idomain = 0; idomain < theElems.size(); idomain++)
{
itface = faceDomains.begin();
std::set<int> oldNodes;
oldNodes.clear();
grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ bool isMultipleDetected = false;
std::set<int>::iterator itn = oldNodes.begin();
for (; itn != oldNodes.end(); ++itn)
{
{
int idom = itdom->first;
//MESSAGE(" domain " << idom);
- if (!nodeDomains[oldId].count(idom))
+ if (!nodeDomains[oldId].count(idom)) // --- node to clone
{
+ if (nodeDomains[oldId].size() >= 2) // a multiple node
+ {
+ vector<int> orderedDoms;
+ //MESSAGE("multiple node " << oldId);
+ isMultipleDetected =true;
+ if (mutipleNodes.count(oldId))
+ orderedDoms = mutipleNodes[oldId];
+ else
+ {
+ map<int,int>::iterator it = nodeDomains[oldId].begin();
+ for (; it != nodeDomains[oldId].end(); ++it)
+ orderedDoms.push_back(it->first);
+ }
+ orderedDoms.push_back(idom); // TODO order ==> push_front or back
+ //stringstream txt;
+ //for (int i=0; i<orderedDoms.size(); i++)
+ // txt << orderedDoms[i] << " ";
+ //MESSAGE("orderedDoms " << txt.str());
+ mutipleNodes[oldId] = orderedDoms;
+ }
double *coords = grid->GetPoint(oldId);
SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]);
int newId = newNode->getVtkId();
nodeDomains[oldId][idom] = newId; // cloned node for other domains
- //MESSAGE(" newNode " << newId);
+ //MESSAGE(" newNode " << newId << " oldNode " << oldId << " size=" <<nodeDomains[oldId].size());
+ }
+ if (nodeDomains[oldId].size() >= 3)
+ {
+ //MESSAGE("confirm multiple node " << oldId);
+ isMultipleDetected =true;
+ }
+ }
+ }
+ if (isMultipleDetected) // check if an edge of the face is shared between 3 or more domains
+ {
+ //MESSAGE("multiple Nodes detected on a shared face");
+ int downId = itface->first.cellId;
+ unsigned char cellType = itface->first.cellType;
+ int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId);
+ const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId);
+ const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId);
+ for (int ie =0; ie < nbEdges; ie++)
+ {
+ int nodes[3];
+ int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes);
+ if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1]))
+ {
+ vector<int> vn0 = mutipleNodes[nodes[0]];
+ vector<int> vn1 = mutipleNodes[nodes[nbNodes - 1]];
+ sort( vn0.begin(), vn0.end() );
+ sort( vn1.begin(), vn1.end() );
+ if (vn0 == vn1)
+ {
+ //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]);
+ double *coords = grid->GetPoint(nodes[0]);
+ gp_Pnt p0(coords[0], coords[1], coords[2]);
+ coords = grid->GetPoint(nodes[nbNodes - 1]);
+ gp_Pnt p1(coords[0], coords[1], coords[2]);
+ gp_Pnt gref;
+ int vtkVolIds[1000]; // an edge can belong to a lot of volumes
+ map<int, SMDS_VtkVolume*> domvol; // domain --> a volume with the edge
+ map<int, double> angleDom; // oriented angles between planes defined by edge and volume centers
+ int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]);
+ for (int id=0; id < vn0.size(); id++)
+ {
+ int idom = vn0[id];
+ for (int ivol=0; ivol<nbvol; ivol++)
+ {
+ int smdsId = meshDS->fromVtkToSmds(vtkVolIds[ivol]);
+ SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId);
+ if (theElems[idom].count(elem))
+ {
+ SMDS_VtkVolume* svol = dynamic_cast<SMDS_VtkVolume*>(elem);
+ domvol[idom] = svol;
+ //MESSAGE(" domain " << idom << " volume " << elem->GetID());
+ double values[3];
+ vtkIdType npts = 0;
+ vtkIdType* pts = 0;
+ grid->GetCellPoints(vtkVolIds[ivol], npts, pts);
+ SMDS_VtkVolume::gravityCenter(grid, pts, npts, values);
+ if (id ==0)
+ {
+ gref.SetXYZ(gp_XYZ(values[0], values[1], values[2]));
+ angleDom[idom] = 0;
+ }
+ else
+ {
+ gp_Pnt g(values[0], values[1], values[2]);
+ angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pi<angle<+pi
+ //MESSAGE(" angle=" << angleDom[idom]);
+ }
+ break;
+ }
+ }
+ }
+ map<double, int> sortedDom; // sort domains by angle
+ for (map<int, double>::iterator ia = angleDom.begin(); ia != angleDom.end(); ++ia)
+ sortedDom[ia->second] = ia->first;
+ vector<int> vnodes;
+ vector<int> vdom;
+ for (map<double, int>::iterator ib = sortedDom.begin(); ib != sortedDom.end(); ++ib)
+ {
+ vdom.push_back(ib->second);
+ //MESSAGE(" ordered domain " << ib->second << " angle " << ib->first);
+ }
+ for (int ino = 0; ino < nbNodes; ino++)
+ vnodes.push_back(nodes[ino]);
+ edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains
+ }
}
}
}
// get node id's of the face (id SMDS = id VTK)
// create flat element with old and new nodes if requested
+ // --- new quad nodes on flat quad elements: oldId --> ((domain1 X domain2) --> newId)
+ // (domain1 X domain2) = domain1 + MAXINT*domain2
+
+ std::map<int, std::map<long,int> > nodeQuadDomains;
+ std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+
if (createJointElems)
{
itface = faceDomains.begin();
- for( ; itface != faceDomains.end();++itface )
+ for (; itface != faceDomains.end(); ++itface)
{
DownIdType face = itface->first;
std::set<int> oldNodes;
oldNodes.clear();
grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
- std::map<int,int> domvol = itface->second;
- std::map<int,int>::iterator itdom = domvol.begin();
+ std::map<int, int> domvol = itface->second;
+ std::map<int, int>::iterator itdom = domvol.begin();
int dom1 = itdom->first;
int vtkVolId = itdom->second;
itdom++;
int dom2 = itdom->first;
- meshDS->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains);
+ SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains,
+ nodeQuadDomains);
+ stringstream grpname;
+ grpname << "j_";
+ if (dom1 < dom2)
+ grpname << dom1 << "_" << dom2;
+ else
+ grpname << dom2 << "_" << dom1;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ }
+ }
+
+ // --- create volumes on multiple domain intersection if requested
+ // iterate on edgesMultiDomains
+
+ if (createJointElems)
+ {
+ std::map<std::vector<int>, std::vector<int> >::iterator ite = edgesMultiDomains.begin();
+ for (; ite != edgesMultiDomains.end(); ++ite)
+ {
+ vector<int> nodes = ite->first;
+ vector<int> orderDom = ite->second;
+ vector<vtkIdType> orderedNodes;
+ if (nodes.size() == 2)
+ {
+ //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]);
+ for (int ino=0; ino < nodes.size(); ino++)
+ if (orderDom.size() == 3)
+ for (int idom = 0; idom <orderDom.size(); idom++)
+ orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+ else
+ for (int idom = orderDom.size()-1; idom >=0; idom--)
+ orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] );
+ SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes);
+
+ stringstream grpname;
+ grpname << "mj_";
+ grpname << 0 << "_" << 0;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ }
+ else
+ {
+ //MESSAGE("Quadratic multiple joints not implemented");
+ // TODO quadratic nodes
+ }
+ }
+ }
+
+ // --- list the explicit faces and edges of the mesh that need to be modified,
+ // i.e. faces and edges built with one or more duplicated nodes.
+ // associate these faces or edges to their corresponding domain.
+ // only the first domain found is kept when a face or edge is shared
+
+ std::map<DownIdType, std::map<int,int>, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell)
+ std::map<int,int> feDom; // vtk id of cell to modify --> id domain
+ faceOrEdgeDom.clear();
+ feDom.clear();
+
+ for (int idomain = 0; idomain < theElems.size(); idomain++)
+ {
+ std::map<int, std::map<int, int> >::const_iterator itnod = nodeDomains.begin();
+ for (; itnod != nodeDomains.end(); ++itnod)
+ {
+ int oldId = itnod->first;
+ //MESSAGE(" node " << oldId);
+ vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId);
+ for (int i = 0; i < l.ncells; i++)
+ {
+ int vtkId = l.cells[i];
+ int vtkType = grid->GetCellType(vtkId);
+ int downId = grid->CellIdToDownId(vtkId);
+ if (downId < 0)
+ continue; // new cells: not to be modified
+ DownIdType aCell(downId, vtkType);
+ int volParents[1000];
+ int nbvol = grid->GetParentVolumes(volParents, vtkId);
+ for (int j = 0; j < nbvol; j++)
+ if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain))
+ if (!feDom.count(vtkId))
+ {
+ feDom[vtkId] = idomain;
+ faceOrEdgeDom[aCell] = emptyMap;
+ faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only
+ //MESSAGE("affect cell " << this->GetMeshDS()->fromVtkToSmds(vtkId) << " domain " << idomain
+ // << " type " << vtkType << " downId " << downId);
+ }
+ }
}
}
// get node id's of the face
// replace old nodes by new nodes in volumes, and update inverse connectivity
- MESSAGE("cellDomains " << cellDomains.size());
- faceDomains.insert(cellDomains.begin(), cellDomains.end());
- itface = faceDomains.begin();
- for( ; itface != faceDomains.end();++itface )
+ std::map<DownIdType, std::map<int,int>, DownIdCompare>* maps[3] = {&faceDomains, &cellDomains, &faceOrEdgeDom};
+ for (int m=0; m<3; m++)
{
- DownIdType face = itface->first;
- std::set<int> oldNodes;
- std::set<int>::iterator itn;
- oldNodes.clear();
- grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
- std::map<int,int> localClonedNodeIds;
-
- std::map<int,int> domvol = itface->second;
- std::map<int,int>::iterator itdom = domvol.begin();
- for(; itdom != domvol.end(); ++itdom)
+ std::map<DownIdType, std::map<int,int>, DownIdCompare>* amap = maps[m];
+ itface = (*amap).begin();
+ for (; itface != (*amap).end(); ++itface)
{
- int idom = itdom->first;
- int vtkVolId = itdom->second;
- localClonedNodeIds.clear();
- for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+ DownIdType face = itface->first;
+ std::set<int> oldNodes;
+ std::set<int>::iterator itn;
+ oldNodes.clear();
+ grid->GetNodeIds(oldNodes, face.cellId, face.cellType);
+ //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType));
+ std::map<int, int> localClonedNodeIds;
+
+ std::map<int, int> domvol = itface->second;
+ std::map<int, int>::iterator itdom = domvol.begin();
+ for (; itdom != domvol.end(); ++itdom)
{
- int oldId = *itn;
- if (nodeDomains[oldId].count(idom))
- localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+ int idom = itdom->first;
+ int vtkVolId = itdom->second;
+ //MESSAGE("modify nodes of cell " << this->GetMeshDS()->fromVtkToSmds(vtkVolId) << " domain " << idom);
+ localClonedNodeIds.clear();
+ for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn)
+ {
+ int oldId = *itn;
+ if (nodeDomains[oldId].count(idom))
+ {
+ localClonedNodeIds[oldId] = nodeDomains[oldId][idom];
+ //MESSAGE(" node " << oldId << " --> " << localClonedNodeIds[oldId]);
+ }
+ }
+ meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
}
- meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds);
}
}
+
+ meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory
grid->BuildLinks();
- // TODO replace also old nodes by new nodes in faces and edges
CHRONOSTOP(50);
counters::stats();
return true;
}
+/*!
+ * \brief Double nodes on some external faces and create flat elements.
+ * Flat elements are mainly used by some types of mechanic calculations.
+ *
+ * Each group of the list must be constituted of faces.
+ * Triangles are transformed in prisms, and quadrangles in hexahedrons.
+ * @param theElems - list of groups of faces, where a group of faces is a set of
+ * SMDS_MeshElements sorted by Id.
+ * @return TRUE if operation has been completed successfully, FALSE otherwise
+ */
+bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector<TIDSortedElemSet>& theElems)
+{
+ MESSAGE("-------------------------------------------------");
+ MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups");
+ MESSAGE("-------------------------------------------------");
+
+ SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS();
+
+ // --- For each group of faces
+ // duplicate the nodes, create a flat element based on the face
+ // replace the nodes of the faces by their clones
+
+ std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> clonedNodes;
+ std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> intermediateNodes;
+ clonedNodes.clear();
+ intermediateNodes.clear();
+ std::map<std::string, SMESH_Group*> mapOfJunctionGroups;
+ mapOfJunctionGroups.clear();
+
+ for (int idom = 0; idom < theElems.size(); idom++)
+ {
+ const TIDSortedElemSet& domain = theElems[idom];
+ TIDSortedElemSet::const_iterator elemItr = domain.begin();
+ for (; elemItr != domain.end(); ++elemItr)
+ {
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr;
+ SMDS_MeshFace* aFace = dynamic_cast<SMDS_MeshFace*> (anElem);
+ if (!aFace)
+ continue;
+ // MESSAGE("aFace=" << aFace->GetID());
+ bool isQuad = aFace->IsQuadratic();
+ vector<const SMDS_MeshNode*> ln0, ln1, ln2, ln3, ln4;
+
+ // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face
+
+ SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator();
+ while (nodeIt->more())
+ {
+ const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*> (nodeIt->next());
+ bool isMedium = isQuad && (aFace->IsMediumNode(node));
+ if (isMedium)
+ ln2.push_back(node);
+ else
+ ln0.push_back(node);
+
+ const SMDS_MeshNode* clone = 0;
+ if (!clonedNodes.count(node))
+ {
+ clone = meshDS->AddNode(node->X(), node->Y(), node->Z());
+ clonedNodes[node] = clone;
+ }
+ else
+ clone = clonedNodes[node];
+
+ if (isMedium)
+ ln3.push_back(clone);
+ else
+ ln1.push_back(clone);
+
+ const SMDS_MeshNode* inter = 0;
+ if (isQuad && (!isMedium))
+ {
+ if (!intermediateNodes.count(node))
+ {
+ inter = meshDS->AddNode(node->X(), node->Y(), node->Z());
+ intermediateNodes[node] = inter;
+ }
+ else
+ inter = intermediateNodes[node];
+ ln4.push_back(inter);
+ }
+ }
+
+ // --- extrude the face
+
+ vector<const SMDS_MeshNode*> ln;
+ SMDS_MeshVolume* vol = 0;
+ vtkIdType aType = aFace->GetVtkType();
+ switch (aType)
+ {
+ case VTK_TRIANGLE:
+ vol = meshDS->AddVolume(ln0[2], ln0[1], ln0[0], ln1[2], ln1[1], ln1[0]);
+ // MESSAGE("vol prism " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ break;
+ case VTK_QUAD:
+ vol = meshDS->AddVolume(ln0[3], ln0[2], ln0[1], ln0[0], ln1[3], ln1[2], ln1[1], ln1[0]);
+ // MESSAGE("vol hexa " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln1[3]);
+ break;
+ case VTK_QUADRATIC_TRIANGLE:
+ vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln0[0], ln0[1], ln0[2], ln3[0], ln3[1], ln3[2],
+ ln2[0], ln2[1], ln2[2], ln4[0], ln4[1], ln4[2]);
+ // MESSAGE("vol quad prism " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln3[0]);
+ ln.push_back(ln3[1]);
+ ln.push_back(ln3[2]);
+ break;
+ case VTK_QUADRATIC_QUAD:
+// vol = meshDS->AddVolume(ln0[0], ln0[1], ln0[2], ln0[3], ln1[0], ln1[1], ln1[2], ln1[3],
+// ln2[0], ln2[1], ln2[2], ln2[3], ln3[0], ln3[1], ln3[2], ln3[3],
+// ln4[0], ln4[1], ln4[2], ln4[3]);
+ vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln1[3], ln0[0], ln0[1], ln0[2], ln0[3],
+ ln3[0], ln3[1], ln3[2], ln3[3], ln2[0], ln2[1], ln2[2], ln2[3],
+ ln4[0], ln4[1], ln4[2], ln4[3]);
+ // MESSAGE("vol quad hexa " << vol->GetID());
+ ln.push_back(ln1[0]);
+ ln.push_back(ln1[1]);
+ ln.push_back(ln1[2]);
+ ln.push_back(ln1[3]);
+ ln.push_back(ln3[0]);
+ ln.push_back(ln3[1]);
+ ln.push_back(ln3[2]);
+ ln.push_back(ln3[3]);
+ break;
+ case VTK_POLYGON:
+ break;
+ default:
+ break;
+ }
+
+ if (vol)
+ {
+ stringstream grpname;
+ grpname << "jf_";
+ grpname << idom;
+ int idg;
+ string namegrp = grpname.str();
+ if (!mapOfJunctionGroups.count(namegrp))
+ mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg);
+ SMESHDS_Group *sgrp = dynamic_cast<SMESHDS_Group*>(mapOfJunctionGroups[namegrp]->GetGroupDS());
+ if (sgrp)
+ sgrp->Add(vol->GetID());
+ }
+
+ // --- modify the face
+
+ aFace->ChangeNodes(&ln[0], ln.size());
+ }
+ }
+ return true;
+}
+
//================================================================================
/*!
* \brief Generates skin mesh (containing 2D cells) from 3D mesh
* \param group - a group to store created boundary elements in
* \param targetMesh - a mesh to store created boundary elements in
* \param toCopyElements - if true, the checked elements will be copied into the targetMesh
- * \param toCopyExistingBondary - if true, not only new but also pre-existing
+ * \param toCopyExistingBoundary - if true, not only new but also pre-existing
* boundary elements will be copied into the targetMesh
+ * \param toAddExistingBondary - if true, not only new but also pre-existing
+ * boundary elements will be added into the new group
+ * \param aroundElements - if true, elements will be created on boundary of given
+ * elements else, on boundary of the whole mesh.
+ * \return nb of added boundary elements
*/
//================================================================================
-void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
- Bnd_Dimension dimension,
- SMESH_Group* group/*=0*/,
- SMESH_Mesh* targetMesh/*=0*/,
- bool toCopyElements/*=false*/,
- bool toCopyExistingBondary/*=false*/)
+int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements,
+ Bnd_Dimension dimension,
+ SMESH_Group* group/*=0*/,
+ SMESH_Mesh* targetMesh/*=0*/,
+ bool toCopyElements/*=false*/,
+ bool toCopyExistingBoundary/*=false*/,
+ bool toAddExistingBondary/*= false*/,
+ bool aroundElements/*= false*/)
{
SMDSAbs_ElementType missType = (dimension == BND_2DFROM3D) ? SMDSAbs_Face : SMDSAbs_Edge;
SMDSAbs_ElementType elemType = (dimension == BND_1DFROM2D) ? SMDSAbs_Face : SMDSAbs_Volume;
throw SALOME_Exception(LOCALIZED("wrong element type"));
if ( !targetMesh )
- toCopyElements = toCopyExistingBondary = false;
+ toCopyElements = toCopyExistingBoundary = false;
SMESH_MeshEditor tgtEditor( targetMesh ? targetMesh : myMesh );
SMESHDS_Mesh* aMesh = GetMeshDS(), *tgtMeshDS = tgtEditor.GetMeshDS();
+ int nbAddedBnd = 0;
+
+ // editor adding present bnd elements and optionally holding elements to add to the group
+ SMESH_MeshEditor* presentEditor;
+ SMESH_MeshEditor tgtEditor2( tgtEditor.GetMesh() );
+ presentEditor = toAddExistingBondary ? &tgtEditor : &tgtEditor2;
SMDS_VolumeTool vTool;
- TIDSortedElemSet emptySet, avoidSet;
+ TIDSortedElemSet avoidSet;
+ const TIDSortedElemSet emptySet, *elemSet = aroundElements ? &elements : &emptySet;
int inode;
typedef vector<const SMDS_MeshNode*> TConnectivity;
const SMDS_MeshElement* elem = eIt->next();
const int iQuad = elem->IsQuadratic();
+ // ------------------------------------------------------------------------------------
// 1. For an elem, get present bnd elements and connectivities of missing bnd elements
+ // ------------------------------------------------------------------------------------
vector<const SMDS_MeshElement*> presentBndElems;
vector<TConnectivity> missingBndElems;
TConnectivity nodes;
if ( vTool.Set(elem) ) // elem is a volume ------------------------------------------
{
vTool.SetExternalNormal();
+ const SMDS_MeshElement* otherVol = 0;
for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
{
- if (!vTool.IsFreeFace(iface))
+ if ( !vTool.IsFreeFace(iface, &otherVol) &&
+ ( !aroundElements || elements.count( otherVol )))
continue;
- int nbFaceNodes = vTool.NbFaceNodes(iface);
+ const int nbFaceNodes = vTool.NbFaceNodes(iface);
const SMDS_MeshNode** nn = vTool.GetFaceNodes(iface);
if ( missType == SMDSAbs_Edge ) // boundary edges
{
presentBndElems.push_back( f );
else
missingBndElems.push_back( nodes );
+
+ if ( targetMesh != myMesh )
+ {
+ // add 1D elements on face boundary to be added to a new mesh
+ const SMDS_MeshElement* edge;
+ for ( inode = 0; inode < nbFaceNodes; inode += 1+iQuad)
+ {
+ if ( iQuad )
+ edge = aMesh->FindEdge( nn[inode], nn[inode+1], nn[inode+2]);
+ else
+ edge = aMesh->FindEdge( nn[inode], nn[inode+1]);
+ if ( edge && avoidSet.insert( edge ).second )
+ presentBndElems.push_back( edge );
+ }
+ }
}
}
}
{
nodes[0] = elem->GetNode(i);
nodes[1] = elem->GetNode((i+1)%nbNodes);
- if ( FindFaceInSet( nodes[0], nodes[1], emptySet, avoidSet))
+ if ( FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet))
continue; // not free link
//if ( iQuad )
}
}
+ // ---------------------------------
// 2. Add missing boundary elements
+ // ---------------------------------
if ( targetMesh != myMesh )
// instead of making a map of nodes in this mesh and targetMesh,
- // we create nodes with same IDs. We can renumber them later, if needed
+ // we create nodes with same IDs.
for ( int i = 0; i < missingBndElems.size(); ++i )
{
TConnectivity& srcNodes = missingBndElems[i];
TConnectivity nodes( srcNodes.size() );
for ( inode = 0; inode < nodes.size(); ++inode )
nodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] );
+ if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+ missType,
+ /*noMedium=*/true))
+ continue;
tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+ ++nbAddedBnd;
}
else
for ( int i = 0; i < missingBndElems.size(); ++i )
{
- TConnectivity& nodes = missingBndElems[i];
+ TConnectivity& nodes = missingBndElems[i];
+ if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes,
+ missType,
+ /*noMedium=*/true))
+ continue;
tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4);
+ ++nbAddedBnd;
}
+ // ----------------------------------
// 3. Copy present boundary elements
- if ( toCopyExistingBondary )
+ // ----------------------------------
+ if ( toCopyExistingBoundary )
for ( int i = 0 ; i < presentBndElems.size(); ++i )
{
const SMDS_MeshElement* e = presentBndElems[i];
TConnectivity nodes( e->NbNodes() );
for ( inode = 0; inode < nodes.size(); ++inode )
nodes[inode] = getNodeWithSameID( tgtMeshDS, e->GetNode(inode) );
- tgtEditor.AddElement(nodes, missType, e->IsPoly());
- // leave only missing elements in tgtEditor.myLastCreatedElems
- tgtEditor.myLastCreatedElems.Remove( tgtEditor.myLastCreatedElems.Size() );
+ presentEditor->AddElement(nodes, e->GetType(), e->IsPoly());
+ }
+ else // store present elements to add them to a group
+ for ( int i = 0 ; i < presentBndElems.size(); ++i )
+ {
+ presentEditor->myLastCreatedElems.Append(presentBndElems[i]);
}
+
} // loop on given elements
- // 4. Fill group with missing boundary elements
+ // ---------------------------------------------
+ // 4. Fill group with boundary elements
+ // ---------------------------------------------
if ( group )
{
if ( SMESHDS_Group* g = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() ))
g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 ));
}
tgtEditor.myLastCreatedElems.Clear();
+ tgtEditor2.myLastCreatedElems.Clear();
+ // -----------------------
// 5. Copy given elements
- if ( toCopyElements )
+ // -----------------------
+ if ( toCopyElements && targetMesh != myMesh )
{
if (elements.empty())
eIt = aMesh->elementsIterator(elemType);
tgtEditor.myLastCreatedElems.Clear();
}
}
- return;
+ return nbAddedBnd;
}
