From: ndjinga Date: Wed, 16 Apr 2008 07:18:08 +0000 (+0000) Subject: Changed abs into fabs X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=daa1420ff4a5e63688327820fa5542215c8317f7;p=tools%2Fmedcoupling.git Changed abs into fabs --- diff --git a/src/INTERP_KERNEL/IntersectorTetra.txx b/src/INTERP_KERNEL/IntersectorTetra.txx index c2ffd2db8..cd4090e6d 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.txx +++ b/src/INTERP_KERNEL/IntersectorTetra.txx @@ -252,7 +252,7 @@ namespace INTERP_KERNEL // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation // that should be used (which is equivalent to dividing by the determinant) - return std::abs(1.0 / _t->determinant() * totalVolume) ; + return std::fabs(1.0 / _t->determinant() * totalVolume) ; } } diff --git a/src/INTERP_KERNEL/PlanarIntersector.txx b/src/INTERP_KERNEL/PlanarIntersector.txx index d9e374e5b..98b94f5ce 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.txx +++ b/src/INTERP_KERNEL/PlanarIntersector.txx @@ -188,12 +188,12 @@ namespace INTERP_KERNEL else { std::cout << " Maille dégénérée " << "epsilon = " << epsilon << std::endl; - std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl; - std::cout << " distance2(&Coords_A[0],&Coords_A[i_A1])= " << distance2(&Coords_A[0],&Coords_A[i_A1]) << std::endl; - std::cout << "abs(normal_A) = " << abs(normal_A[0]) << " ; " <(&Coords_B[0],&Coords_B[i_B1])= " << distance2(&Coords_B[0],&Coords_B[i_B1]) << std::endl; - std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl; + std::cout << " i_A1= " << i_A1 << " i_A2= " << i_A2 << std::endl; + std::cout << " distance2(&Coords_A[0],&Coords_A[i_A1])= " << distance2(&Coords_A[0],&Coords_A[i_A1]) << std::endl; + std::cout << "abs(normal_A) = " << fabs(normal_A[0]) << " ; " <(&Coords_B[0],&Coords_B[i_B1])= " << distance2(&Coords_B[0],&Coords_B[i_B1]) << std::endl; + std::cout << "normal_B = " << normal_B[0] << " ; " << normal_B[1] << " ; " << normal_B[2] << std::endl; } } diff --git a/src/INTERP_KERNEL/PolygonAlgorithms.txx b/src/INTERP_KERNEL/PolygonAlgorithms.txx index 93f8fceee..df8dfde90 100644 --- a/src/INTERP_KERNEL/PolygonAlgorithms.txx +++ b/src/INTERP_KERNEL/PolygonAlgorithms.txx @@ -6,6 +6,8 @@ #include +using namespace std; + namespace INTERP_KERNEL { template @@ -38,7 +40,7 @@ namespace INTERP_KERNEL /******* Resolution of the linear system t1*AB+t2*DC=AC ***********/ det = determinant(AB,DC);//determinant of the first two coordinates - if(abs(det) >_Epsilon) + if(fabs(det) >_Epsilon) { inv_det = 1/det; t1 = determinant(AC,DC)*inv_det;//solves the linear system t1*AB+t2*DC=AC @@ -52,7 +54,7 @@ namespace INTERP_KERNEL return false; case 3://beware AB and CD may belong to a vertical plane det = determinant(&AB[1],&DC[1]);//determinant of the last two coefficients - if(abs(det) > _Epsilon) + if(fabs(det) > _Epsilon) { inv_det = 1/det; t1=(AC[1]*DC[2]-AC[2]*DC[1])*inv_det; @@ -61,7 +63,7 @@ namespace INTERP_KERNEL else //beware AB and CD may belong to a plane y = constant { det = AB[0]*DC[2]-AB[2]*DC[0]; - if(abs(det) > _Epsilon) + if(fabs(det) > _Epsilon) { inv_det = 1/det; t1=(AC[0]*DC[2]-AC[2]*DC[0])*inv_det; @@ -80,7 +82,7 @@ namespace INTERP_KERNEL return true; } } - else if(abs(t1) <= _Precision) + else if(fabs(t1) <= _Precision) { if( t2>_Precision && t2<1-_Precision)//vertex on an edge { @@ -105,9 +107,9 @@ namespace INTERP_KERNEL } } } - else if(abs(t2-1) <= _Precision)//vertex on a vertex (A=D), first run + else if(fabs(t2-1) <= _Precision)//vertex on a vertex (A=D), first run crossprod(A,C,E,_Vdouble);//angle between vectors AC and AE (E=vertex preceding A) - else if(abs(t2) <= _Precision)//vertex on a vertex (A=C), second run + else if(fabs(t2) <= _Precision)//vertex on a vertex (A=C), second run { double Vdoublebis[DIM]; crossprod(A,B,D,Vdoublebis); @@ -117,7 +119,7 @@ namespace INTERP_KERNEL for(int idim=0;idim(Vdoublebis,Vdoublebis) > _Epsilon) + else if(fabs(in_between)<=_Epsilon && dotprod(Vdoublebis,Vdoublebis) > _Epsilon) //ie _Vdouble=0, separation of overlaping edges at a double point { crossprod(A,E,B,_Vdouble); @@ -371,18 +373,20 @@ namespace INTERP_KERNEL _Terminus = N1 < 3 || N2<3; /* list of future events ordered according to their coordinates (x,y,z) (lexicographical order) */ - std::multimap< const double *, int, VertexLess > events; - typename std::multimap< const double *, int, VertexLess >::iterator mi1,mi2; + std::multimap< const double *, int, VertexLess > mmap_events; + typename std::list< pair< const double *, int > >::iterator mi1,mi2; std::multimap< int, std::pair< int,bool> >::iterator mi; /********** Initalisation of events with P1 and P2 vertices ************/ for(i_loc=0;i_loc > events(mmap_events.begin(),mmap_events.end()); + + if(!_Terminus) { /******** Treatment of the first vertex ********/ mi1=events.begin(); @@ -393,18 +397,21 @@ namespace INTERP_KERNEL _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob, false))); _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob, false))); mi1++; - + //std::cout<< "nb_prev= "<< 0 << " i_glob= " << i_glob << std::endl; + /******* Loop until the second polygon is reached *******/ while( !four_neighbours) { i_glob=(* mi1).second;//global index of vertex i nb_prev = _Status.count(i_glob);//counts the number of segments ending at i + + //std::cout<< "nb_prev= "<< nb_prev << " i_glob= " << i_glob << std::endl; switch (nb_prev) { case 1 : mi=_Status.find(i_glob);// pointer to the segment ending at i i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i - i_next= (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1; + i_next= (i_prev_glob - i_glob > 0) == (fabs(i_prev_glob - i_glob) == 1) ? i_glob - 1 : i_glob + 1; if(i_glob < N1) i_next_glob = (i_next +N1)%N1; else i_next_glob = (i_next-N1+N2)%N2 + N1; _Status.erase(mi); @@ -416,19 +423,18 @@ namespace INTERP_KERNEL case 0 : if( (i_glob < N1) != which_start) { - four_neighbours = true; - /* detection of double points */ - const double * Pi1=(* mi1).first; - mi2=mi1++; - const double * Pi2=(* mi2).first; - int i_double= (*mi2).second; - if(distance2(Pi1, Pi2) < _Epsilon && _Status.count(i_double)==1) + mi2=mi1; + mi2++; + /* detection of double points */ + if(distance2((* mi1).first, (*mi2).first) > _Epsilon) + four_neighbours = true; + else /* Rare pothological case: */ { - i_glob=i_double; - nb_prev=1; + //std::cout<<"coucou" << std::endl; + const pair< const double *, int > next_pt= *mi2; events.erase(mi2); - } - mi1--; + mi1=events.insert(mi1,next_pt); + } } break; default: @@ -438,12 +444,13 @@ namespace INTERP_KERNEL /******** Loop until a terminal point or crossing is reached ************/ while( !_Terminus) { + //std::cout<< "nb_prev= "<< nb_prev<< " nb_inter= " << _Inter.size()/DIM << std::endl; switch (nb_prev) { case 1 : mi=_Status.find(i_glob);// pointer to the segment ending at i i_prev_glob = ((*mi).second).first;//starting point of the segment ending at i - sign = (i_prev_glob - i_glob > 0) == (abs(i_prev_glob - i_glob) == 1) ? - 1 : + 1; + sign = (i_prev_glob - i_glob > 0) == (fabs(i_prev_glob - i_glob) == 1) ? - 1 : + 1; i_next_glob = i_glob+sign; _Is_in_intersection = ((*mi).second).second;//boolean that tells if i is in the intersection _Status.erase(mi); @@ -460,7 +467,7 @@ namespace INTERP_KERNEL if(i_glob < N1) for(idim=0;idim _Epsilon)//vertex clearly inside or outside - { + if(fabs(inside) > _Epsilon)//vertex clearly inside or outside + { + //std::cout<<"coucou1" << std::endl; if( _Is_in_intersection) { for(int idim=0;idim(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext); crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev); - inside_next= (dotprod(Vnext,V34)<0); - inside_prev= (dotprod(Vprev,V34)<0); + is_inside_next= (dotprod(Vnext,V34)<0); + is_inside_prev= (dotprod(Vprev,V34)<0); - if(!(inside_next || inside_prev)) return std::deque< double >(); + if(!(is_inside_next || is_inside_prev)) return std::deque< double >(); - if(inside_next) + if(is_inside_next) { _End_segments.push_back(std::make_pair(i_glob,i_next_glob)); add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob, @@ -515,7 +524,7 @@ namespace INTERP_KERNEL mi=_Status.find(j2_glob); ((* mi).second).second= !((* mi).second).second; } - if(inside_prev) + if(is_inside_prev) { _End_segments.push_back(std::make_pair(i_glob,i_prev_glob)); add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob, @@ -533,12 +542,12 @@ namespace INTERP_KERNEL { crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],Vnext); crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],Vprev); - inside_next= dotprod(Vnext,V12)<0; - inside_prev= dotprod(Vprev,V12)<0; + is_inside_next= dotprod(Vnext,V12)<0; + is_inside_prev= dotprod(Vprev,V12)<0; - if(!(inside_next || inside_prev)) return std::deque< double >(); + if(!(is_inside_next || is_inside_prev)) return std::deque< double >(); - if(inside_next) + if(is_inside_next) { _End_segments.push_back(std::make_pair(i_glob,i_next_glob)); add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_next], i_glob, i_next_glob, @@ -551,7 +560,7 @@ namespace INTERP_KERNEL mi=_Status.find(j4_glob); ((* mi).second).second= ! ((* mi).second).second; } - if(inside_prev) + if(is_inside_prev) { _End_segments.push_back(std::make_pair(i_glob,i_prev_glob)); add_crossing0(&Poly1[DIM*i_loc], &Poly1[DIM*i_prev], i_glob, i_prev_glob, @@ -567,32 +576,58 @@ namespace INTERP_KERNEL } else //vertices i, j1 and j3 share the same coordinates { - crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext); + crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_next],Vnext); crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j2], &Poly1[DIM*i_prev],Vprev); crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_next],V12); crossprod(&Poly1[DIM*i_loc], &Poly2[DIM*j4], &Poly1[DIM*i_prev],V34); - - inside_next= dotprod(Vnext,V12)<=_Epsilon; - inside_prev= dotprod(Vprev,V34)<=_Epsilon; - bool inside_j2= dotprod(Vnext,Vprev)<-_Epsilon; - - if(!( inside_next || inside_prev || inside_j2 )) return std::deque< double >(); - if(inside_next) _End_segments.push_back(std::make_pair(i_glob,i_next_glob)); - if(inside_prev) _End_segments.push_back(std::make_pair(i_glob,i_prev_glob)); - if(inside_j2) - { - _End_segments.push_back(std::make_pair(j1_glob,j2_glob)); - mi=_Status.find(j2_glob); - ((* mi).second).second= !((* mi).second).second; - } - if(!( (inside_next && (inside_prev || inside_j2)) || (inside_prev && inside_j2))) - { - _End_segments.push_back(std::make_pair(j3_glob,j4_glob)); - mi=_Status.find(j4_glob); - ((* mi).second).second= !((* mi).second).second; - } - _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob,inside_prev))); - _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob,inside_next))); + + double inside_next= dotprod(Vnext,V12); + double inside_prev= dotprod(Vprev,V34); + double inside_j2 = dotprod(Vnext,Vprev); + double inside_j4 = dotprod(V12,V34); + + std::map > which_is_inside; + which_is_inside[inside_next] = std:: make_pair(i_glob,i_next_glob); + which_is_inside[inside_prev] = std:: make_pair(i_glob,i_prev_glob); + which_is_inside[inside_j2] = std::make_pair(j1_glob,j2_glob); + which_is_inside[inside_j4] = std::make_pair(j3_glob,j4_glob); + + std::map >::iterator min = which_is_inside.begin(); + std::map >::iterator minext = min; + minext++; + std::map >::reverse_iterator max = which_is_inside.rbegin(); + std::multimap< int, std::pair< int,bool> >::iterator j2_in_status = _Status.find(((*min).second).second); + std::multimap< int, std::pair< int,bool> >::iterator j4_in_status = _Status.find(((*minext).second).second); + + if((*min).first < -_Epsilon) //there is someone clearly inside + { + _End_segments.push_back( (*min).second ); + _End_segments.push_back((* minext).second); + if(j2_in_status != _Status.end()) + ((*j2_in_status).second).second = ! ((*j2_in_status).second).second; + if(j4_in_status != _Status.end()) + ((*j4_in_status).second).second = ! ((*j4_in_status).second).second; + is_inside_next = ((*min).second).second == i_next_glob || ((*minext).second).second == i_next_glob; + is_inside_prev = ((*min).second).second == i_prev_glob || ((*minext).second).second == i_prev_glob; + } + else + if(fabs((*min).first) <= _Epsilon) //nobody is clearly inside but two segments are superposed + { + if(fabs((*max).first) > _Epsilon) + return std::deque< double >(); + else //all four segments are superposed + { + _End_segments.push_back(std::make_pair(i_glob,i_next_glob)); + _End_segments.push_back(std::make_pair(i_glob,i_prev_glob)); + is_inside_next= true; + is_inside_prev= true; + } + } + else //there is nobody inside + return std::deque< double >(); + + _Status.insert(std::make_pair(i_prev_glob,std::make_pair(i_glob,is_inside_prev))); + _Status.insert(std::make_pair(i_next_glob,std::make_pair(i_glob,is_inside_next))); } } } @@ -622,8 +657,10 @@ namespace INTERP_KERNEL if(NsubP>3) { std::map< int,int >::iterator mi_prev = subP.begin(); - std::map< int,int >::iterator mi = mi_prev++; - std::map< int,int >::iterator mi_next = mi++; + std::map< int,int >::iterator mi = mi_prev; + mi++; + std::map< int,int >::iterator mi_next = mi; + mi_next++; double directframe=0.; /* Check if the polygon subP is positively oriented */ @@ -631,7 +668,7 @@ namespace INTERP_KERNEL while(mi1 != subP.end() && distance2(&P[DIM*(*subP.begin()).second],&P[DIM*(*mi1).second])< epsilon) mi1++; std::map< int,int >::iterator mi2=mi1; - while(mi2 != subP.end() && abs(directframe)(&P[DIM* (*mi1).second], &P[DIM* (*subP.begin()).second], @@ -647,9 +684,9 @@ namespace INTERP_KERNEL &P[DIM* (*mi_prev).second], &P[DIM* (*mi_next).second], normal); if(directframe > -epsilon){ - mi = mi++; - mi_prev=mi_prev++; - mi_next=mi_next++; + mi ++; + mi_prev++; + mi_next++; } else { @@ -744,7 +781,7 @@ namespace INTERP_KERNEL int i1=1; while(i1(&P[0],&P[i1])< epsilon) i1++; int i2=i1+1; - while(i2(normal,normal))(normal,normal))(&P[i1], &P[0], &P[i2],normal); i2++; diff --git a/src/INTERP_KERNEL/TetraAffineTransform.cxx b/src/INTERP_KERNEL/TetraAffineTransform.cxx index ee31c2e2f..e018176a3 100644 --- a/src/INTERP_KERNEL/TetraAffineTransform.cxx +++ b/src/INTERP_KERNEL/TetraAffineTransform.cxx @@ -257,11 +257,11 @@ namespace INTERP_KERNEL // find pivot int i = k; - double max = std::abs(lu[3*idx[k] + k]); + double max = std::fabs(lu[3*idx[k] + k]); int row = i; while(i < 3) { - if(std::abs(lu[3*idx[i] + k]) > max) + if(std::fabs(lu[3*idx[i] + k]) > max) { max = fabs(lu[3*idx[i] + k]); row = i; diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index 49fbdcea8..f67963f2d 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -116,7 +116,7 @@ namespace INTERP_KERNEL const double term1 = _coords[5*pt1 + off1] * _coords[5*pt2 + off2]; const double term2 = _coords[5*pt1 + off2] * _coords[5*pt2 + off1]; - const long double delta = MULT_PREC_F * ( std::abs(term1) + std::abs(term2) ); + const long double delta = MULT_PREC_F * ( std::fabs(term1) + std::fabs(term2) ); if( epsilonEqual(_doubleProducts[8*seg + dp], 0.0, THRESHOLD_F * delta)) { @@ -125,7 +125,7 @@ namespace INTERP_KERNEL if(_doubleProducts[8*seg + dp] != 0.0) { LOG(5, "Double product for (seg,dp) = (" << seg << ", " << dp << ") = " ); - LOG(5, std::abs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" ); + LOG(5, std::fabs(_doubleProducts[8*seg + dp]) << " is imprecise, reset to 0.0" ); } #endif @@ -547,7 +547,7 @@ namespace INTERP_KERNEL #ifdef FIXED_DELTA const double delta = FIXED_DELTA; #else - const long double delta = MULT_PREC_F * (std::abs(p_term) + std::abs(q_term) + std::abs(r_term)); + const long double delta = MULT_PREC_F * (std::fabs(p_term) + std::fabs(q_term) + std::fabs(r_term)); #endif if( epsilonEqual( p_term + q_term + r_term, 0.0, THRESHOLD_F * delta) ) diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index b6ce5db3a..ad720feb2 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -110,7 +110,7 @@ namespace INTERP_KERNEL inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL) { return y < x ? x - y < errTol : y - x < errTol; - // return std::abs(x - y) < errTol; + // return std::fabs(x - y) < errTol; } /** @@ -128,12 +128,12 @@ namespace INTERP_KERNEL { // necessary for comparing values close to zero // in order to avoid division by very small numbers - if(std::abs(x - y) < absTol) + if(std::fabs(x - y) < absTol) { return true; } - const double relError = std::abs((x - y) / std::max(std::abs(x), std::abs(y))); + const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y))); return relError < relTol; }