From: bph Date: Wed, 27 Jul 2011 14:44:10 +0000 (+0000) Subject: *** empty log message *** X-Git-Tag: EndWorkDidier~7 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=ebd8a5f5cfa1f199ecfd79a3fec7f30cee584d9e;p=tools%2Fmedcoupling.git *** empty log message *** --- diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index b390d445a..33c035227 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -22,6 +22,9 @@ #include "INTERPKERNELDefines.hxx" #include "InterpKernelException.hxx" +#if 1//dp +#include "VectorUtils.hxx" +#endif #include "NormalizedUnstructuredMesh.hxx" @@ -1097,5 +1100,4 @@ namespace INTERP_KERNEL } } - #endif diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx index 04d51338a..ecddf3c6e 100644 --- a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx @@ -81,13 +81,38 @@ namespace INTERP_KERNEL int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)); releaseArrays(); _split.splitTargetCell(targetCell,nbOfNodesT,_tetra); + for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) { double surface = 0.; + std::set listOfTetraFacesTreated; + + // calculate the coordinates of the nodes + const NumberingPolicy numPol=MyMeshType::My_numPol; + typename MyMeshType::MyConnType cellSrc = *iterCellS; + int cellSrcIdx = OTT::indFC(cellSrc); + NormalizedCellType normCellType=Intersector3D::_src_mesh.getTypeOfElement(cellSrcIdx); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + const MyMeshType& _src_mesh = Intersector3D::_src_mesh; + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(cellSrcIdx) : cellModelCell.getNumberOfNodes(); + int *polyNodes=new int[nbOfNodes4Type]; + double **polyCoords = new double*[nbOfNodes4Type]; + for(int i = 0;i<(int)nbOfNodes4Type;++i) + { + // we could store mapping local -> global numbers too, but not sure it is worth it + const int globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(*iterCellS), _src_mesh); + polyNodes[i]=globalNodeNum; + polyCoords[i] = (double*)_src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; + } + for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) - surface += (*iter)->intersectSourceFace(*iterCellS); + surface += (*iter)->intersectSourceFace(normCellType, nbOfNodes4Type, polyNodes, polyCoords, listOfTetraFacesTreated); if(surface!=0.) res[targetCell].insert(std::make_pair(OTT::indFC(*iterCellS), surface)); + listOfTetraFacesTreated.clear(); + delete[] polyNodes; + delete[] polyCoords; + } _split.releaseArrays(); } diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 3fbe1501d..e1bb11d0d 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -67,6 +67,23 @@ namespace INTERP_KERNEL return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2]; } + //dp TODO DP : à commenter + bool operator<(const TriangleFaceKey& key) const + { + for (int i = 0; i < 3; ++i) + { + if (_nodes[i] < key._nodes[i]) + { + return true; + } + else if (_nodes[i] > key._nodes[i]) + { + return false; + } + } + return false; + } + /** * Returns a hash value for the object, based on its three nodes. * This value is not unique for each face. @@ -173,7 +190,11 @@ namespace INTERP_KERNEL ~SplitterTetra(); double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0); - double intersectSourceFace(typename MyMeshType::MyConnType srcCell); + double intersectSourceFace(const NormalizedCellType polyType, + const int polyNodesNbr, + const int *const polyNodes, + const double *const *const polyCoords, + std::set& listOfTetraFacesTreated); double intersectTetra(const double** tetraCorners); @@ -189,11 +210,17 @@ namespace INTERP_KERNEL // member functions inline void createAffineTransform(const double** corners); inline void checkIsOutside(const double* pt, bool* isOutside) const; - inline void checkIsOutsideSurface(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const; //dp à supprimer + inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const; //dp à supprimer inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum); inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key); - + + inline bool isFacesCoplanar(const double *const plane_normal, const double plane_constant, + const double *const *const coordsFace); + inline double calculateIntersectionSurfaceOfCoplanarTriangles(const double *const normal, + const double *const P_1, const double *const P_2, const double *const P_3, + const double *const P_4, const double *const P_5, const double *const P_6, + double dim_caracteristic, double precision); /// disallow copying SplitterTetra(const SplitterTetra& t); @@ -259,20 +286,21 @@ namespace INTERP_KERNEL isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); } +#if 1//dp //dp à supprimer template - inline void SplitterTetra::checkIsOutsideSurface(const double* pt, bool* isOutside, const double errTol) const + inline void SplitterTetra::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const { - isOutside[0] = isOutside[0] && (pt[0] < -errTol); - isOutside[1] = isOutside[1] && (pt[0] > (1.0 + errTol)); - isOutside[2] = isOutside[2] && (pt[1] < -errTol); - isOutside[3] = isOutside[3] && (pt[1] > (1.0 + errTol)); - isOutside[4] = isOutside[4] && (pt[2] < -errTol); - isOutside[5] = isOutside[5] && (pt[2] > (1.0 + errTol)); - isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol); - isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol)); + isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol); + isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol)); + isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol); + isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol)); + isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol); + isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol)); + isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol); + isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol)); } - +#endif /** * Calculates the transformed node with a given global node number. diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 071b7a816..c3748f506 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -28,6 +28,7 @@ #include "CellModel.hxx" #include "Log.hxx" #include "UnitTetraIntersectionBary.hxx" +#include "VolSurfFormulae.hxx" #include #include @@ -258,6 +259,7 @@ namespace INTERP_KERNEL faceNodes=new int[faceModel.getNumberOfNodes()]; cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes); } + // intersect a son with the unit tetra switch(faceType) { @@ -373,6 +375,183 @@ namespace INTERP_KERNEL return std::fabs(1.0 / _t->determinant() * totalVolume) ; } +#if 1//dp TODO DP : à commenter + template + double SplitterTetra::calculateIntersectionSurfaceOfCoplanarTriangles(const double *const plane_normal, + const double *const P_1, const double *const P_2, const double *const P_3, + const double *const P_4, const double *const P_5, const double *const P_6, + double dim_caracteristic, double precision) + { + typedef typename MyMeshType::MyConnType ConnType; + typedef double Vect2[2]; + typedef double Vect3[3]; + typedef double Triangle2[3][2]; + + const double *const tri0[3] = {P_1, P_2, P_3}; + const double *const tri1[3] = {P_4, P_5, P_6}; + + // Plane of the first triangle defined by the normal of the triangle and the constant + // Project triangles onto coordinate plane most aligned with plane normal + int maxNormal = 0; + double fmax = std::abs(plane_normal[0]); + double absMax = std::abs(plane_normal[1]); + if (absMax > fmax) + { + maxNormal = 1; + fmax = absMax; + } + absMax = std::abs(plane_normal[2]); + if (absMax > fmax) + { + maxNormal = 2; + } + + Triangle2 projTri0, projTri1; + int i; + + if (maxNormal == 0) + { + // Project onto yz-plane. + for (i = 0; i < 3; ++i) + { + projTri0[i][0] = tri0[i][1]; + projTri0[i][1] = tri0[i][2]; + projTri1[i][0] = tri1[i][1]; + projTri1[i][1] = tri1[i][2]; + } + } + else if (maxNormal == 1) + { + // Project onto xz-plane. + for (i = 0; i < 3; ++i) + { + projTri0[i][0] = tri0[i][0]; + projTri0[i][1] = tri0[i][2]; + projTri1[i][0] = tri1[i][0]; + projTri1[i][1] = tri1[i][2]; + } + } + else + { + // Project onto xy-plane. + for (i = 0; i < 3; ++i) + { + projTri0[i][0] = tri0[i][0]; + projTri0[i][1] = tri0[i][1]; + projTri1[i][0] = tri1[i][0]; + projTri1[i][1] = tri1[i][1]; + } + } + + // 2D triangle intersection routines require counterclockwise ordering. + Vect2 save; + Vect2 edge0; + Vect2 edge1; + subtract(projTri0[1], projTri0[0], edge0); + subtract(projTri0[2], projTri0[0], edge1); + if ((edge0[0] * edge1[1] - edge0[1] - edge1[0]) < (double) 0.) + { + // Triangle is clockwise, reorder it. + for (int i = 0; i < 2; ++i) + { + save[i] = projTri0[1][i]; + projTri0[1][i] = projTri0[2][i]; + projTri0[2][i] = save[i]; + } + } + + subtract(projTri1[1], projTri1[0], edge0); + subtract(projTri1[2], projTri1[0], edge1); + if ((edge0[0] * edge1[1] - edge0[1] - edge1[0]) < (double) 0.) + { + // Triangle is clockwise, reorder it. + for (int i = 0; i < 2; ++i) + { + save[i] = projTri1[1][i]; + projTri1[1][i] = projTri1[2][i]; + projTri1[2][i] = save[i]; + } + } + + std::vector inter2; + intersec_de_triangle(projTri0[0], projTri0[1], projTri0[2], + projTri1[0], projTri1[1], projTri1[2], + inter2, + dim_caracteristic, precision); + ConnType nb_inter=((ConnType)inter2.size())/2; + if(nb_inter >3) inter2=reconstruct_polygon(inter2); + double surface = 0.; + double normal[3]; + for(ConnType i = 1; i(&inter2[0],&inter2[2*i],&inter2[2*(i+1)],normal); + surface +=0.5*fabs(normal[0]); + } + return surface; +#if 0//dp + int nb_inter = inter2.size(); + if (nb_inter > 0) + { + inter3.resize(3 * nb_inter / 2); + // Map 2D intersections back to the 3D triangle space. + if (maxNormal == 0) + { + double invNX = ((double) 1.) / normal[0]; + for (i = 0; i < nb_inter; i++) + { + inter3[3 * i + 1] = inter2[2 * i]; + inter3[3 * i + 2] = inter2[2 * i + 1]; + inter3[3 * i] = invNX * (constant - normal[1] * inter3[3 * i + 1] - normal[2] * inter3[3 * i + 2]); + } + } + else if (maxNormal == 1) + { + double invNY = ((double) 1.) / normal[1]; + for (i = 0; i < nb_inter; i++) + { + inter3[3 * i] = inter2[2 * i]; + inter3[3 * i + 2] = inter2[2 * i + 1]; + inter3[3 * i + 1] = invNY * (constant - normal[0] * inter3[3 * i] - normal[2] * inter3[3 * i + 2]); + } + } + else + { + double invNZ = ((double) 1.) / normal[2]; + for (i = 0; i < nb_inter; i++) + { + inter3[3 * i] = inter2[2 * i]; + inter3[3 * i + 1] = inter2[2 * i + 1]; + inter3[3 * i + 2] = invNZ * (constant - normal[0] * inter3[3 * i] - normal[1] * inter3[3 * i + 1]); + } + } + } +#endif + } + + template + bool SplitterTetra::isFacesCoplanar(const double *const plane_normal, const double plane_constant, + const double *const *const coordsFace) + { + // Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test. + // For faces not left + int zero = 0; + for (int i = 0; i < 3; ++i) + { + const double distance = dot(plane_normal, coordsFace[i]) - plane_constant; + if (epsilonEqual(distance, 0.)) + { + zero++; + } + } + if (zero == 3) + return true; + else + return false; + } + +#endif + + // TODO DP : adapter les commentaires suivants. _volume => _surface ? /** * Calculates the volume of intersection of an element in the source mesh and the target element. @@ -388,27 +567,18 @@ namespace INTERP_KERNEL * @param element global number of the source element in C mode. */ template - double SplitterTetra::intersectSourceFace(typename MyMeshType::MyConnType element) + double SplitterTetra::intersectSourceFace(const NormalizedCellType polyType, + const int polyNodesNbr, + const int *const polyNodes, + const double *const *const polyCoords, + std::set& listOfTetraFacesTreated) { typedef typename MyMeshType::MyConnType ConnType; - const NumberingPolicy numPol=MyMeshType::My_numPol; - 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(); double totalVolume = 0.0; -#if 1//dp - // calculate the coordinates of the nodes - int *cellNodes=new int[nbOfNodes4Type]; - for(int i = 0;i<(int)nbOfNodes4Type;++i) - { - // we could store mapping local -> global numbers too, but not sure it is worth it - const int globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(element), _src_mesh); - cellNodes[i]=globalNodeNum; - //const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; - } +#if 0 // Is src element coplanar with one of the tetra faces ? double srcNormal[3]; @@ -420,17 +590,19 @@ namespace INTERP_KERNEL } calculateNormalForTria(points[0],points[1],points[2], srcNormal); #else - const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; - _nodes[globalNodeNum] = node; - calculateNormalForPolyg(_nodes, nbOfNodes4Type, srcNormal); + calculateNormalForPolyg((const double **)coordsPoly, nbOfNodes4Type, srcNormal); #endif double faceNormal[3]; double crossNormals[3]; - for (int iFace = 0; iFace < 4; ++iFace) + const int tetraFacesNodesConn[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,2,3}}; + for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace) { - int decal = iFace * 3; - calculateNormalForTria(_coords + decal, _coords + decal + 1, _coords + decal + 2, faceNormal); + const int *const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace]; + const double *const coordsTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM; + const double *const coordsTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM; + const double *const coordsTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM; + calculateNormalForTria(coordsTriNode1, coordsTriNode2, coordsTriNode3, faceNormal); cross(srcNormal, faceNormal, crossNormals); if (epsilonEqual(norm(crossNormals), 0.)) { @@ -440,8 +612,8 @@ namespace INTERP_KERNEL for (int iTri = 0; iTri < nbTria; ++iTri) { std::vector inter; - INTERP_KERNEL::intersec_de_triangle(_nodes[cellNodes[0]], _nodes[cellNodes[1 + iTri]], _nodes[cellNodes[2 + iTri]], - _coords + decal,_coords + decal + 1,_coords + decal + 2, + INTERP_KERNEL::intersec_de_triangle(coordsPoly[0], coordsPoly[1 + iTri], coordsPoly[2 + iTri], + coordsTriNode1, coordsTriNode2, coordsTriNode3, inter, 1., //dp inter, PlanarIntersector::_dim_caracteristic, DEFAULT_ABS_TOL); //dp PlanarIntersector::_precision); @@ -453,7 +625,6 @@ namespace INTERP_KERNEL totalVolume +=0.5*fabs(area[0]); } } - break; } } #endif @@ -469,20 +640,22 @@ namespace INTERP_KERNEL // halfspace filtering bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isStrictlyOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetStrictlyOutside = false; bool isTargetOutside = false; // calculate the coordinates of the nodes #if 0//dp - int *cellNodes=new int[nbOfNodes4Type]; + int *polyNodes=new int[polyNodesNbr]; #endif - for(int i = 0;i<(int)nbOfNodes4Type;++i) + for(int i = 0;i<(int)polyNodesNbr;++i) { #if 0//dp // we could store mapping local -> global numbers too, but not sure it is worth it const int globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(element), _src_mesh); - cellNodes[i]=globalNodeNum; + polyNodes[i]=globalNodeNum; #else - const int globalNodeNum = cellNodes[i]; + const int globalNodeNum = polyNodes[i]; #endif if(_nodes.find(globalNodeNum) == _nodes.end()) { @@ -492,118 +665,189 @@ namespace INTERP_KERNEL calculateNode(globalNodeNum); } - checkIsOutsideSurface(_nodes[globalNodeNum], isOutside); + checkIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside); + checkIsOutside(_nodes[globalNodeNum], isOutside); } // halfspace filtering check // NB : might not be beneficial for caching of triangles for(int i = 0; i < 8; ++i) { - if(isOutside[i]) + if(isStrictlyOutside[i]) + { + isTargetStrictlyOutside = true; + break; + } + else if (isOutside[i]) { isTargetOutside = true; } } - if (!isTargetOutside) + if (!isTargetStrictlyOutside) { - // intersect a son with the unit tetra - switch (normCellType) + + if (isTargetOutside) { - case NORM_TRI3: + // Faces are parallel + const int tetraFacesNodesConn[4][3] = { + { 0, 1, 2 }, + { 0, 2, 3 }, + { 0, 3, 1 }, + { 1, 2, 3 } }; + double plane_normal[3]; + for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace) { - // create the face key - TriangleFaceKey key = TriangleFaceKey(cellNodes[0], cellNodes[1], cellNodes[2]); - - // calculate the triangle if needed - if (_volumes.find(key) == _volumes.end()) - { - TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1]], _nodes[cellNodes[2]]); - calculateSurface(tri, key); - totalVolume += _volumes[key]; - } - else + const int * const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace]; + TriangleFaceKey key = TriangleFaceKey(_conn[tetraFaceNodesConn[0]], + _conn[tetraFaceNodesConn[1]], + _conn[tetraFaceNodesConn[2]]); + if (listOfTetraFacesTreated.find(key) == listOfTetraFacesTreated.end()) { - // count negative as face has reversed orientation - totalVolume -= _volumes[key]; + listOfTetraFacesTreated.insert(key); + const double * const coordsTetraTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM; + const double * const coordsTetraTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM; + const double * const coordsTetraTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM; + calculateNormalForTria(coordsTetraTriNode1, coordsTetraTriNode2, coordsTetraTriNode3, plane_normal); + const double normOfTetraTriNormal = norm(plane_normal); + if (epsilonEqual(normOfTetraTriNormal, 0.)) + { + for (int i = 0; i < 3; ++i) + { + plane_normal[i] = 0.; + } + } + else + { + const double invNormOfTetraTriNormal = 1. / normOfTetraTriNormal; + for (int i = 0; i < 3; ++i) + { + plane_normal[i] *= invNormOfTetraTriNormal; + } + } + double plane_constant = dot(plane_normal, coordsTetraTriNode1); + if (isFacesCoplanar(plane_normal, plane_constant, polyCoords)) + { + int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles + for (int iTri = 0; iTri < nbrPolyTri; ++iTri) + { + totalVolume += calculateIntersectionSurfaceOfCoplanarTriangles(plane_normal, + polyCoords[0], + polyCoords[1 + iTri], + polyCoords[2 + iTri], + coordsTetraTriNode1, + coordsTetraTriNode2, + coordsTetraTriNode3, + 1., //dp PlanarIntersector::_dim_caracteristic, + DEFAULT_ABS_TOL); //dp PlanarIntersector::_precision); + } + } } } - break; - - case NORM_QUAD4: //dp pas d'intérêt vis-à-vis du case suivant (à regrouper les deux cases) - - // simple triangulation of faces along a diagonal : - // - // 2 ------ 3 - // | / | - // | / | - // | / | - // | / | - // | / | - // | / | - // 1 ------ 4 - // - //? not sure if this always works - { - // calculate the triangles if needed - - // local nodes 1, 2, 3 - TriangleFaceKey key1 = TriangleFaceKey(cellNodes[0], cellNodes[1], cellNodes[2]); - if (_volumes.find(key1) == _volumes.end()) - { - TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1]], _nodes[cellNodes[2]]); - calculateSurface(tri, key1); - totalVolume += _volumes[key1]; - } - else - { - // count negative as face has reversed orientation - totalVolume -= _volumes[key1]; - } - - // local nodes 1, 3, 4 - TriangleFaceKey key2 = TriangleFaceKey(cellNodes[0], cellNodes[2], cellNodes[3]); - if (_volumes.find(key2) == _volumes.end()) - { - TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[2]], _nodes[cellNodes[3]]); - calculateSurface(tri, key2); - totalVolume += _volumes[key2]; - } - else + } + else + { + // intersect a son with the unit tetra + switch (polyType) + { + case NORM_TRI3: { - // count negative as face has reversed orientation - totalVolume -= _volumes[key2]; - } - } - break; + // create the face key + TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]); - case NORM_POLYGON: - { - int nbTria = nbOfNodes4Type - 2; // split polygon into nbTria triangles - for (int iTri = 0; iTri < nbTria; ++iTri) - { - TriangleFaceKey key = TriangleFaceKey(cellNodes[0], cellNodes[1 + iTri], cellNodes[2 + iTri]); + // calculate the triangle if needed if (_volumes.find(key) == _volumes.end()) { - TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1 + iTri]], _nodes[cellNodes[2 + iTri]]); + TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]); calculateSurface(tri, key); totalVolume += _volumes[key]; } else { + // count negative as face has reversed orientation totalVolume -= _volumes[key]; } } - } - break; + break; + + case NORM_QUAD4: //dp pas d'intérêt vis-à-vis du case suivant (à regrouper les deux cases) + + // simple triangulation of faces along a diagonal : + // + // 2 ------ 3 + // | / | + // | / | + // | / | + // | / | + // | / | + // | / | + // 1 ------ 4 + // + //? not sure if this always works + { + // calculate the triangles if needed + + // local nodes 1, 2, 3 + TriangleFaceKey key1 = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]); + if (_volumes.find(key1) == _volumes.end()) + { + TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]); + calculateSurface(tri, key1); + totalVolume += _volumes[key1]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key1]; + } + + // local nodes 1, 3, 4 + TriangleFaceKey key2 = TriangleFaceKey(polyNodes[0], polyNodes[2], polyNodes[3]); + if (_volumes.find(key2) == _volumes.end()) + { + TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[2]], _nodes[polyNodes[3]]); + calculateSurface(tri, key2); + totalVolume += _volumes[key2]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key2]; + } + } + break; + + case NORM_POLYGON: + { + int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles + for (int iTri = 0; iTri < nbrPolyTri; ++iTri) + { + TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]); + if (_volumes.find(key) == _volumes.end()) + { + TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1 + iTri]], + _nodes[polyNodes[2 + iTri]]); + calculateSurface(tri, key); + totalVolume += _volumes[key]; + } + else + { + totalVolume -= _volumes[key]; + } + } + } + break; + + default: + std::cout + << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." + << std::endl; + assert(false); + } - default: - std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl; - assert(false); } } - delete [] cellNodes; // reset if it is very small to keep the matrix sparse // is this a good idea? if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT)) diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 5b4025418..7b27d5f5f 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -238,7 +238,9 @@ namespace ParaMEDMEM CPPUNIT_TEST( test2DInterpP1P0PL_2 ); CPPUNIT_TEST( test2DInterpP1P1_1 ); CPPUNIT_TEST( test2DInterpP1P1PL_1 ); +#endif CPPUNIT_TEST( test3DSurfInterpP0P0_1 ); +#if 0 CPPUNIT_TEST( test3DSurfInterpP0P0PL_1 ); CPPUNIT_TEST( test3DSurfInterpP0P1_1 ); CPPUNIT_TEST( test3DSurfInterpP0P1PL_1 ); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx index 456343d41..ad76c84f3 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx @@ -1179,22 +1179,28 @@ double MEDCouplingBasicsTest::sumAll(const std::vector< std::map >& #if 1 MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1() { - double sourceCoords[42]={-12., 6., 10., -12.,10., 6., -16.,10., 10., - -20., 0., 0., -12., 0., 0., -12., 0., -4., -20.,0.,-4., - -20., 0., 10., -12., 0., 10., -20.,10., 10., - -25., 5., -5., 5., 5., -5., 5., 5., 25., -25.,5.,25.}; - int sourceConn[14]={0,1,2, 3,4,5,6, 7,8,9, 10,11,12,13}; + double sourceCoords[63]={-12., 6., 10., -12.,10., 6., -16.,10. , 10., + -20., 0., 0., -12., 0., 0., -12., 0. , -4., -20.,0.,-4., + -20., 0., 10., -12., 0., 10., -20.,10. , 10., + -25., 5., -5., 5., 5., -5., 5., 5. , 25., -25.,5.,25., + -20., 0., 16., -18., 0., 16., -20., 2.5, 16., + -25., 0., -5., 5., 0., -5., 5., 0. , 25., -25.,0.,25. + }; + int sourceConn[25]={0,1,2, 3,4,5,6, 7,8,9, 10,11,12,13, 14,15,16, 3,4,8,7, 17,18,19,20}; MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New(); sourceMesh->setMeshDimension(2); - sourceMesh->allocateCells(1); - //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn); - //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+3); - //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn+7); + sourceMesh->allocateCells(7); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+3); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn+7); sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+10); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn+14); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+17); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+21); sourceMesh->finishInsertingCells(); DataArrayDouble *myCoords=DataArrayDouble::New(); - myCoords->alloc(14,3); - std::copy(sourceCoords,sourceCoords+42,myCoords->getPointer()); + myCoords->alloc(21,3); + std::copy(sourceCoords,sourceCoords+63,myCoords->getPointer()); sourceMesh->setCoords(myCoords); myCoords->decrRef(); return sourceMesh; @@ -1206,8 +1212,8 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh_1() int sourceConn[12]={4,5,7,8, 0,3,2,1,4,7,6,5}; MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New(); sourceMesh->setMeshDimension(3); - sourceMesh->allocateCells(1); - //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,sourceConn); + sourceMesh->allocateCells(2); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,sourceConn); sourceMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,sourceConn + 4); sourceMesh->finishInsertingCells(); DataArrayDouble *myCoords=DataArrayDouble::New(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx index 02d97c2d1..16b053689 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -2304,25 +2304,31 @@ void MEDCouplingBasicsTest::test3D2DInterpP0P0_1() INTERP_KERNEL::Interpolation3D2D myInterpolator; myInterpolator.setPrecision(1e-12); std::vector > res; - INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::GENERAL_48 }; - //INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 }; + //INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::GENERAL_48 }; + INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 }; for ( int i = 0; i < 1; ++i ) { myInterpolator.setSplittingPolicy( sp[i] ); res.clear(); myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); - CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); - - //CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][0],1e-12); - //CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][1],1e-12); - //CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[0][2],1e-12); - //CPPUNIT_ASSERT_DOUBLES_EQUAL(8. ,res[0][3],1e-12); - - //CPPUNIT_ASSERT_DOUBLES_EQUAL(8.*sqrt(3.),res[1][0],1e-12); - //CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[1][1],1e-12); - //CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[1][2],1e-12); - CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[0][0],1e-12); + CPPUNIT_ASSERT_EQUAL(2,(int)res.size()); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8. ,res[0][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.5 ,res[0][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][5],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(32. ,res[0][6],1e-12); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(8.*sqrt(3.),res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[1][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[1][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[1][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][5],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][6],1e-12); //CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7); }