From df2bff64c60f6f4a5af94711b1b9f3b3603173dc Mon Sep 17 00:00:00 2001 From: vbd Date: Tue, 18 Dec 2007 11:14:06 +0000 Subject: [PATCH] optimizing Intersector3D --- src/INTERP_KERNEL/Intersector3D.cxx | 52 ++++++++++++------- src/INTERP_KERNEL/Intersector3D.hxx | 3 +- src/INTERP_KERNEL/MeshUtils.hxx | 3 +- .../TransformedTriangle_math.cxx | 4 +- 4 files changed, 40 insertions(+), 22 deletions(-) diff --git a/src/INTERP_KERNEL/Intersector3D.cxx b/src/INTERP_KERNEL/Intersector3D.cxx index 9fd419442..605d09a84 100644 --- a/src/INTERP_KERNEL/Intersector3D.cxx +++ b/src/INTERP_KERNEL/Intersector3D.cxx @@ -28,11 +28,15 @@ namespace MEDMEM } /* - * Destructor + * Destruct * */ Intersector3D::~Intersector3D() { + for (map::iterator it=_cellmodel_map.begin(); + it!=_cellmodel_map.end(); + it++) + delete it->second; } /* @@ -124,23 +128,23 @@ namespace MEDMEM { // get cell model for the element - CELLMODEL cellModel(type); + CELLMODEL* cellModel=retrieveCellModel(type); - assert(cellModel.getDimension() == 3); + assert(cellModel->getDimension() == 3); assert(element >= 1); assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS)); // try to avoid reallocations with push_back() // 2 * the number of faces should be plenty - triangles.reserve(2 * cellModel.getNumberOfConstituents(1)); + triangles.reserve(2 * cellModel->getNumberOfConstituents(1)); // loop over faces - double** transformedNodes = new double*[cellModel.getNumberOfNodes()]; + double** transformedNodes = new double*[cellModel->getNumberOfNodes()]; bool isOutside[8] = {true, true, true, true, true, true, true, true}; bool isTargetOutside = false; // calculate the coordinates of the nodes - for(int i = 1; i <= cellModel.getNumberOfNodes() ; ++i) + for(int i = 1; i <= cellModel->getNumberOfNodes() ; ++i) { const double* node = getCoordsOfNode(i, element, _srcMesh); transformedNodes[i-1] = new double[3]; @@ -165,9 +169,9 @@ namespace MEDMEM if(!isTargetOutside) { - for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) + for(int i = 1 ; i <= cellModel->getNumberOfConstituents(1) ; ++i) { - const medGeometryElement faceType = cellModel.getConstituentType(1, i); + const medGeometryElement faceType = cellModel->getConstituentType(1, i); CELLMODEL faceModel(faceType); assert(faceModel.getDimension() == 2); @@ -177,10 +181,10 @@ namespace MEDMEM // get the nodes of the face for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) { - nodes[j-1] = cellModel.getNodeConstituent(1, i, j) - 1; + nodes[j-1] = cellModel->getNodeConstituent(1, i, j) - 1; assert(nodes[j-1] +1 >= 0); - assert(nodes[j-1] +1 <= cellModel.getNumberOfNodes()); + assert(nodes[j-1] +1 <= cellModel->getNumberOfNodes()); } // std::cout << "here3 : " << nodes[0] << "," << nodes[1] << ", " << nodes[2] << std::endl; // create transformed triangles from face @@ -226,7 +230,7 @@ namespace MEDMEM ++filtered; } - for(int i = 0 ; i < cellModel.getNumberOfNodes() ; ++i) + for(int i = 0 ; i < cellModel->getNumberOfNodes() ; ++i) { delete[] transformedNodes[i]; } @@ -239,16 +243,16 @@ namespace MEDMEM void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector& triangles, const TetraAffineTransform& T) const { // get cell model for the element - CELLMODEL cellModel(type); + CELLMODEL* cellModel=retrieveCellModel(type); assert(element >= 1); assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS)); - assert(cellModel.getDimension() == 3); + assert(cellModel->getDimension() == 3); // loop over faces - for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i) + for(int i = 1 ; i <= cellModel->getNumberOfConstituents(1) ; ++i) { - medGeometryElement faceType = cellModel.getConstituentType(1, i); + medGeometryElement faceType = cellModel->getConstituentType(1, i); CELLMODEL faceModel(faceType); assert(faceModel.getDimension() == 2); @@ -259,10 +263,10 @@ namespace MEDMEM for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j) { // offset of node from cellIdx - int localNodeNumber = cellModel.getNodeConstituent(1, i, j); + int localNodeNumber = cellModel->getNodeConstituent(1, i, j); assert(localNodeNumber >= 1); - assert(localNodeNumber <= cellModel.getNumberOfNodes()); + assert(localNodeNumber <= cellModel->getNumberOfNodes()); const double* node = getCoordsOfNode(localNodeNumber, element, _srcMesh); @@ -312,7 +316,19 @@ namespace MEDMEM } #endif - +CELLMODEL* retrieveCellModel(const medGeometryElement type) +{ + map::iterator it = _cellmodel_map.find(type); + if (it!=_cellmodel_map.end()) + return iter->second; + else + { + CELLMODEL* cellmodel=new CELLMODEL(type); + _cellmodel_map.insert(make_pair(type,cellmodel)); + return cellmodel; + } +} +} }; diff --git a/src/INTERP_KERNEL/Intersector3D.hxx b/src/INTERP_KERNEL/Intersector3D.hxx index dc8db25c6..3dafafaac 100644 --- a/src/INTERP_KERNEL/Intersector3D.hxx +++ b/src/INTERP_KERNEL/Intersector3D.hxx @@ -42,7 +42,8 @@ namespace MEDMEM const MESH& _srcMesh; const MESH& _targetMesh; - + map _cellmodel_map; + double _dp_accuracy_factor; }; // inline bool Intersector3D::isTriangleOutsideTetra(double* points) const diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx index cf9d2e9b0..6b84fdbdb 100644 --- a/src/INTERP_KERNEL/MeshUtils.hxx +++ b/src/INTERP_KERNEL/MeshUtils.hxx @@ -53,9 +53,10 @@ namespace INTERP_UTILS */ inline int getNumberOfNodesForType(const medGeometryElement type) { +#if DEBUG assert(type > 300); assert(type < 400); - +#endif // int(type) = yxx where y is dimension of entity (here 3) // and xx is the number of nodes of the element return static_cast(type) - 300; diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index ebbf768f3..c9928202e 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -49,7 +49,7 @@ namespace INTERP_UTILS const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; /// Threshold for resetting double and triple products to zero; ( F / f in Grandy ) - const long double TransformedTriangle::THRESHOLD_F = 20.0; + const long double TransformedTriangle::THRESHOLD_F = 500.0; /// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57] const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1; @@ -374,7 +374,7 @@ namespace INTERP_UTILS //? is this more stable? -> no subtraction // return asin( dotProd / ( lenNormal * lenEdgeVec ) ) + 3.141592625358979 / 2.0; - return 3.14159262535 - acos( dotProd / ( lenNormal * lenEdgeVec ) ); + return atan(1.0)*4.0 - acos( dotProd / ( lenNormal * lenEdgeVec ) ); } -- 2.39.2