From 204a1120bdd573522acb77a7282dcec5133e9d03 Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 2 Aug 2013 13:12:19 +0000 Subject: [PATCH] Start debugging 3D interpolation error on OCTA12 in target mesh --- src/INTERP_KERNEL/CMakeLists.txt | 1 + .../Polyhedron3D2DIntersectorP0P0.txx | 80 +++---- src/INTERP_KERNEL/SplitterTetra.cxx | 216 ++++++++++++++++++ src/INTERP_KERNEL/SplitterTetra.hxx | 72 +++++- src/INTERP_KERNEL/SplitterTetra.txx | 47 +--- 5 files changed, 328 insertions(+), 88 deletions(-) create mode 100644 src/INTERP_KERNEL/SplitterTetra.cxx diff --git a/src/INTERP_KERNEL/CMakeLists.txt b/src/INTERP_KERNEL/CMakeLists.txt index 0c29d1e3b..7d2e4a422 100644 --- a/src/INTERP_KERNEL/CMakeLists.txt +++ b/src/INTERP_KERNEL/CMakeLists.txt @@ -38,6 +38,7 @@ SET(interpkernel_SOURCES InterpKernelCellSimplify.cxx InterpKernelMatrixTools.cxx VolSurfUser.cxx + SplitterTetra.cxx Bases/InterpKernelException.cxx Geometric2D/InterpKernelGeo2DAbstractEdge.cxx Geometric2D/InterpKernelGeo2DBounds.cxx diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx index ce3d1a10c..4c1e4fb9e 100644 --- a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx @@ -81,7 +81,6 @@ namespace INTERP_KERNEL * * @param targetCell in C mode. * @param srcCells in C mode. - * */ template void Polyhedron3D2DIntersectorP0P0::intersectCells(ConnType targetCell, @@ -111,7 +110,7 @@ namespace INTERP_KERNEL { // 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; + polyNodes[i] = globalNodeNum; polyCoords[i] = const_cast(src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum); } @@ -125,46 +124,43 @@ namespace INTERP_KERNEL listOfTetraFacesTreated, listOfTetraFacesColinear); - if(surface!=0.) { - - matrix[targetCell].insert(std::make_pair(cellSrcIdx, surface)); - - bool isSrcFaceColinearWithFaceOfTetraTargetCell = false; - std::set::iterator iter; - for (iter = listOfTetraFacesColinear.begin(); iter != listOfTetraFacesColinear.end(); ++iter) - { - if (listOfTetraFacesTreated.count(*iter) != 1) - { - isSrcFaceColinearWithFaceOfTetraTargetCell = false; - break; - } - else - { - isSrcFaceColinearWithFaceOfTetraTargetCell = true; - } - } - - if (isSrcFaceColinearWithFaceOfTetraTargetCell) - { - DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(cellSrcIdx); - if (intersectFacesIter != _intersect_faces.end()) - { - intersectFacesIter->second.insert(targetCell); - } - else - { - std::set targetCellSet; - targetCellSet.insert(targetCell); - _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet)); - } - - } - - } - - delete[] polyNodes; - delete[] polyCoords; - + if(surface!=0.) + { + + matrix[targetCell].insert(std::make_pair(cellSrcIdx, surface)); + + bool isSrcFaceColinearWithFaceOfTetraTargetCell = false; + std::set::iterator iter; + for (iter = listOfTetraFacesColinear.begin(); iter != listOfTetraFacesColinear.end(); ++iter) + { + if (listOfTetraFacesTreated.count(*iter) != 1) + { + isSrcFaceColinearWithFaceOfTetraTargetCell = false; + break; + } + else + { + isSrcFaceColinearWithFaceOfTetraTargetCell = true; + } + } + + if (isSrcFaceColinearWithFaceOfTetraTargetCell) + { + DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(cellSrcIdx); + if (intersectFacesIter != _intersect_faces.end()) + { + intersectFacesIter->second.insert(targetCell); + } + else + { + std::set targetCellSet; + targetCellSet.insert(targetCell); + _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet)); + } + } + } + delete [] polyNodes; + delete [] polyCoords; } _split.releaseArrays(); } diff --git a/src/INTERP_KERNEL/SplitterTetra.cxx b/src/INTERP_KERNEL/SplitterTetra.cxx new file mode 100644 index 000000000..b5e4999af --- /dev/null +++ b/src/INTERP_KERNEL/SplitterTetra.cxx @@ -0,0 +1,216 @@ +// Copyright (C) 2007-2013 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. +// +// 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 +// +// Author : Anthony Geay (CEA/DEN) + +#include "SplitterTetra.hxx" + +namespace INTERP_KERNEL +{ + + void SplitHexa8IntoTetras(SplittingPolicy policy, const int *nodalConnBg, const int *nodalConnEnd, const double *coords, + std::vector& tetrasNodalConn, std::vector& addCoords) throw(INTERP_KERNEL::Exception) + { + if(std::distance(nodalConnBg,nodalConnEnd)!=8) + throw INTERP_KERNEL::Exception("SplitHexa8IntoTetras : input hexa do not have 8 nodes !"); + switch(policy) + { + case PLANAR_FACE_5: + { + tetrasNodalConn.resize(20); + int *conn(&tetrasNodalConn[0]); + conn[0]=nodalConnBg[SPLIT_NODES_5_WO[0]]; conn[1]=nodalConnBg[SPLIT_NODES_5_WO[1]]; conn[2]=nodalConnBg[SPLIT_NODES_5_WO[2]]; conn[3]=nodalConnBg[SPLIT_NODES_5_WO[3]]; + conn[4]=nodalConnBg[SPLIT_NODES_5_WO[4]]; conn[5]=nodalConnBg[SPLIT_NODES_5_WO[5]]; conn[6]=nodalConnBg[SPLIT_NODES_5_WO[6]]; conn[7]=nodalConnBg[SPLIT_NODES_5_WO[7]]; + conn[8]=nodalConnBg[SPLIT_NODES_5_WO[8]]; conn[9]=nodalConnBg[SPLIT_NODES_5_WO[9]]; conn[10]=nodalConnBg[SPLIT_NODES_5_WO[10]]; conn[11]=nodalConnBg[SPLIT_NODES_5_WO[11]]; + conn[12]=nodalConnBg[SPLIT_NODES_5_WO[12]]; conn[13]=nodalConnBg[SPLIT_NODES_5_WO[13]]; conn[14]=nodalConnBg[SPLIT_NODES_5_WO[14]]; conn[15]=nodalConnBg[SPLIT_NODES_5_WO[15]]; + conn[16]=nodalConnBg[SPLIT_NODES_5_WO[16]]; conn[17]=nodalConnBg[SPLIT_NODES_5_WO[17]]; conn[18]=nodalConnBg[SPLIT_NODES_5_WO[18]]; conn[19]=nodalConnBg[SPLIT_NODES_5_WO[19]]; + return ; + } + case PLANAR_FACE_6: + { + tetrasNodalConn.resize(24); + int *conn(&tetrasNodalConn[0]); + conn[0]=nodalConnBg[SPLIT_NODES_6_WO[0]]; conn[1]=nodalConnBg[SPLIT_NODES_6_WO[1]]; conn[2]=nodalConnBg[SPLIT_NODES_6_WO[2]]; conn[3]=nodalConnBg[SPLIT_NODES_6_WO[3]]; + conn[4]=nodalConnBg[SPLIT_NODES_6_WO[4]]; conn[5]=nodalConnBg[SPLIT_NODES_6_WO[5]]; conn[6]=nodalConnBg[SPLIT_NODES_6_WO[6]]; conn[7]=nodalConnBg[SPLIT_NODES_6_WO[7]]; + conn[8]=nodalConnBg[SPLIT_NODES_6_WO[8]]; conn[9]=nodalConnBg[SPLIT_NODES_6_WO[9]]; conn[10]=nodalConnBg[SPLIT_NODES_6_WO[10]]; conn[11]=nodalConnBg[SPLIT_NODES_6_WO[11]]; + conn[12]=nodalConnBg[SPLIT_NODES_6_WO[12]]; conn[13]=nodalConnBg[SPLIT_NODES_6_WO[13]]; conn[14]=nodalConnBg[SPLIT_NODES_6_WO[14]]; conn[15]=nodalConnBg[SPLIT_NODES_6_WO[15]]; + conn[16]=nodalConnBg[SPLIT_NODES_6_WO[16]]; conn[17]=nodalConnBg[SPLIT_NODES_6_WO[17]]; conn[18]=nodalConnBg[SPLIT_NODES_6_WO[18]]; conn[19]=nodalConnBg[SPLIT_NODES_6_WO[19]]; + conn[20]=nodalConnBg[SPLIT_NODES_6_WO[20]]; conn[21]=nodalConnBg[SPLIT_NODES_6_WO[21]]; conn[22]=nodalConnBg[SPLIT_NODES_6_WO[22]]; conn[23]=nodalConnBg[SPLIT_NODES_6_WO[23]]; + return ; + } + case GENERAL_24: + { + addCoords.resize(7*3); + tetrasNodalConn.resize(24*4); + int *conn(&tetrasNodalConn[0]); + double *tmp(&addCoords[18]); + tmp[0]=0.; tmp[1]=0.; tmp[2]=0.; + double *tmp2(&addCoords[0]); + for(int i=0;i<6;i++,tmp2+=3) + { + tmp2[0]=0.; tmp2[1]=0.; tmp2[2]=0.; + for(int j=0;j<4;j++,conn+=4) + { + tmp2[0]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+0]; + tmp2[1]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+1]; + tmp2[2]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+3]; + conn[0]=nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]; + conn[1]=nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+(j+1)%4]]; + conn[2]=-(i+1); conn[3]=-7; + } + tmp2[0]/=4.; tmp2[1]/=4.; tmp2[2]/=4.; + tmp[0]+=tmp2[0]; tmp[1]+=tmp2[1]; tmp[2]+=tmp2[2]; + } + tmp[0]/=6.; tmp[1]/=6.; tmp[2]/=6.; + return ; + } + case GENERAL_48: + { + addCoords.resize(19*3); + tetrasNodalConn.resize(48*4); + double *tmp2(&addCoords[0]),*tmp(&addCoords[0]); + for(int i=0;i<12;i++,tmp2+=3) + { + tmp2[0]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+0]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+0])/2.; + tmp2[1]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+1]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+1])/2.; + tmp2[2]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+2]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+2])/2.; + } + for(int i=0;i<7;i++,tmp2+=3) + { + tmp2[0]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+0]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+0])/2.; + tmp2[1]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+1]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+1])/2.; + tmp2[2]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+2]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+2])/2.; + } + int *conn(&tetrasNodalConn[0]); + std::vector dummy; + for(int i=0;i<8;i++) + { + std::vector c; + SplitHexa8IntoTetras(PLANAR_FACE_6,GENERAL_48_SUBZONES_2+i*8,GENERAL_48_SUBZONES_2+(i+1)*8,coords,c,dummy); + int *conn2(&c[0]); + for(int j=0;j<6;j++,conn+=4,conn2+=4) + { + conn[0]=conn2[0]>=0?nodalConnBg[conn2[0]]:conn2[0]; + conn[1]=conn2[1]>=0?nodalConnBg[conn2[1]]:conn2[1]; + conn[2]=conn2[2]>=0?nodalConnBg[conn2[2]]:conn2[2]; + conn[3]=conn2[3]>=0?nodalConnBg[conn2[3]]:conn2[3]; + } + } + return ; + } + default: + throw INTERP_KERNEL::Exception("SplitHexa8IntoTetras : invalid input policy ! Should be in [PLANAR_FACE_5,PLANAR_FACE_6,GENERAL_24,GENERAL_48] !"); + } + } + + void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords, + std::vector& tetrasNodalConn, std::vector& addCoords) throw(INTERP_KERNEL::Exception) + { + switch(gt) + { + case NORM_TETRA4: + { + std::size_t sz(std::distance(nodalConnBg,nodalConnEnd)); + if(sz!=4) + throw INTERP_KERNEL::Exception("SplitIntoTetras : input tetra do not have 4 nodes !"); + tetrasNodalConn.insert(tetrasNodalConn.end(),nodalConnBg,nodalConnEnd); + return ; + } + case NORM_HEXA8: + { + SplitHexa8IntoTetras(policy,nodalConnBg,nodalConnEnd,coords,tetrasNodalConn,addCoords); + return ; + } + case NORM_PYRA5: + { + std::size_t sz(std::distance(nodalConnBg,nodalConnEnd)); + if(sz!=5) + throw INTERP_KERNEL::Exception("SplitIntoTetras : input pyra5 do not have 5 nodes !"); + tetrasNodalConn.resize(8); + int *conn(&tetrasNodalConn[0]); + conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[4]; + conn[4]=nodalConnBg[0]; conn[5]=nodalConnBg[2]; conn[6]=nodalConnBg[3]; conn[7]=nodalConnBg[4]; + return ; + } + case NORM_PENTA6: + { + std::size_t sz(std::distance(nodalConnBg,nodalConnEnd)); + if(sz!=6) + throw INTERP_KERNEL::Exception("SplitIntoTetras : input penta6 do not have 6 nodes !"); + tetrasNodalConn.resize(12); + int *conn(&tetrasNodalConn[0]); + conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[3]; + conn[4]=nodalConnBg[3]; conn[5]=nodalConnBg[5]; conn[6]=nodalConnBg[4]; conn[7]=nodalConnBg[2]; + conn[8]=nodalConnBg[4]; conn[9]=nodalConnBg[2]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[3]; + return ; + } + case NORM_HEXGP12: + { + std::size_t sz(std::distance(nodalConnBg,nodalConnEnd)); + if(sz!=12) + throw INTERP_KERNEL::Exception("SplitIntoTetras : input octa12 (hexagone prism) do not have 12 nodes !"); + tetrasNodalConn.resize(48); + int *conn(&tetrasNodalConn[0]); + conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[5]; conn[3]=nodalConnBg[6]; + conn[4]=nodalConnBg[6]; conn[5]=nodalConnBg[11]; conn[6]=nodalConnBg[7]; conn[7]=nodalConnBg[5]; + conn[8]=nodalConnBg[7]; conn[9]=nodalConnBg[5]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[6]; + // + conn[12]=nodalConnBg[1]; conn[13]=nodalConnBg[4]; conn[14]=nodalConnBg[5]; conn[15]=nodalConnBg[7]; + conn[16]=nodalConnBg[7]; conn[17]=nodalConnBg[11]; conn[18]=nodalConnBg[10]; conn[19]=nodalConnBg[5]; + conn[20]=nodalConnBg[10]; conn[21]=nodalConnBg[5]; conn[22]=nodalConnBg[4]; conn[23]=nodalConnBg[7]; + // + conn[24]=nodalConnBg[1]; conn[25]=nodalConnBg[2]; conn[26]=nodalConnBg[4]; conn[27]=nodalConnBg[7]; + conn[28]=nodalConnBg[7]; conn[29]=nodalConnBg[10]; conn[30]=nodalConnBg[8]; conn[31]=nodalConnBg[4]; + conn[32]=nodalConnBg[8]; conn[33]=nodalConnBg[4]; conn[34]=nodalConnBg[2]; conn[35]=nodalConnBg[7]; + // + conn[36]=nodalConnBg[2]; conn[37]=nodalConnBg[3]; conn[38]=nodalConnBg[4]; conn[39]=nodalConnBg[8]; + conn[40]=nodalConnBg[8]; conn[41]=nodalConnBg[10]; conn[42]=nodalConnBg[9]; conn[43]=nodalConnBg[4]; + conn[44]=nodalConnBg[9]; conn[45]=nodalConnBg[4]; conn[46]=nodalConnBg[3]; conn[47]=nodalConnBg[8]; + return ; + } + case NORM_POLYHED: + { + std::size_t nbOfFaces(std::count(nodalConnBg,nodalConnEnd,-1)+1); + std::size_t nbOfTetra(std::distance(nodalConnBg,nodalConnEnd)-nbOfFaces+1); + addCoords.resize((nbOfFaces+1)*3); + tetrasNodalConn.resize(nbOfTetra*4); + int *conn(&tetrasNodalConn[0]); + const int *work(nodalConnBg); + double *tmp(&addCoords[0]),*tmp2(&addCoords[3*nbOfFaces]); + tmp2[0]=0.; tmp2[1]=0.; tmp2[2]=0.; + for(std::size_t i=0;i& tetrasNodalConn, std::vector& addCoords) throw(INTERP_KERNEL::Exception); + + void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords, + std::vector& tetrasNodalConn, std::vector& addCoords) throw(INTERP_KERNEL::Exception); + /** * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra. * @@ -280,7 +351,6 @@ namespace INTERP_KERNEL namespace INTERP_KERNEL { - /** * \brief Class calculating the volume of intersection between a tetrahedral target element and * source elements with triangular or quadratilateral faces. diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 9312b5fa6..38aab3c23 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -557,7 +557,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 @@ -1085,25 +1084,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); } } @@ -1238,33 +1222,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]; -- 2.39.2