#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 <TopExp.hxx>
#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
#include <Geom2d_Curve.hxx>
#include <Handle_Geom2d_Curve.hxx>
#include <Handle_Geom_Curve.hxx>
+#include <gp_Pnt2d.hxx>
#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
//=============================================================================
/*!
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]);
}
//=============================================================================
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;
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<const SMDS_FacePosition*>(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
{
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())
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
// - 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());
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))
_cube.V011 = VFirst;
// 1.6 - find remaining faces given 4 vertices
- MESSAGE("---");
+ //MESSAGE("---");
_indY0 = 0;
_cube.quad_Y0 = _quads[_indY0];
_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
// 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
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition* fpos =
- static_cast<const SMDS_FacePosition*>(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;
}
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition* fpos =
- static_cast<const SMDS_FacePosition*>(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;
}
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition * fpos =
- static_cast<const SMDS_FacePosition*>(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;
}
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition* fpos =
- static_cast<const SMDS_FacePosition *>(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;
}
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition * fpos =
- static_cast<const SMDS_FacePosition*>(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;
}
while(itf->more())
{
const SMDS_MeshNode * node = itf->next();
- const SMDS_FacePosition* fpos =
- static_cast<const SMDS_FacePosition*>(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;
}
// - 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;
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;
int n7 = (k + 1) * nbx * nby + (j + 1) * nbx + i + 1;
int n8 = (k + 1) * nbx * nby + (j + 1) * nbx + i;
-// MESSAGE(" "<<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<" "<<n7<<" "<<n8);
- //MESSAGE(" "<<np[n1].nodeId<<" "<<np[n2].nodeId<<" "<<np[n3].nodeId<<" "<<np[n4].nodeId<<" "<<np[n5].nodeId<<" "<<np[n6].nodeId<<" "<<np[n7].nodeId<<" "<<np[n8].nodeId);
-
- SMDS_MeshVolume * 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);
- ;
- 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;
}
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++)
{
}
}
ASSERT(faceIndex > 0);
- SCRUTE(faceIndex);
+ //SCRUTE(faceIndex);
return faceIndex;
}
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);
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;
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;
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);
}
//=============================================================================
{
return hyp.LoadFrom( load );
}
+
+//modified by NIZNHY-PKV Wed Nov 17 15:34:13 2004 f
+///////////////////////////////////////////////////////////////////////////////
+//ZZ
+//#include <stdio.h>
+
+//=======================================================================
+//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;
+}
+
