X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Hexa_3D.cxx;h=a3ed4f11a7142cca9be57bc944b2e49de2af9d57;hp=6c631d0db0e1ac3e6f9af6030a3deeb407f2b610;hb=2387bfa403855b82751bf9f122295b1fc6923a18;hpb=c3bf92bd87b770fd81631a3853f7f5bb1ac6a4e8 diff --git a/src/StdMeshers/StdMeshers_Hexa_3D.cxx b/src/StdMeshers/StdMeshers_Hexa_3D.cxx index 6c631d0db..a3ed4f11a 100644 --- a/src/StdMeshers/StdMeshers_Hexa_3D.cxx +++ b/src/StdMeshers/StdMeshers_Hexa_3D.cxx @@ -32,10 +32,13 @@ using namespace std; #include "StdMeshers_Quadrangle_2D.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_subMesh.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" #include "SMDS_FacePosition.hxx" +#include "SMDS_VolumeTool.hxx" +#include "SMDS_VolumeOfNodes.hxx" #include #include @@ -49,10 +52,16 @@ using namespace std; #include #include #include +#include #include "utilities.h" #include "Utils_ExceptHandlers.hxx" +//modified by NIZNHY-PKV Wed Nov 17 15:31:58 2004 f +#include "StdMeshers_Penta_3D.hxx" + +static bool ComputePentahedralMesh(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape); +//modified by NIZNHY-PKV Wed Nov 17 15:32:00 2004 t //============================================================================= /*! @@ -81,6 +90,8 @@ StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D() { MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D"); + for (int i = 0; i < 6; i++) + StdMeshers_Quadrangle_2D::QuadDelete(_quads[i]); } //============================================================================= @@ -94,7 +105,7 @@ bool StdMeshers_Hexa_3D::CheckHypothesis const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { - MESSAGE("StdMeshers_Hexa_3D::CheckHypothesis"); + //MESSAGE("StdMeshers_Hexa_3D::CheckHypothesis"); bool isOk = true; aStatus = SMESH_Hypothesis::HYP_OK; @@ -104,6 +115,37 @@ bool StdMeshers_Hexa_3D::CheckHypothesis return isOk; } +//======================================================================= +//function : findIJ +//purpose : return i,j of the node +//======================================================================= + +static bool findIJ (const SMDS_MeshNode* node, const FaceQuadStruct * quad, int& I, int& J) +{ + I = J = 0; + const SMDS_FacePosition* fpos = + static_cast(node->GetPosition().get()); + if ( ! fpos ) return false; + gp_Pnt2d uv( fpos->GetUParameter(), fpos->GetVParameter() ); + + double minDist = DBL_MAX; + int nbhoriz = Min(quad->nbPts[0], quad->nbPts[2]); + int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]); + for (int i = 1; i < nbhoriz - 1; i++) { + for (int j = 1; j < nbvertic - 1; j++) { + int ij = j * nbhoriz + i; + gp_Pnt2d uv2( quad->uv_grid[ij].u, quad->uv_grid[ij].v ); + double dist = uv.SquareDistance( uv2 ); + if ( dist < minDist ) { + minDist = dist; + I = i; + J = j; + } + } + } + return true; +} + //============================================================================= /*! * Hexahedron mesh on hexaedron like form @@ -122,15 +164,14 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, { Unexpect aCatch(SalomeException); MESSAGE("StdMeshers_Hexa_3D::Compute"); - - bool isOk = false; + //bool isOk = false; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); + //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape); //const SMESHDS_SubMesh *& subMeshDS = theSubMesh->GetSubMeshDS(); // 0. - shape and face mesh verification // 0.1 - shape must be a solid (or a shell) with 6 faces - MESSAGE("---"); + //MESSAGE("---"); vector < SMESH_subMesh * >meshFaces; for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) @@ -142,56 +183,82 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, if (meshFaces.size() != 6) { SCRUTE(meshFaces.size()); - ASSERT(0); +// ASSERT(0); return false; } // 0.2 - is each face meshed with Quadrangle_2D? (so, with a wire of 4 edges) - MESSAGE("---"); + //MESSAGE("---"); for (int i = 0; i < 6; i++) { - TopoDS_Shape aShape = meshFaces[i]->GetSubShape(); - SMESH_Algo *algo = _gen->GetAlgo(aMesh, aShape); - string algoName = algo->GetName(); - if (algoName != "Quadrangle_2D") - { - // *** delete _quads - SCRUTE(algoName); - ASSERT(0); - return false; - } - StdMeshers_Quadrangle_2D *quadAlgo = - dynamic_cast < StdMeshers_Quadrangle_2D * >(algo); - ASSERT(quadAlgo); - try - { - _quads[i] = quadAlgo->CheckAnd2Dcompute(aMesh, aShape); - // *** to delete after usage - } - catch(SALOME_Exception & S_ex) - { - // *** delete _quads - // *** throw exception - ASSERT(0); - } + TopoDS_Shape aFace = meshFaces[i]->GetSubShape(); + SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace); + string algoName = algo->GetName(); + bool isAllQuad = false; + if (algoName == "Quadrangle_2D") { + SMESHDS_SubMesh * sm = meshDS->MeshElements( aFace ); + if ( sm ) { + isAllQuad = true; + SMDS_ElemIteratorPtr eIt = sm->GetElements(); + while ( isAllQuad && eIt->more() ) + isAllQuad = ( eIt->next()->NbNodes() == 4 ); + } + } + if ( ! isAllQuad ) { + //modified by NIZNHY-PKV Wed Nov 17 15:31:37 2004 f + bool bIsOk; + // + bIsOk=ComputePentahedralMesh(aMesh, aShape); + if (bIsOk) { + return true; + } + //modified by NIZNHY-PKV Wed Nov 17 15:31:42 2004 t + SCRUTE(algoName); + // ASSERT(0); + return false; + } + StdMeshers_Quadrangle_2D *quadAlgo = + dynamic_cast < StdMeshers_Quadrangle_2D * >(algo); + ASSERT(quadAlgo); + try + { + _quads[i] = quadAlgo->CheckAnd2Dcompute(aMesh, aFace); + // *** to delete after usage + } + catch(SALOME_Exception & S_ex) + { + // *** delete _quads + // *** throw exception + // ASSERT(0); + return false; + } + + // 0.2.1 - number of points on the opposite edges must be the same + if (_quads[i]->nbPts[0] != _quads[i]->nbPts[2] || + _quads[i]->nbPts[1] != _quads[i]->nbPts[3]) + { + MESSAGE("different number of points on the opposite edges of face " << i); + // ASSERT(0); + return false; + } } // 1. - identify faces and vertices of the "cube" // 1.1 - ancestor maps vertex->edges in the cube - MESSAGE("---"); + //MESSAGE("---"); TopTools_IndexedDataMapOfShapeListOfShape MS; TopExp::MapShapesAndAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, MS); // 1.2 - first face is choosen as face Y=0 of the unit cube - MESSAGE("---"); + //MESSAGE("---"); const TopoDS_Shape & aFace = meshFaces[0]->GetSubShape(); const TopoDS_Face & F = TopoDS::Face(aFace); // 1.3 - identify the 4 vertices of the face Y=0: V000, V100, V101, V001 - MESSAGE("---"); + //MESSAGE("---"); int i = 0; TopoDS_Edge E = _quads[0]->edge[i]; //edge will be Y=0,Z=0 on unit cube @@ -237,7 +304,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, // - find edge X=1, Z=0 (ancestor of V100 not in face Y=0) // - find edge X=1, Z=1 (ancestor of V101 not in face Y=0) // - find edge X=0, Z=1 (ancestor of V001 not in face Y=0) - MESSAGE("---"); + //MESSAGE("---"); TopoDS_Edge E_0Y0 = EdgeNotInFace(aMesh, aShape, F, _cube.V000, MS); ASSERT(!E_0Y0.IsNull()); @@ -252,7 +319,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, ASSERT(!E_0Y1.IsNull()); // 1.5 - identify the 4 vertices in face Y=1: V010, V110, V111, V011 - MESSAGE("---"); + //MESSAGE("---"); TopExp::Vertices(E_0Y0, VFirst, VLast); if (VFirst.IsSame(_cube.V000)) @@ -279,7 +346,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, _cube.V011 = VFirst; // 1.6 - find remaining faces given 4 vertices - MESSAGE("---"); + //MESSAGE("---"); _indY0 = 0; _cube.quad_Y0 = _quads[_indY0]; @@ -304,7 +371,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, _cube.V100, _cube.V101, _cube.V110, _cube.V111); _cube.quad_X1 = _quads[_indX1]; - MESSAGE("---"); + //MESSAGE("---"); // 1.7 - get convertion coefs from face 2D normalized to 3D normalized @@ -330,18 +397,17 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, // 1.8 - create a 3D structure for normalized values - MESSAGE("---"); - int nbx = _cube.quad_Y0->nbPts[0]; - int nby = _cube.quad_Y0->nbPts[1]; - int nbz; - if (cx0.a1 != 0) - nbz = _cube.quad_X0->nbPts[1]; - else - nbz = _cube.quad_X0->nbPts[0]; - //SCRUTE(nbx); - //SCRUTE(nby); - //SCRUTE(nbz); - int nbxyz = nbx * nby * nbz; + //MESSAGE("---"); + int nbx = _cube.quad_Z0->nbPts[0]; + if (cz0.a1 == 0.) nbx = _cube.quad_Z0->nbPts[1]; + + int nby = _cube.quad_X0->nbPts[0]; + if (cx0.a1 == 0.) nby = _cube.quad_X0->nbPts[1]; + + int nbz = _cube.quad_Y0->nbPts[0]; + if (cy0.a1 != 0.) nbz = _cube.quad_Y0->nbPts[1]; + + int i1, j1, nbxyz = nbx * nby * nbz; Point3DStruct *np = new Point3DStruct[nbxyz]; // 1.9 - store node indexes of faces @@ -360,12 +426,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition* fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -396,12 +457,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition* fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -432,12 +488,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition * fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -468,12 +519,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition* fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -504,12 +550,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition * fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -540,12 +581,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, while(itf->more()) { const SMDS_MeshNode * node = itf->next(); - const SMDS_FacePosition* fpos = - static_cast(node->GetPosition().get()); - double ri = fpos->GetUParameter(); - double rj = fpos->GetVParameter(); - int i1 = int (ri); - int j1 = int (rj); + findIJ( node, quad, i1, j1 ); int ij1 = j1 * nbdown + i1; quad->uv_grid[ij1].node = node; } @@ -570,17 +606,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, // - compute the point 3D // - store the point 3D in SMESHDS, store its ID in 3D structure - TopoDS_Shell aShell; - TopExp_Explorer exp(aShape, TopAbs_SHELL); - if (exp.More()) - { - aShell = TopoDS::Shell(exp.Current()); - } - else - { - MESSAGE("no shell..."); - ASSERT(0); - } + int shapeID = meshDS->ShapeToIndex( aShape ); Pt3 p000, p001, p010, p011, p100, p101, p110, p111; Pt3 px00, px01, px10, px11; @@ -654,17 +680,39 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMDS_MeshNode * node = meshDS->AddNode(X[0], X[1], X[2]); np[ijk].node = node; - //meshDS->SetNodeInVolume(node, TopoDS::Solid(aShape)); - meshDS->SetNodeInVolume(node, aShell); + meshDS->SetNodeInVolume(node, shapeID); } } } + // find orientation of furute volumes according to MED convention + vector< bool > forward( nbx * nby ); + SMDS_VolumeTool vTool; + for (int i = 0; i < nbx - 1; i++) + for (int j = 0; j < nby - 1; j++) + { + int n1 = j * nbx + i; + int n2 = j * nbx + i + 1; + int n3 = (j + 1) * nbx + i + 1; + int n4 = (j + 1) * nbx + i; + int n5 = nbx * nby + j * nbx + i; + int n6 = nbx * nby + j * nbx + i + 1; + int n7 = nbx * nby + (j + 1) * nbx + i + 1; + int n8 = nbx * nby + (j + 1) * nbx + i; + + SMDS_VolumeOfNodes tmpVol (np[n1].node,np[n2].node,np[n3].node,np[n4].node, + np[n5].node,np[n6].node,np[n7].node,np[n8].node); + vTool.Set( &tmpVol ); + forward[ n1 ] = vTool.IsForward(); + } + //2.1 - for each node of the cube (less 3 *1 Faces): // - store hexahedron in SMESHDS MESSAGE("Storing hexahedron into the DS"); for (int i = 0; i < nbx - 1; i++) for (int j = 0; j < nby - 1; j++) + { + bool isForw = forward.at( j * nbx + i ); for (int k = 0; k < nbz - 1; k++) { int n1 = k * nbx * nby + j * nbx + i; @@ -676,47 +724,31 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, int n7 = (k + 1) * nbx * nby + (j + 1) * nbx + i + 1; int n8 = (k + 1) * nbx * nby + (j + 1) * nbx + i; -// MESSAGE(" "<AddVolume(np[n1].node, - np[n2].node, - np[n3].node, - np[n4].node, - np[n5].node, - np[n6].node, - np[n7].node, - np[n8].node); - ; - meshDS->SetMeshElementOnShape(elt, aShell); - - // *** 5 tetrahedres ... verifier orientations, - // mettre en coherence &vec quadrangles-> triangles - // choisir afficher 1 parmi edges, face et volumes -// int tetra1 = meshDS->AddVolume(np[n1].nodeId, -// np[n2].nodeId, -// np[n4].nodeId, -// np[n5].nodeId); -// int tetra2 = meshDS->AddVolume(np[n2].nodeId, -// np[n3].nodeId, -// np[n4].nodeId, -// np[n7].nodeId); -// int tetra3 = meshDS->AddVolume(np[n5].nodeId, -// np[n6].nodeId, -// np[n7].nodeId, -// np[n2].nodeId); -// int tetra4 = meshDS->AddVolume(np[n5].nodeId, -// np[n7].nodeId, -// np[n8].nodeId, -// np[n4].nodeId); -// int tetra5 = meshDS->AddVolume(np[n5].nodeId, -// np[n7].nodeId, -// np[n2].nodeId, -// np[n4].nodeId); - - } - - MESSAGE("End of StdMeshers_Hexa_3D::Compute()"); + SMDS_MeshVolume * elt; + if ( isForw ) + elt = meshDS->AddVolume(np[n1].node, + np[n2].node, + np[n3].node, + np[n4].node, + np[n5].node, + np[n6].node, + np[n7].node, + np[n8].node); + else + elt = meshDS->AddVolume(np[n1].node, + np[n4].node, + np[n3].node, + np[n2].node, + np[n5].node, + np[n8].node, + np[n7].node, + np[n6].node); + + meshDS->SetMeshElementOnShape(elt, shapeID); + } + } + if ( np ) delete [] np; + //MESSAGE("End of StdMeshers_Hexa_3D::Compute()"); return true; } @@ -750,7 +782,7 @@ int StdMeshers_Hexa_3D::GetFaceIndex(SMESH_Mesh & aMesh, const TopoDS_Vertex & V1, const TopoDS_Vertex & V2, const TopoDS_Vertex & V3) { - MESSAGE("StdMeshers_Hexa_3D::GetFaceIndex"); + //MESSAGE("StdMeshers_Hexa_3D::GetFaceIndex"); int faceIndex = -1; for (int i = 1; i < 6; i++) { @@ -771,7 +803,7 @@ int StdMeshers_Hexa_3D::GetFaceIndex(SMESH_Mesh & aMesh, } } ASSERT(faceIndex > 0); - SCRUTE(faceIndex); + //SCRUTE(faceIndex); return faceIndex; } @@ -788,13 +820,13 @@ TopoDS_Edge const TopoDS_Vertex & aVertex, const TopTools_IndexedDataMapOfShapeListOfShape & MS) { - MESSAGE("StdMeshers_Hexa_3D::EdgeNotInFace"); + //MESSAGE("StdMeshers_Hexa_3D::EdgeNotInFace"); TopTools_IndexedDataMapOfShapeListOfShape MF; TopExp::MapShapesAndAncestors(aFace, TopAbs_VERTEX, TopAbs_EDGE, MF); const TopTools_ListOfShape & ancestorsInSolid = MS.FindFromKey(aVertex); const TopTools_ListOfShape & ancestorsInFace = MF.FindFromKey(aVertex); - SCRUTE(ancestorsInSolid.Extent()); - SCRUTE(ancestorsInFace.Extent()); +// SCRUTE(ancestorsInSolid.Extent()); +// SCRUTE(ancestorsInFace.Extent()); ASSERT(ancestorsInSolid.Extent() == 6); // 6 (edges doublees) ASSERT(ancestorsInFace.Extent() == 2); @@ -836,7 +868,7 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad, const TopoDS_Vertex & V1, const TopoDS_Vertex & V2, const TopoDS_Vertex & V3, Conv2DStruct & conv) { - MESSAGE("StdMeshers_Hexa_3D::GetConv2DCoefs"); +// MESSAGE("StdMeshers_Hexa_3D::GetConv2DCoefs"); const TopoDS_Face & F = TopoDS::Face(aShape); TopoDS_Edge E = quad.edge[0]; double f, l; @@ -936,8 +968,8 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad, b2 = -1; c2 = 1; // 1-y } - MESSAGE("X = " << c1 << "+ " << a1 << "*x + " << b1 << "*y"); - MESSAGE("Y = " << c2 << "+ " << a2 << "*x + " << b2 << "*y"); +// MESSAGE("X = " << c1 << "+ " << a1 << "*x + " << b1 << "*y"); +// MESSAGE("Y = " << c2 << "+ " << a2 << "*x + " << b2 << "*y"); conv.a1 = a1; conv.b1 = b1; conv.c1 = c1; @@ -955,8 +987,8 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad, conv.jb = int (b2); conv.jc = int (c2 * a2 * a2) * (nbdown - 1) + int (c2 * b2 * b2) * (nbright - 1); - MESSAGE("I " << conv.ia << " " << conv.ib << " " << conv.ic); - MESSAGE("J " << conv.ja << " " << conv.jb << " " << conv.jc); +// MESSAGE("I " << conv.ia << " " << conv.ib << " " << conv.ic); +// MESSAGE("J " << conv.ja << " " << conv.jb << " " << conv.jc); } //============================================================================= @@ -1002,3 +1034,35 @@ istream & operator >>(istream & load, StdMeshers_Hexa_3D & hyp) { return hyp.LoadFrom( load ); } + +//modified by NIZNHY-PKV Wed Nov 17 15:34:13 2004 f +/////////////////////////////////////////////////////////////////////////////// +//ZZ +//#include + +//======================================================================= +//function : ComputePentahedralMesh +//purpose : +//======================================================================= +bool ComputePentahedralMesh(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) +{ + //printf(" ComputePentahedralMesh HERE\n"); + // + bool bOK; + //int iErr; + StdMeshers_Penta_3D anAlgo; + // + bOK=anAlgo.Compute(aMesh, aShape); + /* + iErr=anAlgo.ErrorStatus(); + + if (iErr) { + printf(" *** Error# %d\n", iErr); + } + else { + printf(" *** No errors# %d\n", iErr); + } + */ + return bOK; +} +