Skip to content
Snippets Groups Projects
Commit ced5d60d authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

added Shewchuk's robust geometrical predicates (since they can be duplicated in
Tetgen/triangle I've put them in the gmsh namespace)
parent f1543de5
No related branches found
No related tags found
No related merge requests found
// $Id: 2D_DivAndConq.cpp,v 1.22 2006-09-11 15:23:54 remacle Exp $ // $Id: 2D_DivAndConq.cpp,v 1.23 2006-09-11 17:58:19 geuzaine Exp $
// //
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
// //
...@@ -183,53 +183,22 @@ Segment UpperCommonTangent(DT vl, DT vr) ...@@ -183,53 +183,22 @@ Segment UpperCommonTangent(DT vl, DT vr)
/* return 1 if the point k is NOT in the circumcircle of /* return 1 if the point k is NOT in the circumcircle of
triangle hij */ triangle hij */
#ifndef HAVE_TETGEN
int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k) int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
{ {
double xc, yc, rcarre, distca;
if((h == i) && (h == j) && (h == k)) {
Msg(GERROR, "3 identical points in Qtest");
return (0); /* returning 1 will cause looping for ever */
}
if(CircumCircle((double)pPointArray[h].where.h,
(double)pPointArray[h].where.v,
(double)pPointArray[i].where.h,
(double)pPointArray[i].where.v,
(double)pPointArray[j].where.h,
(double)pPointArray[j].where.v, &xc, &yc)) {
rcarre = square(xc - (double)pPointArray[i].where.h) +
square(yc - (double)pPointArray[i].where.v);
distca = square(xc - (double)pPointArray[k].where.h) +
square(yc - (double)pPointArray[k].where.v);
return (distca > rcarre);
}
else
return (0); /* point not in circle, because no circle ! */
}
#else
#include "tetgen.h"
int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
{
exactinit() ;
if((h == i) && (h == j) && (h == k)) { if((h == i) && (h == j) && (h == k)) {
Msg(GERROR, "3 identical points in Qtest"); Msg(GERROR, "3 identical points in Qtest");
return (0); /* returning 1 will cause looping for ever */ return (0); /* returning 1 will cause looping for ever */
} }
REAL pa[2] = {(REAL)pPointArray[h].where.h,(REAL)pPointArray[h].where.v}; double pa[2] = {(double)pPointArray[h].where.h, (double)pPointArray[h].where.v};
REAL pb[2] = {(REAL)pPointArray[i].where.h,(REAL)pPointArray[i].where.v}; double pb[2] = {(double)pPointArray[i].where.h, (double)pPointArray[i].where.v};
REAL pc[2] = {(REAL)pPointArray[j].where.h,(REAL)pPointArray[j].where.v}; double pc[2] = {(double)pPointArray[j].where.h, (double)pPointArray[j].where.v};
REAL pd[2] = {(REAL)pPointArray[k].where.h,(REAL)pPointArray[k].where.v}; double pd[2] = {(double)pPointArray[k].where.h, (double)pPointArray[k].where.v};
REAL result = incircle ( pa, pb, pc, pd ) * orient2d (pa,pb,pc); double result = gmsh::incircle(pa, pb, pc, pd) * gmsh::orient2d(pa, pb, pc);
// Msg(INFO, "Qtest %12.5E",result);
return (result < 0)?1:0; return (result < 0) ? 1 : 0;
} }
#endif
int merge(DT vl, DT vr) int merge(DT vl, DT vr)
......
# $Id: Makefile,v 1.28 2006-08-04 14:28:03 geuzaine Exp $ # $Id: Makefile,v 1.29 2006-09-11 17:58:19 geuzaine Exp $
# #
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
# #
...@@ -27,6 +27,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE} ...@@ -27,6 +27,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Numeric.cpp\ SRC = Numeric.cpp\
EigSolve.cpp\ EigSolve.cpp\
predicates.cpp\
gsl_newt.cpp\ gsl_newt.cpp\
gsl_brent.cpp gsl_brent.cpp
......
// $Id: Numeric.cpp,v 1.27 2006-08-16 05:25:22 geuzaine Exp $ // $Id: Numeric.cpp,v 1.28 2006-09-11 17:58:19 geuzaine Exp $
// //
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
// //
...@@ -58,6 +58,11 @@ int check_gsl() ...@@ -58,6 +58,11 @@ int check_gsl()
} }
// set new error handler // set new error handler
gsl_set_error_handler(&new_handler); gsl_set_error_handler(&new_handler);
// initilize robust geometric predicates
gmsh::exactinit() ;
return 1; return 1;
} }
...@@ -65,6 +70,10 @@ int check_gsl() ...@@ -65,6 +70,10 @@ int check_gsl()
int check_gsl() int check_gsl()
{ {
// initilize robust geometric predicates
gmsh::exactinit() ;
return 1; return 1;
} }
......
...@@ -97,4 +97,13 @@ void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, ...@@ -97,4 +97,13 @@ void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
void newt(double x[], int n, int *check, void newt(double x[], int n, int *check,
void (*vecfunc)(int, double [], double [])); void (*vecfunc)(int, double [], double []));
/* Robust geometrical predicates */
namespace gmsh{
double exactinit();
double incircle(double *pa, double *pb, double *pc, double *pd);
double insphere(double *pa, double *pb, double *pc, double *pd, double *pe);
double orient2d(double *pa, double *pb, double *pc);
double orient3d(double *pa, double *pb, double *pc, double *pd);
}
#endif #endif
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment