X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2FINTERP_KERNEL%2FSplitterTetra.txx;h=ce0136d99360f1b9c708b70eadd3b703e77ac738;hb=b219559763498c4bd10c730cd3d2c62b1eed45db;hp=9312b5fa6eb8d31dfc81259d5e712d3eab7825d2;hpb=5312c32a71d5e4f50cd8160e1de1a58903d5b428;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 9312b5fa6..ce0136d99 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -1,9 +1,9 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D +// Copyright (C) 2007-2019 CEA/DEN, EDF R&D // // 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. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -84,6 +84,32 @@ namespace INTERP_KERNEL // create the affine transform _t=new TetraAffineTransform(_coords); } + + /** + * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a + * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish + * if a partial computation of a triangle has already been performed (to increase performance). + * + * The \a srcMesh can contain polyhedron cells. + * + * + * Constructor creating object from the four corners of the tetrahedron. + * + * \param [in] srcMesh mesh containing the source elements + * \param [in] tetraCorners array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z). + */ + template + SplitterTetra::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn): _t(0),_src_mesh(srcMesh) + { + if(!conn) + { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; } + else + { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; } + _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5]; + _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11]; + // create the affine transform + _t=new TetraAffineTransform(_coords); + } /** * Destructor @@ -95,7 +121,7 @@ namespace INTERP_KERNEL SplitterTetra::~SplitterTetra() { delete _t; - for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter) + for(typename HashMap< ConnType, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter) delete[] iter->second; } @@ -165,17 +191,17 @@ namespace INTERP_KERNEL // get type of cell NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(element)); const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); - unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(element)) : cellModelCell.getNumberOfNodes(); + ConnType nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(element)) : cellModelCell.getNumberOfNodes(); // halfspace filtering bool isOutside[8] = {true, true, true, true, true, true, true, true}; bool isTargetOutside = false; // calculate the coordinates of the nodes - int *cellNodes=new int[nbOfNodes4Type]; - for(int i = 0;i<(int)nbOfNodes4Type;++i) + ConnType *cellNodes=new ConnType[nbOfNodes4Type]; + for(ConnType i = 0;i global numbers too, but not sure it is worth it - const int globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(element), _src_mesh); + const ConnType globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(element), _src_mesh); cellNodes[i]=globalNodeNum; if(_nodes.find(globalNodeNum) == _nodes.end()) { @@ -184,8 +210,7 @@ namespace INTERP_KERNEL //std::cout << std::endl << "*** " << globalNodeNum << std::endl; calculateNode(globalNodeNum); } - - checkIsOutside(_nodes[globalNodeNum], isOutside); + CheckIsOutside(_nodes[globalNodeNum], isOutside); } // halfspace filtering check @@ -207,19 +232,19 @@ namespace INTERP_KERNEL // get nb of sons of a cell const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]); - const int rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ]; + const ConnType rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ]; unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes); for(unsigned ii = 0 ; ii < nbOfSons; ++ii) { // get sons connectivity NormalizedCellType faceType; - int *faceNodes, nbFaceNodes=-1; + ConnType *faceNodes, nbFaceNodes=-1; if ( cellModelCell.isDynamic() ) { - faceNodes=new int[nbOfNodes4Type]; + faceNodes=new ConnType[nbOfNodes4Type]; nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType); - for ( int i = 0; i < nbFaceNodes; ++i ) + for ( ConnType i = 0; i < nbFaceNodes; ++i ) faceNodes[i] = OTT::coo2C(faceNodes[i]); } else @@ -227,7 +252,8 @@ namespace INTERP_KERNEL faceType = cellModelCell.getSonType(ii); const CellModel& faceModel=CellModel::GetCellModel(faceType); assert(faceModel.getDimension() == 2); - faceNodes=new int[faceModel.getNumberOfNodes()]; + nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii); + faceNodes = new ConnType[nbFaceNodes]; cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes); } // intersect a son with the unit tetra @@ -300,8 +326,8 @@ namespace INTERP_KERNEL case NORM_POLYGON: { - int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles - for ( int iTri = 0; iTri < nbTria; ++iTri ) + ConnType nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles + for ( ConnType iTri = 0; iTri < nbTria; ++iTri ) { TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]); if(_volumes.find(key) == _volumes.end()) @@ -368,7 +394,6 @@ namespace INTERP_KERNEL { typedef typename MyMeshType::MyConnType ConnType; typedef double Vect2[2]; - typedef double Vect3[3]; typedef double Triangle2[3][2]; const double *const tri0[3] = {p1, p2, p3}; @@ -557,7 +582,6 @@ namespace INTERP_KERNEL * @param polyNodesNbr number of the nodes of the polygon source face * @param polyNodes numbers of the nodes of the polygon source face * @param polyCoords coordinates of the nodes of the polygon source face - * @param polyCoords coordinates of the nodes of the polygon source face * @param dimCaracteristic characteristic size of the meshes containing the triangles * @param precision precision for double float data used for comparison * @param listOfTetraFacesTreated list of tetra faces treated @@ -565,16 +589,14 @@ namespace INTERP_KERNEL */ template double SplitterTetra::intersectSourceFace(const NormalizedCellType polyType, - const int polyNodesNbr, - const int *const polyNodes, + const ConnType polyNodesNbr, + const ConnType *const polyNodes, const double *const *const polyCoords, const double dimCaracteristic, const double precision, std::multiset& listOfTetraFacesTreated, std::set& listOfTetraFacesColinear) { - typedef typename MyMeshType::MyConnType ConnType; - double totalSurface = 0.0; // check if we have planar tetra element @@ -592,16 +614,16 @@ namespace INTERP_KERNEL bool isTargetOutside = false; // calculate the coordinates of the nodes - for(int i = 0;i<(int)polyNodesNbr;++i) + for(ConnType i = 0;iapply(nodes[i], tetraCorners[i]); - checkIsOutside(nodes[i], isOutside); + CheckIsOutside(nodes[i], isOutside); } // halfspace filtering check @@ -844,7 +866,7 @@ namespace INTERP_KERNEL if(!isTargetOutside) { const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4); - int cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3]; + ConnType cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3]; for(unsigned ii = 0 ; ii < 4 ; ++ii) { @@ -887,9 +909,9 @@ namespace INTERP_KERNEL void SplitterTetra2::releaseArrays() { // free potential sub-mesh nodes that have been allocated - typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. - if((int)_nodes.size()>=/*8*/nbOfNodesT) + if(_nodes.size()>=/*8*/_node_ids.size()) { + typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast(_node_ids.size()); std::vector::iterator iter = _nodes.begin() + /*8*/nbOfNodesT; while(iter != _nodes.end()) { @@ -899,6 +921,47 @@ namespace INTERP_KERNEL } _nodes.clear(); } + + /*! + * \param [in] targetCell in C mode. + * \param [out] tetra is the output result tetra containers. + */ + template + void SplitterTetra2::splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra* >& tetra) + { + typedef typename MyMeshTypeT::MyConnType TConnType; + const TConnType *refConn(_target_mesh.getConnectivityPtr()); + const TConnType *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]); + INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell)); + std::vector tetrasNodalConn; + std::vector addCoords; + const double *coords(_target_mesh.getCoordinatesPtr()); + SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords); + std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras); + double tmp[12]; + typename MyMeshTypeS::MyConnType tmp2[4]; + for(std::size_t i=0;i=0) + { + tmp[j*3+0]=coords[3*cellId+0]; + tmp[j*3+1]=coords[3*cellId+1]; + tmp[j*3+2]=coords[3*cellId+2]; + } + else + { + tmp[j*3+0]=addCoords[3*(-cellId-1)+0]; + tmp[j*3+1]=addCoords[3*(-cellId-1)+1]; + tmp[j*3+2]=addCoords[3*(-cellId-1)+2]; + } + } + tetra[i]=new SplitterTetra(_src_mesh,tmp,tmp2); + } + } /*! * @param targetCell in C mode. @@ -918,7 +981,7 @@ namespace INTERP_KERNEL _node_ids.resize(8); tetra.reserve(1); const double *nodes[4]; - int conn[4]; + ConnType conn[4]; for(int node = 0; node < 4 ; ++node) { nodes[node]=getCoordsOfNode2(node, OTT::indFC(targetCell),_target_mesh,conn[node]); @@ -999,7 +1062,7 @@ namespace INTERP_KERNEL for(int i = 0; i < 5; ++i) { const double* nodes[4]; - int conn[4]; + typename MyMeshTypeS::MyConnType conn[4]; for(int j = 0; j < 4; ++j) { conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ]; @@ -1023,7 +1086,7 @@ namespace INTERP_KERNEL for(int i = 0; i < 6; ++i) { const double* nodes[4]; - int conn[4]; + typename MyMeshTypeS::MyConnType conn[4]; for(int j = 0; j < 4; ++j) { conn[j] = subZone[SPLIT_NODES_6[4*i+j]]; @@ -1046,11 +1109,11 @@ namespace INTERP_KERNEL { // The two nodes of the original mesh cell used in each tetrahedron. // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2) - // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes + // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes // nodes to use for tetrahedron const double* nodes[4]; - int conn[4]; + typename MyMeshTypeS::MyConnType conn[4]; // get the cell center conn[0] = 14; nodes[0] = getCoordsOfSubNode(conn[0]); @@ -1085,25 +1148,10 @@ namespace INTERP_KERNEL */ template void SplitterTetra2::calculateGeneral48Tetra(typename std::vector< SplitterTetra* >& tetra) - { - // Define 8 hexahedral subzones as in Grandy, p449 - // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell - // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes - static const int subZones[64] = - { - 0,8,21,12,9,20,26,22, - 8,1,13,21,20,10,23,26, - 12,21,16,3,22,26,25,17, - 21,13,2,16,26,23,18,25, - 9,20,26,22,4,11,24,14, - 20,10,23,26,11,5,15,24, - 22,26,25,17,14,24,19,7, - 26,23,18,25,24,15,6,19 - }; - + { for(int i = 0; i < 8; ++i) { - sixSplit(&subZones[8*i],tetra); + sixSplit(GENERAL_48_SUBZONES+8*i,tetra); } } @@ -1125,7 +1173,7 @@ namespace INTERP_KERNEL // create tetrahedra const double* nodes[4]; - int conn[4]; + typename MyMeshTypeS::MyConnType conn[4]; for(int i = 0; i < 2; ++i) { for(int j = 0; j < 4; ++j) @@ -1151,22 +1199,24 @@ namespace INTERP_KERNEL // get type of cell and nb of cell nodes NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT::indFC(targetCell)); const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); - unsigned nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)) : cellModelCell.getNumberOfNodes(); + ConnType nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)) : cellModelCell.getNumberOfNodes(); // get nb of cell sons (faces) const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]); - const int rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ]; + const ConnType rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ]; unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes); // indices of nodes of a son - static std::vector allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1 - while ( allNodeIndices.size() < nbOfCellNodes ) - allNodeIndices.push_back( allNodeIndices.size() ); - std::vector classicFaceNodes(4); - int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0]; + static std::vector allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1 + while ( allNodeIndices.size() < (std::size_t)nbOfCellNodes ) + allNodeIndices.push_back( static_cast(allNodeIndices.size()) ); + std::vector classicFaceNodes(4); + if(cellModelCell.isQuadratic()) + throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !"); + ConnType* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0]; // nodes of tetrahedron - int conn[4]; + typename MyMeshTypeS::MyConnType conn[4]; const double* nodes[4]; nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter @@ -1211,7 +1261,7 @@ namespace INTERP_KERNEL { // retrieve real mesh nodes - typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. _node_ids.resize(8); + typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast(_node_ids.size());// Issue 0020634. _node_ids.resize(8); for(int node = 0; node < nbOfNodesT ; ++node) { // calculate only normal nodes @@ -1238,33 +1288,6 @@ namespace INTERP_KERNEL case GENERAL_48: { - // Each sub-node is the barycenter of two other nodes. - // For the edges, these lie on the original mesh. - // For the faces, these are the edge sub-nodes. - // For the cell these are two face sub-nodes. - static const int GENERAL_48_SUB_NODES[38] = - { - 0,1, // sub-node 9 (edge) - 0,4, // sub-node 10 (edge) - 1,5, // sub-node 11 (edge) - 4,5, // sub-node 12 (edge) - 0,3, // sub-node 13 (edge) - 1,2, // sub-node 14 (edge) - 4,7, // sub-node 15 (edge) - 5,6, // sub-node 16 (edge) - 2,3, // sub-node 17 (edge) - 3,7, // sub-node 18 (edge) - 2,6, // sub-node 19 (edge) - 6,7, // sub-node 20 (edge) - 8,11, // sub-node 21 (face) - 12,13, // sub-node 22 (face) - 9,17, // sub-node 23 (face) - 10,18, // sub-node 24 (face) - 14,15, // sub-node 25 (face) - 16,19, // sub-node 26 (face) - 20,25 // sub-node 27 (cell) - }; - for(int i = 0; i < 19; ++i) { double* barycenter = new double[3];