From 9be637ee6ebc843f3d44052b2e650176e3c686bc Mon Sep 17 00:00:00 2001 From: eap Date: Fri, 9 Jan 2009 16:43:50 +0000 Subject: [PATCH] MEDMEM Industrialization 2008 new version --- src/INTERP_KERNEL/Interpolation3D.txx | 135 ++++++++++++++++++---- src/INTERP_KERNEL/InterpolationPlanar.txx | 77 ++++++++++-- 2 files changed, 181 insertions(+), 31 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx index b0fb83433..214d1a609 100644 --- a/src/INTERP_KERNEL/Interpolation3D.txx +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -80,7 +80,10 @@ namespace INTERP_KERNEL * */ template - void Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result) + void Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, + const MyMeshType& targetMesh, + MatrixType& result, + const std::string& method) { typedef typename MyMeshType::MyConnType ConnType; const NumberingPolicy numPol=MyMeshType::My_numPol; @@ -107,9 +110,22 @@ namespace INTERP_KERNEL targetElems[i] = new MeshElement(OTT::indFC(i), targetMesh); } - // create empty maps for all source elements - result.resize(numTargetElems); + bool srcCornerCoeff = ( method[1] == '1' ); // "P1P0" or "P1dP0" + bool tgtCornerCoeff = ( method[3] == '1' ); // "P0P1" or "P0P1d" + bool nodeCoeff = ( method.size() < 5 ); // no P1d + bool cornerCoeff = ( srcCornerCoeff || tgtCornerCoeff ); + const int nbCorners = 4; // coeffs for corners are calculated for tetras only + double baryCentreBuf[3], baryCoordBuf[nbCorners]; + double* baryCentre = cornerCoeff ? baryCentreBuf : 0; + double* baryCoords = cornerCoeff ? baryCoordBuf : 0; + // create empty maps for all source elements + if ( tgtCornerCoeff && nodeCoeff) + result.resize(targetMesh.getNumberOfNodes()); + else if ( tgtCornerCoeff && !nodeCoeff) + result.resize(numTargetElems * nbCorners); + else + result.resize(numTargetElems); #ifdef USE_RECURSIVE_BBOX_FILTER @@ -196,20 +212,56 @@ namespace INTERP_KERNEL iter != currNode->getSrcRegion().getEndElements() ; ++iter) { const ConnType srcIdx = (*iter)->getIndex(); - double* baryCenre = - result.getBaryCentre(OTT::ind2C(targetIdx),srcIdx); - const double vol = intersector->intersectSourceCell(srcIdx,baryCenre); + const double vol = intersector->intersectSourceCell(srcIdx,baryCentre); if(vol != 0.0) { - result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx,vol)); - LOG3(2, "Result : V (" << srcIdx << "- " << targetIdx << ") = " << result[ OTT::ind2C(srcIdx) ][targetIdx] ); + if ( !baryCentre ) + { + result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx,vol)); + LOG3(2, "Result : V (" << srcIdx << "- " << targetIdx << ") = " << result[ OTT::ind2C(srcIdx) ][targetIdx] ); + } + else if ( tgtCornerCoeff ) + { + // this code intended for tetras only -> 4 nodes + getBarycentricCoordinates( baryCentre, targetIdx, + targetMesh, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per target node + for ( int n = 0; n < 4; ++n ) { + ConnType node = getGlobalNumberOfNode(n, targetIdx, targetMesh); + result[node].insert(make_pair(srcIdx,vol*baryCoords[n])); + } + } + else { // several coeffs per target node + int tgtIdx0 = nbCorners * OTT::ind2C( targetIdx ); + result[tgtIdx0+0].insert(make_pair(srcIdx,vol*baryCoords[0])); + result[tgtIdx0+1].insert(make_pair(srcIdx,vol*baryCoords[1])); + result[tgtIdx0+2].insert(make_pair(srcIdx,vol*baryCoords[2])); + result[tgtIdx0+3].insert(make_pair(srcIdx,vol*baryCoords[3])); + } + } + else + { + getBarycentricCoordinates( baryCentre, srcIdx, + srcMesh, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per source node + for ( int n = 0; n < 4; ++n ) { + ConnType node = getGlobalNumberOfNode(n, srcIdx, srcMesh); + result[OTT::ind2C(targetIdx)].insert(make_pair(OTT::indFC(node),vol*baryCoords[n])); + } + } + else { // several coeffs per source node + int srcIdx0 = OTT::indFC( nbCorners * OTT::ind2C(srcIdx) ); + for ( int n = 0; n < 4; ++n ) + result[OTT::ind2C(targetIdx)].insert(std::make_pair(srcIdx0+n,vol*baryCoords[0])); + } + } } } delete intersector; - } + } else // recursion { @@ -334,32 +386,67 @@ namespace INTERP_KERNEL IntersectorTetra intersector(srcMesh, targetMesh, targetIdx); for(vector::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter) - { + { - const int srcIdx = *iter + 1; - double baryCenre* = result.getBaryCentre(targetIdx-1,srcIdx); - const double vol = intersector.intersectSourceCell(srcIdx, baryCenre); + const int srcIdx = *iter + 1; + const double vol = intersector.intersectSourceCell(srcIdx, baryCenre); - if(vol != 0.0) - { - result[targetIdx - 1].insert(make_pair(srcIdx, vol)); + if(vol != 0.0) + { + if ( !baryCentre ) + { + result[targetIdx - 1].insert(make_pair(srcIdx, vol)); + } + else if ( tgtCornerCoeff ) + { + // this code intended for tetras only -> 4 nodes + getBarycentricCoordinates( baryCentre, targetIdx, targetMesh, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per target node + for ( int n = 0; n < 4; ++n ) { + ConnType node = getGlobalNumberOfNode(n, targetIdx, myMesh_P); + result[node].insert(make_pair(srcIdx,vol*baryCoords[n])); + } } - + else { // several coeffs per target node + int tgtIdx0 = nbCorners * OTT::ind2C( targetIdx ); + result[tgtIdx0+0].insert(make_pair(srcIdx,vol*baryCoords[0])); + result[tgtIdx0+1].insert(make_pair(srcIdx,vol*baryCoords[1])); + result[tgtIdx0+2].insert(make_pair(srcIdx,vol*baryCoords[2])); + result[tgtIdx0+3].insert(make_pair(srcIdx,vol*baryCoords[3])); + } + } + else + { + getBarycentricCoordinates( baryCentre, srcIdx, srcMesh, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per source node + for ( n = 0; n < 4; ++n ) { + ConnType node = getGlobalNumberOfNode(n, srcIdx, myMesh_S); + result[OTT::ind2C(targetIdx)].insert(make_pair(OTT::indFC(node),vol*baryCoords[n])); + } + } + else { // several coeffs per source node + int srcIdx0 = OTT::indFC( nbCorners * OTT::ind2C(srcIdx) ); + result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx0+0,vol*baryCoords[0])); + result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx0+1,vol*baryCoords[1])); + result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx0+2,vol*baryCoords[2])); + result[OTT::ind2C(targetIdx)].insert(make_pair(srcIdx0+3,vol*baryCoords[3])); + } + } } + } } - #endif // free allocated memory for(unsigned long i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } + { + delete srcElems[i]; + } for(unsigned long i = 0 ; i < numTargetElems ; ++i) - { - delete targetElems[i]; - } + { + delete targetElems[i]; + } } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index e0559be37..f9061110f 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -32,6 +32,7 @@ #include "VectorUtils.hxx" #include "BBTree.txx" #include +#include "MeshUtils.hxx" namespace INTERP_KERNEL { @@ -104,7 +105,10 @@ namespace INTERP_KERNEL */ template template - void InterpolationPlanar::interpolateMeshes(const MyMeshType& myMesh_S, const MyMeshType& myMesh_P, MatrixType& result) + void InterpolationPlanar::interpolateMeshes(const MyMeshType& myMesh_S, + const MyMeshType& myMesh_P, + MatrixType& result, + const std::string& method) { static const int SPACEDIM=MyMeshType::MY_SPACEDIM; typedef typename MyMeshType::MyConnType ConnType; @@ -184,7 +188,21 @@ namespace INTERP_KERNEL long end_filtering=clock(); - result.resize(nbMaille_P);//on initialise. + bool srcCornerCoeff = ( method[1] == '1' ); // "P1P0" or "P1dP0" + bool tgtCornerCoeff = ( method[3] == '1' ); // "P0P1" or "P0P1d" + bool nodeCoeff = ( method.size() < 5 ); // no P1d + bool cornerCoeff = ( srcCornerCoeff || tgtCornerCoeff ); + const int nbCorners = 3; // coeffs for corners are calculated for triangles only + double baryCentreBuf[SPACEDIM], baryCoordBuf[nbCorners]; + double* baryCentre = cornerCoeff ? baryCentreBuf : 0; + double* baryCoords = cornerCoeff ? baryCoordBuf : 0; + + if ( tgtCornerCoeff && nodeCoeff) + result.resize(myMesh_P.getNumberOfNodes()); + else if ( tgtCornerCoeff && !nodeCoeff) + result.resize(nbMaille_P * nbCorners); + else + result.resize(nbMaille_P);//on initialise. /****************************************************/ /* Loop on the target cells - core of the algorithm */ @@ -210,10 +228,9 @@ namespace INTERP_KERNEL //BBTree structure returns numbers between 0 and n-1 int i_S=intersecting_elems[ielem]; //MN: Global number of cell ? int nb_nodesS=connIndxS[i_S+1]-connIndxS[i_S]; - double* baryCenre = result.getBaryCentre( i_P, OTT::indFC(i_S) ); double surf=intersector->intersectCells(OTT::indFC(i_P), OTT::indFC(i_S), - nb_nodesP,nb_nodesS,baryCenre); + nb_nodesP,nb_nodesS,baryCentre); //filtering out zero surfaces and badly oriented surfaces // orientation = -1,0,1 @@ -221,20 +238,66 @@ namespace INTERP_KERNEL // 0 : the intersection is always taken into account // 1 : the intersection is taken into account if target and cells have the same orientation int orientation=InterpolationOptions::getOrientation(); - if (( surf > 0.0 && orientation >=0 ) || ( surf < 0.0 && orientation <=0 )) - result[i_P].insert(std::make_pair(OTT::indFC(i_S),surf)); + if (( surf > 0.0 && orientation >=0 ) || ( surf < 0.0 && orientation <=0 )) { + if ( !baryCentre ) + { + result[i_P].insert(std::make_pair(OTT::indFC(i_S),surf)); + } + else if ( tgtCornerCoeff ) + { + // this code intended for triangles only -> 3 nodes + getBarycentricCoordinates( baryCentre, OTT::indFC(i_P), + myMesh_P, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per target node + cout << "P elem " << OTT::indFC(i_P) + << ", S elem " << OTT::indFC(i_S) + << ", surf: " << surf + << endl; + for ( int n = 0; n < 3; ++n ) { + ConnType node = getGlobalNumberOfNode(n, OTT::indFC(i_P), myMesh_P); + result[node].insert(make_pair(OTT::indFC(i_S),surf*baryCoords[n])); + cout << "node " << node << " baryCoord[" << n << "]=" << baryCoords[n] + << ", surf*bc = " << surf*baryCoords[n] << endl; + } + } + else { // several coeffs per target node + int tgtIdx0 = nbCorners * i_P; + result[tgtIdx0+0].insert(make_pair(OTT::indFC(i_S),surf*baryCoords[0])); + result[tgtIdx0+1].insert(make_pair(OTT::indFC(i_S),surf*baryCoords[1])); + result[tgtIdx0+2].insert(make_pair(OTT::indFC(i_S),surf*baryCoords[2])); + } + } + else + { + getBarycentricCoordinates( baryCentre, OTT::indFC(i_S), + myMesh_S, baryCoords ); + if ( nodeCoeff ) { // 1 coeff per source node + for ( int n = 0; n < 3; ++n ) { + ConnType node = getGlobalNumberOfNode(n, OTT::indFC(i_S), myMesh_S); + result[i_P].insert(make_pair(OTT::indFC(node),surf*baryCoords[n])); + } + } + else { // several coeffs per source node + int srcIdx0 = OTT::indFC(nbCorners * i_S); + result[i_P].insert(make_pair(srcIdx0+0,surf*baryCoords[0])); + result[i_P].insert(make_pair(srcIdx0+1,surf*baryCoords[1])); + result[i_P].insert(make_pair(srcIdx0+2,surf*baryCoords[2])); + } + } + } counter++; } intersecting_elems.clear(); } delete intersector; + /***********************************/ /* DEBUG prints */ /***********************************/ if (InterpolationOptions::getPrintLevel() >=1) - { + { long end_intersection=clock(); // if (_printLevel >=2) // { -- 2.39.2