-// SMESH SMESH : implementaion of SMESH idl descriptions
+// Copyright (C) 2007-2008 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
+//
+// 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.
//
-// Copyright (C) 2003 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 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
+// 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
//
+// SMESH SMESH : implementaion of SMESH idl descriptions
// File : StdMeshers_Quadrangle_2D.cxx
// Moved here from SMESH_Quadrangle_2D.cxx
// Author : Paul RASCLE, EDF
// Module : SMESH
-// $Header$
-
+//
#include "StdMeshers_Quadrangle_2D.hxx"
#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_QuadrangleParams.hxx"
+
#include "SMESH_Gen.hxx"
#include "SMESH_Mesh.hxx"
#include "SMESH_subMesh.hxx"
#include "SMDS_EdgePosition.hxx"
#include "SMDS_FacePosition.hxx"
-#include <BRepTools.hxx>
#include <BRepTools_WireExplorer.hxx>
#include <BRep_Tool.hxx>
#include <Geom_Surface.hxx>
#include <TColStd_SequenceOfReal.hxx>
#include <TColgp_SequenceOfXY.hxx>
#include <TopExp.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopoDS.hxx>
#include "utilities.h"
*/
//=============================================================================
-StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId, SMESH_Gen* gen)
+StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
+ SMESH_Gen* gen)
: SMESH_2D_Algo(hypId, studyId, gen)
{
MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
_name = "Quadrangle_2D";
_shapeType = (1 << TopAbs_FACE);
+ _compatibleHypothesis.push_back("QuadrangleParams");
_compatibleHypothesis.push_back("QuadranglePreference");
+ _compatibleHypothesis.push_back("TrianglePreference");
myTool = 0;
}
bool isOk = true;
aStatus = SMESH_Hypothesis::HYP_OK;
- // there is only one compatible Hypothesis so far
- const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape, false);
- myQuadranglePreference = hyps.size() > 0;
+ const list <const SMESHDS_Hypothesis * >&hyps =
+ GetUsedHypothesis(aMesh, aShape, false);
+ const SMESHDS_Hypothesis *theHyp = 0;
+
+ if( hyps.size() == 1 ) {
+ myTriaVertexID = -1;
+ theHyp = hyps.front();
+ if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) {
+ const StdMeshers_QuadrangleParams* theHyp1 =
+ (const StdMeshers_QuadrangleParams*)theHyp;
+ myTriaVertexID = theHyp1->GetTriaVertex();
+ myQuadranglePreference= false;
+ myTrianglePreference= false;
+ }
+ if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
+ myQuadranglePreference= true;
+ myTrianglePreference= false;
+ myTriaVertexID = -1;
+ }
+ else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
+ myQuadranglePreference= false;
+ myTrianglePreference= true;
+ myTriaVertexID = -1;
+ }
+ }
+
+ else if( hyps.size() > 1 ) {
+ theHyp = hyps.front();
+ if(strcmp("QuadrangleParams", theHyp->GetName()) == 0) {
+ const StdMeshers_QuadrangleParams* theHyp1 =
+ (const StdMeshers_QuadrangleParams*)theHyp;
+ myTriaVertexID = theHyp1->GetTriaVertex();
+ theHyp = hyps.back();
+ if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
+ myQuadranglePreference= true;
+ myTrianglePreference= false;
+ }
+ else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
+ myQuadranglePreference= false;
+ myTrianglePreference= true;
+ }
+ }
+ else {
+ if(strcmp("QuadranglePreference", theHyp->GetName()) == 0) {
+ myQuadranglePreference= true;
+ myTrianglePreference= false;
+ }
+ else if(strcmp("TrianglePreference", theHyp->GetName()) == 0){
+ myQuadranglePreference= false;
+ myTrianglePreference= true;
+ }
+ const StdMeshers_QuadrangleParams* theHyp2 =
+ (const StdMeshers_QuadrangleParams*)hyps.back();
+ myTriaVertexID = theHyp2->GetTriaVertex();
+ }
+ }
+
+ else {
+ myQuadranglePreference = false;
+ myTrianglePreference = false;
+ myTriaVertexID = -1;
+ }
return isOk;
}
c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) {
+ meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
}
}
}
if (near == g) { // make triangle
- //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
SMDS_MeshFace* face = myTool->AddFace(a, b, c);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else { // make quadrangle
if (near - 1 < ilow)
else
d = quad->uv_grid[nbhoriz + near - 1].node;
//SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
- SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+
+ if(!myTrianglePreference){
+ SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else {
+ SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ }
// if node d is not at position g - make additional triangles
if (near - 1 > g) {
d = uv_e3[1].node;
else
d = quad->uv_grid[nbhoriz + k - 1].node;
- //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
SMDS_MeshFace* face = myTool->AddFace(a, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
}
g = near;
}
if (near == g) { // make triangle
- //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
SMDS_MeshFace* face = myTool->AddFace(a, b, c);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else { // make quadrangle
if (near + 1 > iup)
else
d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
//SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
- SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(!myTrianglePreference){
+ SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else {
+ SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ }
if (near + 1 < g) { // if d not is at g - make additional triangles
for (int k = near + 1; k < g; k++) {
d = uv_e1[nbright - 2].node;
else
d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
- //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
SMDS_MeshFace* face = myTool->AddFace(a, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
}
g = near;
}
if (near == g) { // make triangle
- //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
SMDS_MeshFace* face = myTool->AddFace(a, b, c);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else { // make quadrangle
if (near - 1 < jlow)
else
d = quad->uv_grid[nbhoriz*near - 2].node;
//SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
- SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+
+ if(!myTrianglePreference){
+ SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else {
+ SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ }
if (near - 1 > g) { // if d not is at g - make additional triangles
for (int k = near - 1; k > g; k--) {
d = uv_e0[nbdown - 2].node;
else
d = quad->uv_grid[nbhoriz*k - 2].node;
- //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
SMDS_MeshFace* face = myTool->AddFace(a, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
}
g = near;
}
if (near == g) { // make triangle
- //SMDS_MeshFace* face = meshDS->AddFace(a, b, c);
SMDS_MeshFace* face = myTool->AddFace(a, b, c);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
else { // make quadrangle
if (near + 1 > jup)
else
d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
//SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
- SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(!myTrianglePreference){
+ SMDS_MeshFace* face = myTool->AddFace(a, b, c, d);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
+ }
+ else {
+ SplitQuad(meshDS, geomFaceID, a, b, c, d);
+ }
if (near + 1 < g) { // if d not is at g - make additional triangles
for (int k = near + 1; k < g; k++) {
d = uv_e2[1].node;
else
d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
- //SMDS_MeshFace* face = meshDS->AddFace(a, c, d);
SMDS_MeshFace* face = myTool->AddFace(a, c, d);
- meshDS->SetMeshElementOnShape(face, geomFaceID);
+ if(face) meshDS->SetMeshElementOnShape(face, geomFaceID);
}
}
g = near;
return isOk;
}
+
+//=============================================================================
+/*!
+ * Evaluate
+ */
+//=============================================================================
+
+bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape,
+ MapShapeNbElems& aResMap)
+
+{
+ aMesh.GetSubMesh(aShape);
+
+ std::vector<int> aNbNodes(4);
+ bool IsQuadratic = false;
+ if( !CheckNbEdgesForEvaluate( aMesh, aShape, aResMap, aNbNodes, IsQuadratic ) ) {
+ std::vector<int> aResVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aResVec));
+ SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+ smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+ return false;
+ }
+
+ if(myQuadranglePreference) {
+ int n1 = aNbNodes[0];
+ int n2 = aNbNodes[1];
+ int n3 = aNbNodes[2];
+ int n4 = aNbNodes[3];
+ int nfull = n1+n2+n3+n4;
+ int ntmp = nfull/2;
+ ntmp = ntmp*2;
+ if( nfull==ntmp && ( (n1!=n3) || (n2!=n4) ) ) {
+ // special path for using only quandrangle faces
+ return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
+ //return true;
+ }
+ }
+
+ int nbdown = aNbNodes[0];
+ int nbup = aNbNodes[2];
+
+ int nbright = aNbNodes[1];
+ int nbleft = aNbNodes[3];
+
+ int nbhoriz = Min(nbdown, nbup);
+ int nbvertic = Min(nbright, nbleft);
+
+ int dh = Max(nbdown, nbup) - nbhoriz;
+ int dv = Max(nbright, nbleft) - nbvertic;
+
+ //int kdh = 0;
+ //if(dh>0) kdh = 1;
+ //int kdv = 0;
+ //if(dv>0) kdv = 1;
+
+ int nbNodes = (nbhoriz-2)*(nbvertic-2);
+ //int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2;
+ int nbFaces3 = dh + dv;
+ //if( kdh==1 && kdv==1 ) nbFaces3 -= 2;
+ //if( dh>0 && dv>0 ) nbFaces3 -= 2;
+ //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
+ int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
+
+ std::vector<int> aVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
+ if(IsQuadratic) {
+ aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
+ aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
+ int nbbndedges = nbdown + nbup + nbright + nbleft -4;
+ int nbintedges = ( nbFaces4*4 + nbFaces3*3 - nbbndedges ) / 2;
+ aVec[SMDSEntity_Node] = nbNodes + nbintedges;
+ if( aNbNodes.size()==5 ) {
+ aVec[SMDSEntity_Quad_Triangle] = nbFaces3 + aNbNodes[3] -1;
+ aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4 - aNbNodes[3] +1;
+ }
+ }
+ else {
+ aVec[SMDSEntity_Node] = nbNodes;
+ aVec[SMDSEntity_Triangle] = nbFaces3;
+ aVec[SMDSEntity_Quadrangle] = nbFaces4;
+ if( aNbNodes.size()==5 ) {
+ aVec[SMDSEntity_Triangle] = nbFaces3 + aNbNodes[3] - 1;
+ aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
+ }
+ }
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aVec));
+
+ return true;
+}
+
+
+//================================================================================
+/*!
+ * \brief Return true if only two given edges meat at their common vertex
+ */
+//================================================================================
+
+static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
+ const TopoDS_Edge& e2,
+ SMESH_Mesh & mesh)
+{
+ TopoDS_Vertex v;
+ if ( !TopExp::CommonVertex( e1, e2, v ))
+ return false;
+ TopTools_ListIteratorOfListOfShape ancestIt( mesh.GetAncestors( v ));
+ for ( ; ancestIt.More() ; ancestIt.Next() )
+ if ( ancestIt.Value().ShapeType() == TopAbs_EDGE )
+ if ( !e1.IsSame( ancestIt.Value() ) && !e2.IsSame( ancestIt.Value() ))
+ return false;
+ return true;
+}
+
//=============================================================================
/*!
*
int nbSides = 0;
list< TopoDS_Edge >::iterator edgeIt = edges.begin();
- if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
+ if ( nbEdgesInWire.front() == 3 ) { // exactly 3 edges
+ if(myTriaVertexID>0) {
+ SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
+ TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
+ if(!V.IsNull()) {
+ TopoDS_Edge E1,E2,E3;
+ for(; edgeIt != edges.end(); ++edgeIt) {
+ TopoDS_Edge E = TopoDS::Edge(*edgeIt);
+ TopoDS_Vertex VF, VL;
+ TopExp::Vertices(E, VF, VL, true);
+ if( VF.IsSame(V) )
+ E1 = E;
+ else if( VL.IsSame(V) )
+ E3 = E;
+ else
+ E2 = E;
+ }
+ quad->side.reserve(4);
+ quad->side.push_back( new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes));
+ quad->side.push_back( new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes));
+ quad->side.push_back( new StdMeshers_FaceSide(F, E3, &aMesh, false, ignoreMediumNodes));
+ std::vector<UVPtStruct> UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
+ std::vector<UVPtStruct> UVPStop = quad->side[1]->GetUVPtStruct(false,1);
+ std::vector<UVPtStruct> UVPSright = quad->side[2]->GetUVPtStruct(true,1);
+ const SMDS_MeshNode* aNode = UVPSleft[0].node;
+ gp_Pnt2d aPnt2d( UVPSleft[0].u, UVPSleft[0].v );
+ StdMeshers_FaceSide* VertFS =
+ new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]);
+ quad->side.push_back(VertFS);
+ return quad;
+ }
+ }
+ return 0;
+ }
+ else if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
for ( ; edgeIt != edges.end(); ++edgeIt, nbSides++ )
quad->side.push_back( new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
nbSides<TOP_SIDE, ignoreMediumNodes));
nbSides<TOP_SIDE, ignoreMediumNodes));
++nbSides;
}
+ // issue 20222. Try to unite only edges shared by two same faces
+ if (nbSides < 4) {
+ // delete found sides
+ { FaceQuadStruct cleaner( *quad ); }
+ quad->side.clear();
+ quad->side.reserve(nbEdgesInWire.front());
+ nbSides = 0;
+
+ SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ while ( !edges.empty()) {
+ sideEdges.clear();
+ sideEdges.splice( sideEdges.end(), edges, edges.begin());
+ bool sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide =
+ SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
+ twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.end(), edges, edges.begin());
+ }
+ if ( nbSides == 0 ) { // go backward from the first edge
+ sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide =
+ SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
+ twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.begin(), edges, --edges.end());
+ }
+ }
+ quad->side.push_back( new StdMeshers_FaceSide(F, sideEdges, &aMesh,
+ nbSides<TOP_SIDE, ignoreMediumNodes));
+ ++nbSides;
+ }
+ }
}
if (nbSides != 4) {
#ifdef _DEBUG_
- cout << endl << "StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:";
+ MESSAGE ( "StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n" );
for ( int i = 0; i < nbSides; ++i ) {
- cout << " ( ";
+ MESSAGE ( " ( " );
for ( int e = 0; e < quad->side[i]->NbEdges(); ++e )
- cout << myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " ";
- cout << ")";
+ MESSAGE ( myTool->GetMeshDS()->ShapeToIndex( quad->side[i]->Edge( e )) << " " );
+ MESSAGE ( ")\n" );
}
- cout << endl;
+ //cout << endl;
#endif
if ( !nbSides )
nbSides = nbEdgesInWire.front();
return quad;
}
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+
+bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
+ const TopoDS_Shape & aShape,
+ MapShapeNbElems& aResMap,
+ std::vector<int>& aNbNodes,
+ bool& IsQuadratic)
+
+{
+ const TopoDS_Face & F = TopoDS::Face(aShape);
+
+ // verify 1 wire only, with 4 edges
+ TopoDS_Vertex V;
+ list< TopoDS_Edge > edges;
+ list< int > nbEdgesInWire;
+ int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ if (nbWire != 1) {
+ return false;
+ }
+
+ aNbNodes.resize(4);
+
+ int nbSides = 0;
+ list< TopoDS_Edge >::iterator edgeIt = edges.begin();
+ SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) {
+ return false;
+ }
+ std::vector<int> aVec = (*anIt).second;
+ IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+ if ( nbEdgesInWire.front() == 3 ) { // exactly 3 edges
+ if(myTriaVertexID>0) {
+ SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
+ TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
+ if(!V.IsNull()) {
+ TopoDS_Edge E1,E2,E3;
+ for(; edgeIt != edges.end(); ++edgeIt) {
+ TopoDS_Edge E = TopoDS::Edge(*edgeIt);
+ TopoDS_Vertex VF, VL;
+ TopExp::Vertices(E, VF, VL, true);
+ if( VF.IsSame(V) )
+ E1 = E;
+ else if( VL.IsSame(V) )
+ E3 = E;
+ else
+ E2 = E;
+ }
+ SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) return false;
+ std::vector<int> aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
+ else
+ aNbNodes[0] = aVec[SMDSEntity_Node] + 2;
+ sm = aMesh.GetSubMesh(E2);
+ anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) return false;
+ aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[1] = (aVec[SMDSEntity_Node]-1)/2 + 2;
+ else
+ aNbNodes[1] = aVec[SMDSEntity_Node] + 2;
+ sm = aMesh.GetSubMesh(E3);
+ anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) return false;
+ aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[2] = (aVec[SMDSEntity_Node]-1)/2 + 2;
+ else
+ aNbNodes[2] = aVec[SMDSEntity_Node] + 2;
+ aNbNodes[3] = aNbNodes[1];
+ aNbNodes.resize(5);
+ nbSides = 4;
+ }
+ }
+ }
+ if ( nbEdgesInWire.front() == 4 ) { // exactly 4 edges
+ for(; edgeIt != edges.end(); edgeIt++) {
+ SMESH_subMesh * sm = aMesh.GetSubMesh( *edgeIt );
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) {
+ return false;
+ }
+ std::vector<int> aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
+ else
+ aNbNodes[nbSides] = aVec[SMDSEntity_Node] + 2;
+ nbSides++;
+ }
+ }
+ else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
+ list< TopoDS_Edge > sideEdges;
+ while ( !edges.empty()) {
+ sideEdges.clear();
+ sideEdges.splice( sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
+ bool sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide = SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.end(), edges, edges.begin());
+ }
+ if ( nbSides == 0 ) { // go backward from the first edge
+ sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide = SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.begin(), edges, --edges.end());
+ }
+ }
+ list<TopoDS_Edge>::iterator ite = sideEdges.begin();
+ aNbNodes[nbSides] = 1;
+ for(; ite!=sideEdges.end(); ite++) {
+ SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) {
+ return false;
+ }
+ std::vector<int> aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
+ else
+ aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
+ }
+ ++nbSides;
+ }
+ // issue 20222. Try to unite only edges shared by two same faces
+ if (nbSides < 4) {
+ nbSides = 0;
+ SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
+ while ( !edges.empty()) {
+ sideEdges.clear();
+ sideEdges.splice( sideEdges.end(), edges, edges.begin());
+ bool sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide =
+ SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() ) &&
+ twoEdgesMeatAtVertex( sideEdges.back(), edges.front(), aMesh );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.end(), edges, edges.begin());
+ }
+ if ( nbSides == 0 ) { // go backward from the first edge
+ sameSide = true;
+ while ( !edges.empty() && sameSide ) {
+ sameSide =
+ SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() ) &&
+ twoEdgesMeatAtVertex( sideEdges.front(), edges.back(), aMesh );
+ if ( sameSide )
+ sideEdges.splice( sideEdges.begin(), edges, --edges.end());
+ }
+ }
+ list<TopoDS_Edge>::iterator ite = sideEdges.begin();
+ aNbNodes[nbSides] = 1;
+ for(; ite!=sideEdges.end(); ite++) {
+ SMESH_subMesh * sm = aMesh.GetSubMesh( *ite );
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if(anIt==aResMap.end()) {
+ return false;
+ }
+ std::vector<int> aVec = (*anIt).second;
+ if(IsQuadratic)
+ aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
+ else
+ aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
+ }
+ ++nbSides;
+ }
+ }
+ }
+ if (nbSides != 4) {
+ if ( !nbSides )
+ nbSides = nbEdgesInWire.front();
+ error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
+ return false;
+ }
+
+ return true;
+}
+
+
//=============================================================================
/*!
* CheckAnd2Dcompute
}
// normalized 2d values on grid
- for (int i = 0; i < nbhoriz; i++)
- {
- for (int j = 0; j < nbvertic; j++)
- {
+ for (int i = 0; i < nbhoriz; i++) {
+ for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
// --- droite i cste : x = x0 + y(x1-x0)
double x0 = uv_e0[i].normParam; // bas - sud
gp_UV a2( uv_e2.back().u, uv_e2.back().v );
gp_UV a3( uv_e2.front().u, uv_e2.front().v );
- for (int i = 0; i < nbhoriz; i++)
- {
- for (int j = 0; j < nbvertic; j++)
- {
+ for (int i = 0; i < nbhoriz; i++) {
+ for (int j = 0; j < nbvertic; j++) {
int ij = j * nbhoriz + i;
double x = uv_grid[ij].x;
double y = uv_grid[ij].y;
return uv;
}
+//=======================================================================
+//function : CalcUV2
+//purpose : auxilary function for ComputeQuadPref
+//=======================================================================
+
+static gp_UV CalcUV2(double x, double y,
+ FaceQuadStruct* quad,
+ const gp_UV& a0, const gp_UV& a1,
+ const gp_UV& a2, const gp_UV& a3)
+{
+ const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0 );
+ const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
+ const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1 );
+ const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
+
+ //double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
+ //double y = y0 + x * (y1 - y0);
+
+ double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
+ double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
+ double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
+ double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
+
+ gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
+ gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
+ gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(param_t).XY();
+ gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(param_l).XY();
+
+ gp_UV uv = p0 * (1 - y) + p1 * x + p2 * y + p3 * (1 - x);
+
+ uv -= (1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3;
+
+ return uv;
+}
+
+
//=======================================================================
/*!
* Create only quandrangle faces
const TopoDS_Shape& aShape,
FaceQuadStruct* quad)
{
+ // Auxilary key in order to keep old variant
+ // of meshing after implementation new variant
+ // for bug 0016220 from Mantis.
+ bool OldVersion = false;
+
SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
const TopoDS_Face& F = TopoDS::Face(aShape);
Handle(Geom_Surface) S = BRep_Tool::Surface(F);
int addh = 0;
int addv = 0;
+ // ----------- Old version ---------------
// orientation of face and 3 main domain for future faces
// 0 top 1
// 1------------1
// 0------------0
// 0 bottom 1
+ // ----------- New version ---------------
+ // orientation of face and 3 main domain for future faces
+ // 0 top 1
+ // 1------------1
+ // | |____| |
+ // | / \ |
+ // | / C \ |
+ // left |/________\| rigth
+ // | |
+ // | |
+ // | |
+ // 0------------0
+ // 0 bottom 1
+
if(dh>dv) {
addv = (dh-dv)/2;
nbv = nbv + addv;
npl.Append(uv_el[i].normParam);
}
- // add some params to right and left after the first param
- // insert to right
- int dr = nbv - nr;
- double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
- for(i=1; i<=dr; i++) {
- npr.InsertAfter(1,npr.Value(2)-dpr);
- }
- // insert to left
- int dl = nbv - nl;
- dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
- for(i=1; i<=dl; i++) {
- npl.InsertAfter(1,npl.Value(2)-dpr);
+ int dl,dr;
+ if(OldVersion) {
+ // add some params to right and left after the first param
+ // insert to right
+ dr = nbv - nr;
+ double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
+ for(i=1; i<=dr; i++) {
+ npr.InsertAfter(1,npr.Value(2)-dpr);
+ }
+ // insert to left
+ dl = nbv - nl;
+ dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
+ for(i=1; i<=dl; i++) {
+ npl.InsertAfter(1,npl.Value(2)-dpr);
+ }
}
//cout<<"npb:";
//for(i=1; i<=npb.Length(); i++) {
TColgp_SequenceOfXY UVL;
TColgp_SequenceOfXY UVR;
- // step1: create faces for left domain
- StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
- // add left nodes
- for(j=1; j<=nl; j++)
- NodesL.SetValue(1,j,uv_el[j-1].node);
- if(dl>0) {
- // add top nodes
- for(i=1; i<=dl; i++)
- NodesL.SetValue(i+1,nl,uv_et[i].node);
+ if(OldVersion) {
+ // step1: create faces for left domain
+ StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
+ // add left nodes
+ for(j=1; j<=nl; j++)
+ NodesL.SetValue(1,j,uv_el[j-1].node);
+ if(dl>0) {
+ // add top nodes
+ for(i=1; i<=dl; i++)
+ NodesL.SetValue(i+1,nl,uv_et[i].node);
+ // create and add needed nodes
+ TColgp_SequenceOfXY UVtmp;
+ for(i=1; i<=dl; i++) {
+ double x0 = npt.Value(i+1);
+ double x1 = x0;
+ // diagonal node
+ double y0 = npl.Value(i+1);
+ double y1 = npr.Value(i+1);
+ gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesL.SetValue(i+1,1,N);
+ if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
+ // internal nodes
+ for(j=2; j<nl; j++) {
+ double y0 = npl.Value(dl+j);
+ double y1 = npr.Value(dl+j);
+ gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesL.SetValue(i+1,j,N);
+ if( i==dl ) UVtmp.Append(UV);
+ }
+ }
+ for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
+ UVL.Append(UVtmp.Value(i));
+ }
+ //cout<<"Dump NodesL:"<<endl;
+ //for(i=1; i<=dl+1; i++) {
+ // cout<<"i="<<i;
+ // for(j=1; j<=nl; j++) {
+ // cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
+ // }
+ // cout<<endl;
+ //}
+ // create faces
+ for(i=1; i<=dl; i++) {
+ for(j=1; j<nl; j++) {
+ if(WisF) {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
+ NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ else {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
+ NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ }
+ }
+ }
+ else {
+ // fill UVL using c2d
+ for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
+ UVL.Append( gp_UV ( uv_el[i].u, uv_el[i].v ));
+ }
+ }
+
+ // step2: create faces for right domain
+ StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
+ // add right nodes
+ for(j=1; j<=nr; j++)
+ NodesR.SetValue(1,j,uv_er[nr-j].node);
+ if(dr>0) {
+ // add top nodes
+ for(i=1; i<=dr; i++)
+ NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
+ // create and add needed nodes
+ TColgp_SequenceOfXY UVtmp;
+ for(i=1; i<=dr; i++) {
+ double x0 = npt.Value(nt-i);
+ double x1 = x0;
+ // diagonal node
+ double y0 = npl.Value(i+1);
+ double y1 = npr.Value(i+1);
+ gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesR.SetValue(i+1,nr,N);
+ if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
+ // internal nodes
+ for(j=2; j<nr; j++) {
+ double y0 = npl.Value(nbv-j+1);
+ double y1 = npr.Value(nbv-j+1);
+ gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesR.SetValue(i+1,j,N);
+ if( i==dr ) UVtmp.Prepend(UV);
+ }
+ }
+ for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
+ UVR.Append(UVtmp.Value(i));
+ }
+ // create faces
+ for(i=1; i<=dr; i++) {
+ for(j=1; j<nr; j++) {
+ if(WisF) {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
+ NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ else {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
+ NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ }
+ }
+ }
+ else {
+ // fill UVR using c2d
+ for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
+ UVR.Append( gp_UV( uv_er[i].u, uv_er[i].v ));
+ }
+ }
+
+ // step3: create faces for central domain
+ StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
+ // add first string using NodesL
+ for(i=1; i<=dl+1; i++)
+ NodesC.SetValue(1,i,NodesL(i,1));
+ for(i=2; i<=nl; i++)
+ NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
+ // add last string using NodesR
+ for(i=1; i<=dr+1; i++)
+ NodesC.SetValue(nb,i,NodesR(i,nr));
+ for(i=1; i<nr; i++)
+ NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
+ // add top nodes (last columns)
+ for(i=dl+2; i<nbh-dr; i++)
+ NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
+ // add bottom nodes (first columns)
+ for(i=2; i<nb; i++)
+ NodesC.SetValue(i,1,uv_eb[i-1].node);
+
// create and add needed nodes
- TColgp_SequenceOfXY UVtmp;
- for(i=1; i<=dl; i++) {
- double x0 = npt.Value(i+1);
+ // add linear layers
+ for(i=2; i<nb; i++) {
+ double x0 = npt.Value(dl+i);
double x1 = x0;
- // diagonal node
- double y0 = npl.Value(i+1);
- double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
- gp_Pnt P = S->Value(UV.X(),UV.Y());
- SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
- NodesL.SetValue(i+1,1,N);
- if(UVL.Length()<nbv-nnn-1) UVL.Append(UV);
- // internal nodes
- for(j=2; j<nl; j++) {
- double y0 = npl.Value(dl+j);
- double y1 = npr.Value(dl+j);
+ for(j=1; j<nnn; j++) {
+ double y0 = npl.Value(nbv-nnn+j);
+ double y1 = npr.Value(nbv-nnn+j);
gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
gp_Pnt P = S->Value(UV.X(),UV.Y());
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
- NodesL.SetValue(i+1,j,N);
- if( i==dl ) UVtmp.Append(UV);
+ NodesC.SetValue(i,nbv-nnn+j,N);
}
}
- for(i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
- UVL.Append(UVtmp.Value(i));
- }
- //cout<<"Dump NodesL:"<<endl;
- //for(i=1; i<=dl+1; i++) {
- // cout<<"i="<<i;
- // for(j=1; j<=nl; j++) {
- // cout<<" ("<<NodesL.Value(i,j)->X()<<","<<NodesL.Value(i,j)->Y()<<","<<NodesL.Value(i,j)->Z()<<")";
- // }
- // cout<<endl;
+ // add diagonal layers
+ //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
+ //cout<<"Dump UVL:"<<endl;
+ //for(i=1; i<=UVL.Length(); i++) {
+ // cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
//}
+ //cout<<endl;
+ for(i=1; i<nbv-nnn; i++) {
+ double du = UVR.Value(i).X() - UVL.Value(i).X();
+ double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
+ for(j=2; j<nb; j++) {
+ double u = UVL.Value(i).X() + du*npb.Value(j);
+ double v = UVL.Value(i).Y() + dv*npb.Value(j);
+ gp_Pnt P = S->Value(u,v);
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, u, v);
+ NodesC.SetValue(j,i+1,N);
+ }
+ }
// create faces
- for(i=1; i<=dl; i++) {
- for(j=1; j<nl; j++) {
+ for(i=1; i<nb; i++) {
+ for(j=1; j<nbv; j++) {
if(WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
- NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
+ myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
+ NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
else {
SMDS_MeshFace* F =
- myTool->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
- NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
+ myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
+ NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
}
}
}
- else {
- // fill UVL using c2d
- for(i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
- UVL.Append( gp_UV ( uv_el[i].u, uv_el[i].v ));
- }
- }
- // step2: create faces for right domain
- StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
- // add right nodes
- for(j=1; j<=nr; j++)
- NodesR.SetValue(1,j,uv_er[nr-j].node);
- if(dr>0) {
- // add top nodes
- for(i=1; i<=dr; i++)
- NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
- // create and add needed nodes
- TColgp_SequenceOfXY UVtmp;
- for(i=1; i<=dr; i++) {
- double x0 = npt.Value(nt-i);
- double x1 = x0;
- // diagonal node
- double y0 = npl.Value(i+1);
- double y1 = npr.Value(i+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
- gp_Pnt P = S->Value(UV.X(),UV.Y());
- SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
- NodesR.SetValue(i+1,nr,N);
- if(UVR.Length()<nbv-nnn-1) UVR.Append(UV);
- // internal nodes
- for(j=2; j<nr; j++) {
- double y0 = npl.Value(nbv-j+1);
- double y1 = npr.Value(nbv-j+1);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
- gp_Pnt P = S->Value(UV.X(),UV.Y());
+ else { // New version (!OldVersion)
+ // step1: create faces for bottom rectangle domain
+ StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
+ // fill UVL and UVR using c2d
+ for(j=0; j<nb; j++) {
+ NodesBRD.SetValue(j+1,1,uv_eb[j].node);
+ }
+ for(i=1; i<nnn-1; i++) {
+ NodesBRD.SetValue(1,i+1,uv_el[i].node);
+ NodesBRD.SetValue(nb,i+1,uv_er[i].node);
+ double du = uv_er[i].u - uv_el[i].u;
+ double dv = uv_er[i].v - uv_el[i].v;
+ for(j=2; j<nb; j++) {
+ double u = uv_el[i].u + du*npb.Value(j);
+ double v = uv_el[i].v + dv*npb.Value(j);
+ gp_Pnt P = S->Value(u,v);
SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
- NodesR.SetValue(i+1,j,N);
- if( i==dr ) UVtmp.Prepend(UV);
+ meshDS->SetNodeOnFace(N, geomFaceID, u, v);
+ NodesBRD.SetValue(j,i+1,N);
+
}
}
- for(i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
- UVR.Append(UVtmp.Value(i));
+ int nbf=0;
+ for(j=1; j<nnn-1; j++) {
+ for(i=1; i<nb; i++) {
+ nbf++;
+ if(WisF) {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
+ NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ else {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
+ NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ }
}
- // create faces
- for(i=1; i<=dr; i++) {
- for(j=1; j<nr; j++) {
+ int drl = abs(nr-nl);
+ // create faces for region C
+ StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
+ // add nodes from previous region
+ for(j=1; j<=nb; j++) {
+ NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
+ }
+ if( (drl+addv) > 0 ) {
+ int n1,n2;
+ if(nr>nl) {
+ n1 = 1;
+ n2 = drl + 1;
+ TColgp_SequenceOfXY UVtmp;
+ double drparam = npr.Value(nr) - npr.Value(nnn-1);
+ double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
+ double y0,y1;
+ for(i=1; i<=drl; i++) {
+ // add existed nodes from right edge
+ NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
+ //double dtparam = npt.Value(i+1);
+ y1 = npr.Value(nnn+i-1); // param on right edge
+ double dpar = (y1 - npr.Value(nnn-1))/drparam;
+ y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
+ double dy = y1 - y0;
+ for(j=1; j<nb; j++) {
+ double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
+ double y = y0 + dy*x;
+ gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesC.SetValue(j,i+1,N);
+ }
+ }
+ double dy0 = (1-y0)/(addv+1);
+ double dy1 = (1-y1)/(addv+1);
+ for(i=1; i<=addv; i++) {
+ double yy0 = y0 + dy0*i;
+ double yy1 = y1 + dy1*i;
+ double dyy = yy1 - yy0;
+ for(j=1; j<=nb; j++) {
+ double x = npt.Value(i+1+drl) +
+ npb.Value(j) * ( npt.Value(nt-i) - npt.Value(i+1+drl) );
+ double y = yy0 + dyy*x;
+ gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesC.SetValue(j,i+drl+1,N);
+ }
+ }
+ }
+ else { // nr<nl
+ n2 = 1;
+ n1 = drl + 1;
+ TColgp_SequenceOfXY UVtmp;
+ double dlparam = npl.Value(nl) - npl.Value(nnn-1);
+ double drparam = npr.Value(nnn) - npr.Value(nnn-1);
+ double y0 = npl.Value(nnn-1);
+ double y1 = npr.Value(nnn-1);
+ for(i=1; i<=drl; i++) {
+ // add existed nodes from right edge
+ NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
+ y0 = npl.Value(nnn+i-1); // param on left edge
+ double dpar = (y0 - npl.Value(nnn-1))/dlparam;
+ y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
+ double dy = y1 - y0;
+ for(j=2; j<=nb; j++) {
+ double x = npb.Value(j)*npt.Value(nt-i);
+ double y = y0 + dy*x;
+ gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesC.SetValue(j,i+1,N);
+ }
+ }
+ double dy0 = (1-y0)/(addv+1);
+ double dy1 = (1-y1)/(addv+1);
+ for(i=1; i<=addv; i++) {
+ double yy0 = y0 + dy0*i;
+ double yy1 = y1 + dy1*i;
+ double dyy = yy1 - yy0;
+ for(j=1; j<=nb; j++) {
+ double x = npt.Value(i+1) +
+ npb.Value(j) * ( npt.Value(nt-i-drl) - npt.Value(i+1) );
+ double y = yy0 + dyy*x;
+ gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
+ gp_Pnt P = S->Value(UV.X(),UV.Y());
+ SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
+ meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
+ NodesC.SetValue(j,i+drl+1,N);
+ }
+ }
+ }
+ // create faces
+ for(j=1; j<=drl+addv; j++) {
+ for(i=1; i<nb; i++) {
+ nbf++;
+ if(WisF) {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
+ NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ else {
+ SMDS_MeshFace* F =
+ myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
+ NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
+ }
+ }
+ } // end nr<nl
+
+ StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
+ for(i=1; i<=nt; i++) {
+ NodesLast.SetValue(i,2,uv_et[i-1].node);
+ }
+ int nnn=0;
+ for(i=n1; i<drl+addv+1; i++) {
+ nnn++;
+ NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
+ }
+ for(i=1; i<=nb; i++) {
+ nnn++;
+ NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
+ }
+ for(i=drl+addv; i>=n2; i--) {
+ nnn++;
+ NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
+ }
+ for(i=1; i<nt; i++) {
+ nbf++;
if(WisF) {
SMDS_MeshFace* F =
- myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
- NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
+ myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
+ NodesLast.Value(i+1,2), NodesLast.Value(i,2));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
else {
SMDS_MeshFace* F =
- myTool->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
- NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
+ myTool->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
+ NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
+ if(F) meshDS->SetMeshElementOnShape(F, geomFaceID);
}
}
+ } // if( (drl+addv) > 0 )
+
+ } // end new version implementation
+
+ bool isOk = true;
+ return isOk;
+}
+
+
+//=======================================================================
+/*!
+ * Evaluate only quandrangle faces
+ */
+//=======================================================================
+
+bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh,
+ const TopoDS_Shape& aShape,
+ std::vector<int>& aNbNodes,
+ MapShapeNbElems& aResMap,
+ bool IsQuadratic)
+{
+ // Auxilary key in order to keep old variant
+ // of meshing after implementation new variant
+ // for bug 0016220 from Mantis.
+ bool OldVersion = false;
+
+ const TopoDS_Face& F = TopoDS::Face(aShape);
+ Handle(Geom_Surface) S = BRep_Tool::Surface(F);
+
+ int nb = aNbNodes[0];
+ int nr = aNbNodes[1];
+ int nt = aNbNodes[2];
+ int nl = aNbNodes[3];
+ int dh = abs(nb-nt);
+ int dv = abs(nr-nl);
+
+ if( dh>=dv ) {
+ if( nt>nb ) {
+ // it is a base case => not shift
+ }
+ else {
+ // we have to shift on 2
+ nb = aNbNodes[2];
+ nr = aNbNodes[3];
+ nt = aNbNodes[0];
+ nl = aNbNodes[1];
}
}
else {
- // fill UVR using c2d
- for(i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
- UVR.Append( gp_UV( uv_er[i].u, uv_er[i].v ));
+ if( nr>nl ) {
+ // we have to shift quad on 1
+ nb = aNbNodes[3];
+ nr = aNbNodes[0];
+ nt = aNbNodes[1];
+ nl = aNbNodes[2];
+ }
+ else {
+ // we have to shift quad on 3
+ nb = aNbNodes[1];
+ nr = aNbNodes[2];
+ nt = aNbNodes[3];
+ nl = aNbNodes[0];
}
}
- // step3: create faces for central domain
- StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
- // add first string using NodesL
- for(i=1; i<=dl+1; i++)
- NodesC.SetValue(1,i,NodesL(i,1));
- for(i=2; i<=nl; i++)
- NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
- // add last string using NodesR
- for(i=1; i<=dr+1; i++)
- NodesC.SetValue(nb,i,NodesR(i,nr));
- for(i=1; i<nr; i++)
- NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
- // add top nodes (last columns)
- for(i=dl+2; i<nbh-dr; i++)
- NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
- // add bottom nodes (first columns)
- for(i=2; i<nb; i++)
- NodesC.SetValue(i,1,uv_eb[i-1].node);
-
- // create and add needed nodes
- // add linear layers
- for(i=2; i<nb; i++) {
- double x0 = npt.Value(dl+i);
- double x1 = x0;
- for(j=1; j<nnn; j++) {
- double y0 = npl.Value(nbv-nnn+j);
- double y1 = npr.Value(nbv-nnn+j);
- gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
- gp_Pnt P = S->Value(UV.X(),UV.Y());
- SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
- NodesC.SetValue(i,nbv-nnn+j,N);
+ dh = abs(nb-nt);
+ dv = abs(nr-nl);
+ int nbh = Max(nb,nt);
+ int nbv = Max(nr,nl);
+ int addh = 0;
+ int addv = 0;
+
+ if(dh>dv) {
+ addv = (dh-dv)/2;
+ nbv = nbv + addv;
+ }
+ else { // dv>=dh
+ addh = (dv-dh)/2;
+ nbh = nbh + addh;
+ }
+
+ int dl,dr;
+ if(OldVersion) {
+ // add some params to right and left after the first param
+ // insert to right
+ dr = nbv - nr;
+ // insert to left
+ dl = nbv - nl;
+ }
+
+ int nnn = Min(nr,nl);
+
+ int nbNodes = 0;
+ int nbFaces = 0;
+ if(OldVersion) {
+ // step1: create faces for left domain
+ if(dl>0) {
+ nbNodes += dl*(nl-1);
+ nbFaces += dl*(nl-1);
+ }
+ // step2: create faces for right domain
+ if(dr>0) {
+ nbNodes += dr*(nr-1);
+ nbFaces += dr*(nr-1);
}
+ // step3: create faces for central domain
+ nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
+ nbFaces += (nb-1)*(nbv-1);
}
- // add diagonal layers
- //cout<<"UVL.Length()="<<UVL.Length()<<" UVR.Length()="<<UVR.Length()<<endl;
- //cout<<"Dump UVL:"<<endl;
- //for(i=1; i<=UVL.Length(); i++) {
- // cout<<" ("<<UVL.Value(i).X()<<","<<UVL.Value(i).Y()<<")";
- //}
- //cout<<endl;
- for(i=1; i<nbv-nnn; i++) {
- double du = UVR.Value(i).X() - UVL.Value(i).X();
- double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
- for(j=2; j<nb; j++) {
- double u = UVL.Value(i).X() + du*npb.Value(j);
- double v = UVL.Value(i).Y() + dv*npb.Value(j);
- gp_Pnt P = S->Value(u,v);
- SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
- meshDS->SetNodeOnFace(N, geomFaceID, u, v);
- NodesC.SetValue(j,i+1,N);
+ else { // New version (!OldVersion)
+ nbNodes += (nnn-2)*(nb-2);
+ nbFaces += (nnn-2)*(nb-1);
+ int drl = abs(nr-nl);
+ nbNodes += drl*(nb-1) + addv*nb;
+ nbFaces += (drl+addv)*(nb-1) + (nt-1);
+ } // end new version implementation
+
+ std::vector<int> aVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
+ if(IsQuadratic) {
+ aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
+ aVec[SMDSEntity_Node] = nbNodes + nbFaces*4;
+ if( aNbNodes.size()==5 ) {
+ aVec[SMDSEntity_Quad_Triangle] = aNbNodes[3] - 1;
+ aVec[SMDSEntity_Quad_Quadrangle] = nbFaces - aNbNodes[3] + 1;
}
}
- // create faces
- for(i=1; i<nb; i++) {
- for(j=1; j<nbv; j++) {
- if(WisF) {
- SMDS_MeshFace* F =
- myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
- NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
- }
- else {
- SMDS_MeshFace* F =
- myTool->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
- NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
- meshDS->SetMeshElementOnShape(F, geomFaceID);
- }
+ else {
+ aVec[SMDSEntity_Node] = nbNodes;
+ aVec[SMDSEntity_Quadrangle] = nbFaces;
+ if( aNbNodes.size()==5 ) {
+ aVec[SMDSEntity_Triangle] = aNbNodes[3] - 1;
+ aVec[SMDSEntity_Quadrangle] = nbFaces - aNbNodes[3] + 1;
}
}
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aVec));
- bool isOk = true;
- return isOk;
+ return true;
}
+
+//=============================================================================
+/*! Split quadrangle in to 2 triangles by smallest diagonal
+ *
+ */
+//=============================================================================
+void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
+ int theFaceID,
+ const SMDS_MeshNode* theNode1,
+ const SMDS_MeshNode* theNode2,
+ const SMDS_MeshNode* theNode3,
+ const SMDS_MeshNode* theNode4)
+{
+ gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
+ gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
+ gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
+ gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
+ SMDS_MeshFace* face;
+ if(a.Distance(c) > b.Distance(d)){
+ face = myTool->AddFace(theNode2, theNode4 , theNode1);
+ if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
+ face = myTool->AddFace(theNode2, theNode3, theNode4);
+ if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
+
+ }
+ else{
+ face = myTool->AddFace(theNode1, theNode2 ,theNode3);
+ if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
+ face = myTool->AddFace(theNode1, theNode3, theNode4);
+ if(face) theMeshDS->SetMeshElementOnShape(face, theFaceID );
+ }
+}
+
+
