From: vbd Date: Fri, 3 Aug 2007 12:52:06 +0000 (+0000) Subject: staffan : X-Git-Tag: trio_trio_coupling~81 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9cba1d446831840539ad0ee2519a2dcf11b4a76c;p=tools%2Fmedcoupling.git staffan : * correction of the MeshElement::getCoordsOfNode() which was totally wrong - it worked well with only one element because of the correspondance between the arrays connectivityIndex and connectivity that exists in this case --- diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx index 8e538178d..9355fbd27 100644 --- a/src/INTERP_KERNEL/Interpolation3D.cxx +++ b/src/INTERP_KERNEL/Interpolation3D.cxx @@ -280,11 +280,16 @@ namespace MEDMEM double Interpolation3D::calculateIntersectionVolume(const MeshElement& srcElement, const MeshElement& targetElement) { - // std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl; + //std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl; // (a), (b) and (c) not yet implemented // (d) : without fine-level filtering (a) - (c) for the time being - + + std::cout << "Source : "; + srcElement.dumpCoords(); + std::cout << "Target : "; + targetElement.dumpCoords(); + // get array of points of target tetraeder const double* tetraCorners[4]; for(int i = 0 ; i < 4 ; ++i) @@ -294,10 +299,8 @@ namespace MEDMEM // create AffineTransform TetraAffineTransform T( tetraCorners ); - std::cout << "Transform : " << std::endl; - T.dump(); - std::cout << std::endl; - + + // check if we have planar tetra element if(T.determinant() == 0.0) { // tetra is planar @@ -305,6 +308,10 @@ namespace MEDMEM return 0.0; } + std::cout << "Transform : " << std::endl; + T.dump(); + std::cout << std::endl; + // triangulate source element faces (assumed tetraeder for the time being) // do nothing vector triangles; @@ -328,7 +335,7 @@ namespace MEDMEM //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that // should be used - return 1.0 / T.determinant() * volume ; + return std::abs(1.0 / T.determinant() * volume) ; } } diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx index f57e20f27..1501b5322 100644 --- a/src/INTERP_KERNEL/MeshElement.cxx +++ b/src/INTERP_KERNEL/MeshElement.cxx @@ -99,9 +99,12 @@ namespace INTERP_UTILS */ const double* MeshElement::getCoordsOfNode(int node) const { + assert(node >= 1); + assert(node <= getNumberNodes()); const int nodeOffset = node - 1; const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1; - return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + 3*nodeOffset]); + const int connIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1; + return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connIdx]); } /** @@ -130,9 +133,10 @@ namespace INTERP_UTILS CELLMODEL faceModel(faceType); assert(faceModel.getDimension() == 2); + assert(faceModel.getNumberOfNodes() == 3); double transformedNodes[3 * faceModel.getNumberOfNodes()]; - const double* coords = _mesh->getCoordinates(MED_FULL_INTERLACE); + // loop over nodes of face for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) @@ -149,11 +153,13 @@ namespace INTERP_UTILS //{ not totally efficient since we transform each node once per face T.apply(&transformedNodes[3*(j-1)], node); + std::cout << "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to " << vToStr(&transformedNodes[3*(j-1)]) << std::endl; + } assert(faceType == MED_TRIA3); - // create triangles from face + // create transformed triangles from face switch(faceType) { case MED_TRIA3: @@ -218,11 +224,22 @@ namespace INTERP_UTILS } } #endif + int MeshElement::getIndex() const { return _index + 1; } + void MeshElement::dumpCoords() const + { + std::cout << "Element " << _index + 1 << " has nodes " << std::endl; + for(int i = 1 ; i <= getNumberNodes() ; ++i) + { + std::cout << vToStr(getCoordsOfNode(i)) << ", "; + } + std::cout << std::endl; + } + /// ElementBBoxOrder bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2) diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx index ab339da23..f58c36d0f 100644 --- a/src/INTERP_KERNEL/MeshElement.hxx +++ b/src/INTERP_KERNEL/MeshElement.hxx @@ -97,6 +97,8 @@ namespace INTERP_UTILS int getIndex() const; + void dumpCoords() const; + private: const int _index; const MEDMEM::MESH* _mesh; diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx index c887c75d0..3da899fd9 100644 --- a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx +++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx @@ -105,8 +105,6 @@ void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const ch m = interpolator->interpol_maillages(sMesh, tMesh); - dumpIntersectionMatrix(m); - std::cout << "Intersection calculation done. " << std::endl << std::endl; } diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx index c564438d0..159b973ac 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.hxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx @@ -17,8 +17,8 @@ namespace INTERP_UTILS TetraAffineTransform(const double** pts) { - std::cout << "Creating transform from tetraeder : " << std::endl; - std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl; + // std::cout << "Creating transform from tetraeder : " << std::endl; + // std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl; #if 0 do { diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx index f91241aaf..b601b7fe2 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx @@ -109,12 +109,12 @@ namespace INTERP_UTILS // calculate point for(int i = 0; i < 3; ++i) { - std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl; + // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl; pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] + alpha * COORDS_TET_CORNER[3*corners[1] + i]; - std::cout << pt[i] << std::endl; - //assert(pt[i] >= 0.0); - //assert(pt[i] <= 1.0); + // std::cout << pt[i] << std::endl; + assert(pt[i] >= 0.0); + assert(pt[i] <= 1.0); } } @@ -429,8 +429,8 @@ namespace INTERP_UTILS // special case : facet H = 0 bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet); - std::cout << "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) << std::endl; - std::cout << "c2 = " << cVals[2] << ", c3 = " << cVals[3] << std::endl; + // std::cout << "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) << std::endl; + // std::cout << "c2 = " << cVals[2] << ", c3 = " << cVals[3] << std::endl; return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0); } @@ -657,7 +657,7 @@ namespace INTERP_UTILS const double cQR = calcStableC(QR, DoubleProduct(edge)); const double cRP = calcStableC(RP, DoubleProduct(edge)); - std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; + // std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; // if two or more c-values are zero we disallow x-edge intersection // Grandy, p.446 @@ -665,7 +665,7 @@ namespace INTERP_UTILS if(numZeros >= 2 ) { - std::cout << "TriangleSurroundsEdge test fails due to too many 0 dp" << std::endl; + // std::cout << "TriangleSurroundsEdge test fails due to too many 0 dp" << std::endl; } return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2; @@ -700,7 +700,7 @@ namespace INTERP_UTILS const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]); //? should equality with zero use epsilon? - std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl; + // std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl; return (t1*t2 <= 0.0) && (t1 - t2 != 0.0); } @@ -807,7 +807,7 @@ namespace INTERP_UTILS //? NB here we have no correction for precision - is this good? // Our authority Grandy says nothing - std::cout << "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; + // std::cout << "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl; return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 ); }