From a351e519f2afd9a39a8ae5617cdd2d8ddceba9f6 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Thu, 31 Jul 2014 15:30:53 +0000
Subject: [PATCH] levelset yarn
---
Geo/gmshLevelset.cpp | 210 +++++++++++---------
Geo/gmshLevelset.h | 458 ++++++++++++++++++++++---------------------
Numeric/Numeric.cpp | 224 ++++++++++++++++-----
Numeric/Numeric.h | 13 +-
4 files changed, 543 insertions(+), 362 deletions(-)
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 96599d5b7f..83842d2380 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -27,13 +27,13 @@ void insertActiveCells(double x, double y, double z, double rmax,
{
int id1 = box.getCellContainingPoint(x - rmax, y - rmax, z - rmax);
int id2 = box.getCellContainingPoint(x + rmax, y + rmax, z + rmax);
- int i1, j1 ,k1;
+ int i1, j1, k1;
box.getCellIJK(id1, i1, j1, k1);
int i2, j2, k2;
box.getCellIJK(id2, i2, j2, k2);
- for (int i = i1; i <= i2; i++)
- for (int j = j1; j <= j2; j++)
- for (int k = k1; k <= k2; k++)
+ for(int i = i1; i <= i2; i++)
+ for(int j = j1; j <= j2; j++)
+ for(int k = k1; k <= k2; k++)
box.insertActiveCell(box.getCellIndex(i, j, k));
}
@@ -116,15 +116,15 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
{
std::vector<SPoint3> nodes;
std::vector<int> indices;
- for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
- it != box.nodalValuesEnd(); ++it){
+ for(cartesianBox<double>::valIter it = box.nodalValuesBegin();
+ it != box.nodalValuesEnd(); ++it){
nodes.push_back(box.getNodeCoordinates(it->first));
indices.push_back(it->first);
}
//Msg::Info(" %d nodes in the grid at level %d", (int)nodes.size(), box.getLevel());
std::vector<double> dist, localdist;
std::vector<SPoint3> dummy;
- for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
+ for(GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
if((*fit)->geomType() == GEntity::DiscreteSurface){
for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
std::vector<double> iDistances;
@@ -144,7 +144,7 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
if(dist.empty())
dist = localdist;
else{
- for (unsigned int j = 0; j < localdist.size(); j++){
+ for(unsigned int j = 0; j < localdist.size(); j++){
dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
}
}
@@ -158,7 +158,7 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
}
}
- for (unsigned int j = 0; j < dist.size(); j++)
+ for(unsigned int j = 0; j < dist.size(); j++)
box.setNodalValue(indices[j], dist[j]);
if(box.getChildBox()) computeLevelset(gm, *box.getChildBox());
@@ -203,16 +203,16 @@ inline double evalRadialFnDer(int p, int index, double dx, double dy, double dz,
double ep)
{
double r2 = dx*dx+dy*dy+dz*dz; //r^2
- switch (index) {
+ switch(index) {
case 0 : // GA
- switch (p) {
+ switch(p) {
case 0: return exp(-ep*ep*r2);
case 1: return -2*ep*ep*dx*exp(-ep*ep*r2); // _x
case 2: return -2*ep*ep*dy*exp(-ep*ep*r2); // _y
case 3: return -2*ep*ep*dz*exp(-ep*ep*r2); // _z
}
case 1 : //MQ
- switch (p) {
+ switch(p) {
case 0: return sqrt(1+ep*ep*r2);
case 1: return ep*ep*dx/sqrt(1+ep*ep*r2);
case 2: return ep*ep*dy/sqrt(1+ep*ep*r2);
@@ -396,16 +396,16 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv)
// level set defined by points (RBF interpolation)
fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
- const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2) const
+ const fullMatrix<double> &nodes1,
+ const fullMatrix<double> &nodes2) const
{
int m = nodes2.size1();
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
double eps =0.5/delta;
- for (int i = 0; i < m; i++) {
- for (int j = 0; j < n; j++) {
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < n; j++) {
double dx = nodes2(i,0)-nodes1(j,0);
double dy = nodes2(i,1)-nodes1(j,1);
double dz = nodes2(i,2)-nodes1(j,2);
@@ -417,15 +417,15 @@ fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
}
void gLevelsetPoints::RbfOp(int p, int index,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &D,
- bool isLocal) const
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &nodes,
+ fullMatrix<double> &D,
+ bool isLocal) const
{
fullMatrix<double> rbfMatB = generateRbfMat(p,index, cntrs,nodes);
fullMatrix<double> rbfInvA;
- if (isLocal){
+ if(isLocal){
rbfInvA = generateRbfMat(0,index, cntrs,cntrs);
rbfInvA.invertInPlace();
}
@@ -450,8 +450,8 @@ void gLevelsetPoints::evalRbfDer(int p, int index,
}
void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
- fullMatrix<double> &level_set_nodes,
- fullMatrix<double> &level_set_funvals)
+ fullMatrix<double> &level_set_nodes,
+ fullMatrix<double> &level_set_funvals)
{
int numNodes = cntrs.size1();
int nTot = 3*numNodes;
@@ -464,17 +464,17 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
// Computes the normal vectors to the surface at each node
double dist_min = 1.e6;
double dist_max = 1.e-6;
- for (int i = 0; i < numNodes; ++i){
+ for(int i = 0; i < numNodes; ++i){
ONES(i,0)=1.0;
cntrsPlus(i,0) = cntrs(i,0);
cntrsPlus(i,1) = cntrs(i,1);
cntrsPlus(i,2) = cntrs(i,2);
for(int j = i+1; j < numNodes; j++){
double dist = sqrt((cntrs(i,0)-cntrs(j,0))*(cntrs(i,0)-cntrs(j,0))+
- (cntrs(i,1)-cntrs(j,1))*(cntrs(i,1)-cntrs(j,1))+
- (cntrs(i,2)-cntrs(j,2))*(cntrs(i,2)-cntrs(j,2)));
- if (dist<dist_min) dist_min = dist;
- if (dist>dist_max) dist_max = dist;
+ (cntrs(i,1)-cntrs(j,1))*(cntrs(i,1)-cntrs(j,1))+
+ (cntrs(i,2)-cntrs(j,2))*(cntrs(i,2)-cntrs(j,2)));
+ if(dist<dist_min) dist_min = dist;
+ if(dist>dist_max) dist_max = dist;
}
}
ONES(numNodes,0) = -1.0;
@@ -486,7 +486,7 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
evalRbfDer(1, 1, cntrsPlus,cntrs,ONES,sx, true);
evalRbfDer(2, 1, cntrsPlus,cntrs,ONES,sy, true);
evalRbfDer(3, 1, cntrsPlus,cntrs,ONES,sz, true);
- for (int i = 0; i < numNodes; ++i){
+ for(int i = 0; i < numNodes; ++i){
normFactor = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
sx(i,0) = sx(i,0)/normFactor;
sy(i,0) = sy(i,0)/normFactor;
@@ -494,8 +494,8 @@ void gLevelsetPoints::setup_level_set(const fullMatrix<double> &cntrs,
norms(i,0) = sx(i,0);norms(i,1) = sy(i,0);norms(i,2) = sz(i,0);
}
- for (int i = 0; i < numNodes; ++i){
- for (int j = 0; j < 3; ++j){
+ for(int i = 0; i < numNodes; ++i){
+ for(int j = 0; j < 3; ++j){
level_set_nodes(i,j) = cntrs(i,j);
level_set_nodes(i+numNodes,j) = cntrs(i,j)-delta*norms(i,j);
level_set_nodes(i+2*numNodes,j) = cntrs(i,j)+delta*norms(i,j);
@@ -530,7 +530,7 @@ gLevelsetPoints::gLevelsetPoints(const gLevelsetPoints &lv)
double gLevelsetPoints::operator()(double x, double y, double z) const
{
if(mapP.empty())
- Msg::Info("Levelset Points: call computeLS() before calling operator()");
+ Msg::Info("Levelset Points : call computeLS() before calling operator()\n");
SPoint3 sp(x,y,z);
std::map<SPoint3,double>::const_iterator it = mapP.find(sp);
@@ -549,13 +549,13 @@ double gLevelsetPoints::operator()(double x, double y, double z) const
void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert)
{
fullMatrix<double> xyz_eval(vert.size(), 3), surf_eval(vert.size(), 1);
- for (unsigned int i = 0; i < vert.size(); i++){
+ for(unsigned int i = 0; i < vert.size(); i++){
xyz_eval(i,0) = vert[i]->x();
xyz_eval(i,1) = vert[i]->y();
xyz_eval(i,2) = vert[i]->z();
}
evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval);
- for (unsigned int i = 0; i < vert.size(); i++){
+ for(unsigned int i = 0; i < vert.size(); i++){
mapP[SPoint3(vert[i]->x(), vert[i]->y(),vert[i]->z())] = surf_eval(i,0);
}
}
@@ -718,7 +718,7 @@ gLevelsetShamrock::gLevelsetShamrock(double _xmid, double _ymid, double _zmid,
// creating the iso-zero
double angle = 0.;
double r;
- while (angle<=2.*M_PI){
+ while(angle<=2.*M_PI){
r = a+b*sin(c*angle);
iso_x.push_back(r*sin(angle)+xmid);
iso_y.push_back(r*cos(angle)+xmid);
@@ -739,13 +739,13 @@ double gLevelsetShamrock::operator()(double x, double y, double z) const
std::vector<double>::const_iterator itminx = iso_x.begin();
std::vector<double>::const_iterator itminy = iso_y.begin();
itx++;ity++;
- for (;itx!=itxen;itx++,ity++){
+ for(;itx!=itxen;itx++,ity++){
xi = *itx;
yi = *ity;
dx = x-xi;
dy = y-yi;
d = sqrt(dx*dx+dy*dy);
- if (distance>d){
+ if(distance>d){
distance = d;
itminx = itx;
itminy = ity;
@@ -756,7 +756,7 @@ double gLevelsetShamrock::operator()(double x, double y, double z) const
SVector3 t,p;
p(0) = (*itminx) - x;
p(1) = (*itminy) - y;
- if (itminx==(itxen-1)){
+ if(itminx==(itxen-1)){
t(0) = iso_x[0] - (*itminx);
t(1) = iso_y[0] - (*itminy);
}
@@ -767,7 +767,7 @@ double gLevelsetShamrock::operator()(double x, double y, double z) const
double vectprod = (p(0)*t(1) - p(1)*t(0));
double sign = (vectprod < 0.) ? -1. : 1.;
- return (sign*distance);
+ return sign * distance;
}
gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r0,
@@ -782,19 +782,19 @@ gLevelsetPopcorn::gLevelsetPopcorn(double _xc, double _yc, double _zc, double _r
zc = _zc;
}
-double gLevelsetPopcorn::operator() (double x, double y, double z) const
+double gLevelsetPopcorn::operator()(double x, double y, double z) const
{
- double s2 = (sigma)*(sigma);
+ double s2 = sigma * sigma;
double r = sqrt((x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc)*(z-zc));
//printf("z=%g zc=%g r=%g \n", z, zc, r);
double val = r - r0;
- for (int k = 0; k< 5; k++){
+ for(int k = 0; k< 5; k++){
double xk = r0/sqrt(5.0)*(2.*cos(2*k*M_PI/5.0));
double yk = r0/sqrt(5.0)*(2.*sin(2*k*M_PI/5.0));
double zk = r0/sqrt(5.0);
val -= A*exp(-((x-xc-xk)*(x-xc-xk)+(y-yc-yk)*(y-yc-yk)+(z-zc-zk)*(z-zc-zk))/s2);
}
- for (int k = 5; k< 10; k++){
+ for(int k = 5; k< 10; k++){
double xk = r0/sqrt(5.0)*(2.*cos((2.*(k-5.)-1.)*M_PI/5.0));
double yk = r0/sqrt(5.0)*(2.*sin((2.*(k-5.)-1.)*M_PI/5.0));
double zk = -r0/sqrt(5.0);
@@ -816,7 +816,7 @@ gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag)
_expr = new mathEvaluator(expressions, variables);
}
-double gLevelsetMathEval::operator() (double x, double y, double z) const
+double gLevelsetMathEval::operator()(double x, double y, double z) const
{
std::vector<double> values(3), res(1);
values[0] = x;
@@ -838,7 +838,7 @@ gLevelsetMathEvalAll::gLevelsetMathEvalAll(std::vector<std::string> expressions,
_expr = new mathEvaluator(expressions, variables);
}
-double gLevelsetMathEvalAll::operator() (double x, double y, double z) const
+double gLevelsetMathEvalAll::operator()(double x, double y, double z) const
{
std::vector<double> values(3), res(13);
values[0] = x;
@@ -891,26 +891,26 @@ gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClo
{
std::map<int, std::vector<GEntity*> > groups [4];
gm->getPhysicalGroups(groups);
- for (GModel::piter it = gm->firstPhysicalName() ;
+ for(GModel::piter it = gm->firstPhysicalName();
it != gm->lastPhysicalName() ; ++it){
- if (it->second == physical){
+ if(it->second == physical){
_entities = groups[it->first.first][it->first.second];
}
}
- if (_entities.size() == 0){
+ if(_entities.size() == 0){
Msg::Error("distanceToMesh: the physical name '%s' does not exist in the GModel",
physical.c_str());
}
//setup
std::set<MVertex *> _all;
- for (unsigned int i=0;i<_entities.size();i++){
- for (unsigned int k = 0; k < _entities[i]->getNumMeshElements(); k++) {
+ for(unsigned int i = 0; i < _entities.size(); i++){
+ for(unsigned int k = 0; k < _entities[i]->getNumMeshElements(); k++) {
MElement *e = _entities[i]->getMeshElement(k);
- for (int j = 0; j< e->getNumVertices();j++){
- MVertex *v = _entities[i]->getMeshElement(k)->getVertex(j);
- _all.insert(v);
- _v2e.insert (std::make_pair(v,e));
+ for(int j = 0; j< e->getNumVertices();j++){
+ MVertex *v = _entities[i]->getMeshElement(k)->getVertex(j);
+ _all.insert(v);
+ _v2e.insert(std::make_pair(v, e));
}
}
}
@@ -918,7 +918,7 @@ gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClo
nodes = annAllocPts(_all.size(), 3);
std::set<MVertex*>::iterator itp = _all.begin();
int ind = 0;
- while (itp != _all.end()){
+ while(itp != _all.end()){
MVertex* pt = *itp;
nodes[ind][0] = pt->x();
nodes[ind][1] = pt->y();
@@ -931,55 +931,54 @@ gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClo
gLevelsetDistMesh::~gLevelsetDistMesh()
{
- if (_kdtree) {
+ if(_kdtree) {
ANNpointArray nodes = _kdtree->thePoints();
annDeallocPts(nodes);
- delete _kdtree;
+ delete _kdtree;
}
-
}
-double gLevelsetDistMesh::operator () (double x, double y, double z) const
+double gLevelsetDistMesh::operator()(double x, double y, double z) const
{
- std::vector<ANNidx> _index(_nbClose);
- std::vector<ANNdist> _dist(_nbClose);
- double point[3] = {x,y,z};
- _kdtree->annkSearch(point, _nbClose, &_index[0], &_dist[0]);
+ std::vector<ANNidx> index(_nbClose);
+ std::vector<ANNdist> dist(_nbClose);
+ double point[3] = {x, y, z};
+ _kdtree->annkSearch(point, _nbClose, &index[0], &dist[0]);
std::set<MElement*> elements;
- for (int i=0;i<_nbClose;i++){
- int iVertex = _index[i];
+ for(int i = 0; i < _nbClose; i++){
+ int iVertex = index[i];
MVertex *v = _vertices[iVertex];
- for (std::multimap<MVertex*,MElement*>::const_iterator itm =
- _v2e.lower_bound(v); itm != _v2e.upper_bound(v); ++itm)
- elements.insert (itm->second);
+ for(std::multimap<MVertex*,MElement*>::const_iterator itm =
+ _v2e.lower_bound(v); itm != _v2e.upper_bound(v); ++itm)
+ elements.insert(itm->second);
}
double minDistance = 1.e22;
SPoint3 closest;
- for (std::set<MElement*>::iterator it = elements.begin();
- it != elements.end();++it){
+ for(std::set<MElement*>::iterator it = elements.begin();
+ it != elements.end(); ++it){
double distance = 0.;
MVertex *v1 = (*it)->getVertex(0);
MVertex *v2 = (*it)->getVertex(1);
SPoint3 p1(v1->x(), v1->y(), v1->z());
SPoint3 p2(v2->x(), v2->y(), v2->z());
SPoint3 pt;
- if ((*it)->getDim() == 1){
- signedDistancePointLine(p1, p2,SPoint3(x,y,z),distance,pt);
+ if((*it)->getDim() == 1){
+ signedDistancePointLine(p1, p2, SPoint3(x, y, z), distance, pt);
}
- else if ((*it)->getDim() == 2){
+ else if((*it)->getDim() == 2){
MVertex *v3 = (*it)->getVertex(2);
SPoint3 p3(v3->x(), v3->y(), v3->z());
- signedDistancePointTriangle(p1, p2, p3,SPoint3(x,y,z),distance,pt);
+ signedDistancePointTriangle(p1, p2, p3, SPoint3(x, y, z), distance, pt);
}
else{
- Msg::Error("Cannot compute a dsitance to an entity of dimension %d",
+ Msg::Error("Cannot compute a distance to an entity of dimension %d\n",
(*it)->getDim());
}
- if (fabs(distance) < fabs(minDistance)){
- minDistance=distance;
+ if(fabs(distance) < fabs(minDistance)){
+ minDistance = distance;
}
}
- return -1.0*minDistance;
+ return -1.0 * minDistance;
}
#endif
@@ -997,7 +996,7 @@ gLevelsetPostView::gLevelsetPostView(int index, int tag)
}
}
-double gLevelsetPostView::operator () (double x, double y, double z) const
+double gLevelsetPostView::operator()(double x, double y, double z) const
{
if(!_octree) return 1.;
double val = 1.;
@@ -1019,7 +1018,7 @@ gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
translate(pt);
}
-gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv)
+gLevelsetGenCylinder::gLevelsetGenCylinder(const gLevelsetGenCylinder& lv)
: gLevelsetQuadric(lv){}
gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir,
@@ -1085,6 +1084,42 @@ gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
children[i] = _children[i]->clone();
}
+gLevelsetYarn::gLevelsetYarn(int dim, int phys, double minA, double majA, int type, int tag)
+ : gLevelsetPrimitive(tag), minorAxis(minA), majorAxis(majA), typeLs(type)
+{
+ std::map<int, std::vector<GEntity*> > groups[4];
+ GModel::current()->getPhysicalGroups(groups);
+ entities = groups[dim][phys];
+ if(!entities.size())
+ printf("No physical %d found for levelset yarn!\n", phys);
+}
+
+double gLevelsetYarn::operator()(double x, double y, double z) const
+{
+ double dist = 0.0;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity *g = entities[i];
+ for(unsigned int j = 0; j < g->getNumMeshElements(); j++){
+ MElement *e = g->getMeshElement(j);
+ MVertex *v1 = e->getVertex(0);
+ MVertex *v2 = e->getVertex(1);
+ SPoint3 p1(v1->x(), v1->y(), v1->z());
+ SPoint3 p2(v2->x(), v2->y(), v2->z());
+ if(e->getType() == TYPE_LIN){
+ signedDistancesPointsEllipseLine(iDistances, iDistancesE, iIsInYarn, iClosePts,
+ pts, p1, p2, majorAxis, minorAxis,
+ majorAxis, minorAxis, typeLs);
+ }
+ else if(e->getType() == TYPE_TRI){
+ MVertex *v3 = e->getVertex(2);
+ SPoint3 p3(v3->x(),v3->y(),v3->z());
+ signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
+ }
+ }
+ }
+ return dist;
+}
+
gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv)
: gLevelset(lv)
{
@@ -1248,7 +1283,6 @@ void gLevelsetNACA00::getClosestBndPoint(double x, double y, double z,
double &xb, double &yb, double &curvRad,
bool &in) const
{
-
static const int maxIter = 100;
static const double tol = 1.e-10;
@@ -1258,7 +1292,7 @@ void gLevelsetNACA00::getClosestBndPoint(double x, double y, double z,
// Point translated according to airfoil origin and symmetry
const double xt = x-_x0, yt = fabs(y-_y0);
- if (xt-_c > 1.21125*_t*yt) {
+ if(xt-_c > 1.21125*_t*yt) {
// Behind line normal to airfoil at trailing edge, closest boundary point is
// trailing edge...
xb = _x0+_c;
@@ -1269,7 +1303,7 @@ void gLevelsetNACA00::getClosestBndPoint(double x, double y, double z,
const double fact = 5.*_t*_c;
double xtb = std::max(xt,tolr), ytb;
double dyb, ddyb;
- for (int it=0; it<maxIter; it++) {
+ for(int it=0; it<maxIter; it++) {
const double xbr = xtb/_c, sxbr = sqrt(xbr), xbr32 = xbr*sxbr,
xbr2 = xbr*xbr, xbr3 = xbr2*xbr, xbr4 = xbr2*xbr2;
ytb = fact*(0.2969*sxbr-0.1260*xbr-0.3516*xbr2+0.2843*xbr3-0.1036*xbr4);
@@ -1280,9 +1314,9 @@ void gLevelsetNACA00::getClosestBndPoint(double x, double y, double z,
const double dDistSq = -2.*(xx+dyb*yy);
const double ddDistSq = 2.*(1.-ddyb*yy+dyb*dyb);
const double mIncr = dDistSq/ddDistSq;
- if (fabs(mIncr) < tolr) break;
+ if(fabs(mIncr) < tolr) break;
else xtb -= mIncr;
- if (xtb < tolr) xtb = tolr;
+ if(xtb < tolr) xtb = tolr;
if (xtb > _c-tolr) xtb = _c-tolr;
}
xb = _x0+xtb;
@@ -1292,7 +1326,7 @@ void gLevelsetNACA00::getClosestBndPoint(double x, double y, double z,
}
}
-double gLevelsetNACA00::operator() (double x, double y, double z) const
+double gLevelsetNACA00::operator()(double x, double y, double z) const
{
double xb, yb, curvRadb;
bool in;
@@ -1303,8 +1337,8 @@ double gLevelsetNACA00::operator() (double x, double y, double z) const
return in ? -sqrt(distSq) : sqrt(distSq);
}
-void gLevelsetNACA00::gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const
+void gLevelsetNACA00::gradient(double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
{
double xb, yb, curvRadb;
bool in;
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 581148c36d..9f5ebbd8fe 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -54,11 +54,11 @@ class ANNkd_tree;
class gLevelset : public simpleFunction<double>
{
-protected:
- // negative values of the levelset are inside the domain.
+ protected:
+ // negative values of the levelset are inside the domain.
static const short insideDomain = -1;
int tag_; // must be greater than 0
-public:
+ public:
gLevelset() : tag_(-1) {}
gLevelset(const gLevelset &);
virtual ~gLevelset(){}
@@ -69,15 +69,15 @@ public:
virtual double operator() (double x, double y, double z) const = 0;
bool isInsideDomain(const double &x, const double &y, const double &z) const
{
- return this->operator()(x,y,z) * insideDomain > 0.;
+ return this->operator()(x, y, z) * insideDomain > 0.;
}
bool isOutsideDomain(const double &x, const double &y, const double &z) const
{
- return this->operator()(x,y,z) * insideDomain < 0.;
+ return this->operator()(x, y, z) * insideDomain < 0.;
}
bool isOnBorder(const double &x, const double &y, const double &z) const
{
- return this->operator()(x,y,z) == 0.;
+ return this->operator()(x, y, z) == 0.;
}
virtual std::vector<gLevelset *> getChildren() const = 0;
virtual double choose (double d1, double d2) const = 0;
@@ -85,12 +85,12 @@ public:
virtual bool isPrimitive() const = 0;
void setTag(int t) { tag_ = t; }
virtual int getTag() const { return tag_; }
- void getPrimitives(std::vector<gLevelset *> &primitives) ;
- void getPrimitivesPO(std::vector<gLevelset *> &primitives) ;
- void getRPN(std::vector<gLevelset *> &gLsRPN) ;
- double H (const double &x, const double &y, const double &z) const
+ void getPrimitives(std::vector<gLevelset *> &primitives);
+ void getPrimitivesPO(std::vector<gLevelset *> &primitives);
+ void getRPN(std::vector<gLevelset *> &gLsRPN);
+ double H(const double &x, const double &y, const double &z) const
{
- if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
+ if(isInsideDomain(x, y, z) || isOnBorder(x, y, z)) return 1.0;
return 0.0;
}
void print() const
@@ -120,24 +120,22 @@ public:
class gLevelsetPrimitive : public gLevelset
{
-public:
+ public:
gLevelsetPrimitive() : gLevelset() {}
gLevelsetPrimitive(const gLevelsetPrimitive &lv) : gLevelset(lv) {}
gLevelsetPrimitive(int tag) {
- if (tag < 1) {
+ if(tag < 1) {
printf("Tag of the levelset (%d) must be greater than 0.\n", tag);
tag = abs(tag);
}
tag_ = tag;
}
- virtual double operator () (double x, double y, double z) const = 0;
- std::vector<gLevelset *> getChildren() const
- {
+ virtual double operator()(double x, double y, double z) const = 0;
+ std::vector<gLevelset *> getChildren() const {
std::vector<gLevelset *> p; return p;
}
- double choose (double d1, double d2) const
- {
- Msg::Error("Cannot use function \"choose\" with a primitive!");
+ double choose(double d1, double d2) const {
+ Msg::Error("Cannot use function \"choose\" with a primitive!\n");
return d1;
}
virtual int type() const = 0;
@@ -146,12 +144,12 @@ public:
class gLevelsetSphere : public gLevelsetPrimitive
{
-protected:
+ protected:
double xc, yc, zc, r;
-public:
- gLevelsetSphere (const double &x, const double &y, const double &z,
- const double &R, int tag=1);
- virtual double operator () (double x, double y, double z) const
+ public:
+ gLevelsetSphere(const double &x, const double &y, const double &z,
+ const double &R, int tag = 1);
+ virtual double operator()(double x, double y, double z) const
{
if(r >= 0.)
return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) +
@@ -159,12 +157,12 @@ public:
return (- r - sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) +
(zc - z) * (zc - z)));
}
- void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const;
- void hessian (double x, double y, double z,
- double & dfdxx, double & dfdxy, double & dfdxz,
- double & dfdyx, double & dfdyy, double & dfdyz,
- double & dfdzx, double & dfdzy, double & dfdzz ) const;
+ void gradient(double x, double y, double z,
+ double &dfdx, double &dfdy, double &dfdz) const;
+ void hessian(double x, double y, double z,
+ double &dfdxx, double &dfdxy, double &dfdxz,
+ double &dfdyx, double &dfdyy, double &dfdyz,
+ double &dfdzx, double &dfdzy, double &dfdzz) const;
int type() const { return SPHERE; }
};
@@ -172,24 +170,24 @@ class gLevelsetPlane : public gLevelsetPrimitive
{
protected:
double a, b, c, d;
-public:
+ public:
// define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
- gLevelsetPlane (const double _a, const double _b, const double _c,
- const double _d, int tag=1)
+ gLevelsetPlane(const double _a, const double _b, const double _c,
+ const double _d, int tag = 1)
: gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
// define the plane passing through the point pt and with outward normal norm
- gLevelsetPlane (const std::vector<double> &pt, const std::vector<double> &norm,
- int tag=1);
- gLevelsetPlane (const double *pt, const double *norm, int tag=1);
+ gLevelsetPlane(const std::vector<double> &pt, const std::vector<double> &norm,
+ int tag = 1);
+ gLevelsetPlane(const double *pt, const double *norm, int tag = 1);
// define the plane passing through the 3 points pt1,pt2,pt3 and with outward
// normal (pt1,pt2)x(pt1,pt3)
- gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3,
- int tag=1);
+ gLevelsetPlane(const double *pt1, const double *pt2, const double *pt3,
+ int tag = 1);
// copy constructor
gLevelsetPlane(const gLevelsetPlane &lv);
- virtual gLevelset * clone() const{ return new gLevelsetPlane(*this); }
+ virtual gLevelset * clone() const { return new gLevelsetPlane(*this); }
// return negative value inward and positive value outward
- virtual double operator() (double x, double y, double z) const
+ virtual double operator()(double x, double y, double z) const
{
return a * x + b * y + c * z + d;
}
@@ -198,99 +196,99 @@ public:
class gLevelsetPoints : public gLevelsetPrimitive
{
-protected:
+ protected:
fullMatrix<double> points;
fullMatrix<double> surf;
fullMatrix<double> matAInv;
double delta;
- std::map<SPoint3,double> mapP;
+ std::map<SPoint3, double> mapP;
fullMatrix<double> generateRbfMat(int p, int index,
- const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2) const;
+ const fullMatrix<double> &nodes1,
+ const fullMatrix<double> &nodes2) const;
void RbfOp(int p, int index,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &D,
- bool isLocal = false) const;
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &nodes,
+ fullMatrix<double> &D,
+ bool isLocal = false) const;
void evalRbfDer(int p, int index,
- const fullMatrix<double> &cntrs,
- const fullMatrix<double> &nodes,
- const fullMatrix<double> &fValues,
- fullMatrix<double> &fApprox, bool isLocal = false) const;
+ const fullMatrix<double> &cntrs,
+ const fullMatrix<double> &nodes,
+ const fullMatrix<double> &fValues,
+ fullMatrix<double> &fApprox, bool isLocal = false) const;
void setup_level_set(const fullMatrix<double> &cntrs,
- fullMatrix<double> &level_set_nodes,
- fullMatrix<double> &level_set_funvals);
+ fullMatrix<double> &level_set_nodes,
+ fullMatrix<double> &level_set_funvals);
-public:
+ public:
// define the data points
- gLevelsetPoints(fullMatrix<double> &_centers, int tag=1);
+ gLevelsetPoints(fullMatrix<double> &_centers, int tag = 1);
// copy constructor
gLevelsetPoints(const gLevelsetPoints &lv);
- virtual gLevelset * clone() const{return new gLevelsetPoints(*this);}
+ virtual gLevelset * clone() const { return new gLevelsetPoints(*this); }
// return negative value inward and positive value outward
- virtual double operator() (double x, double y, double z) const;
+ virtual double operator()(double x, double y, double z) const;
void computeLS(std::vector<MVertex*> &vert);
- int type() const {return LSPOINTS;}
+ int type() const { return LSPOINTS; }
};
class gLevelsetQuadric : public gLevelsetPrimitive
{
-protected:
+ protected:
double A[3][3], B[3], C;
- void translate (const double transl[3]);
- void rotate (const double rotate[3][3]);
- void computeRotationMatrix (const double dir[3], double t[3][3]);
- void computeRotationMatrix (const double dir1[3], const double dir2[3] ,
- double t[3][3]);
- void Ax (const double x[3], double res[3], double fact=1.0);
- void xAx (const double x[3], double &res, double fact=1.0);
- void init ();
-public:
- gLevelsetQuadric(int tag=1) : gLevelsetPrimitive(tag) {init(); }
+ void translate(const double transl[3]);
+ void rotate(const double rotate[3][3]);
+ void computeRotationMatrix(const double dir[3], double t[3][3]);
+ void computeRotationMatrix(const double dir1[3], const double dir2[3],
+ double t[3][3]);
+ void Ax(const double x[3], double res[3], double fact = 1.0);
+ void xAx(const double x[3], double &res, double fact = 1.0);
+ void init();
+ public:
+ gLevelsetQuadric(int tag = 1) : gLevelsetPrimitive(tag) { init(); }
gLevelsetQuadric(const gLevelsetQuadric &);
virtual ~gLevelsetQuadric() {}
- double operator () (double x, double y, double z) const;
+ double operator()(double x, double y, double z) const;
virtual int type() const = 0;
};
class gLevelsetGenCylinder : public gLevelsetQuadric
{
-public:
- gLevelsetGenCylinder (const double *pt, const double *dir, const double &R,
- int tag=1);
- gLevelsetGenCylinder (const gLevelsetGenCylinder& );
- virtual gLevelset * clone() const{ return new gLevelsetGenCylinder(*this); }
- int type() const {return GENCYLINDER;}
+ public:
+ gLevelsetGenCylinder(const double *pt, const double *dir, const double &R,
+ int tag = 1);
+ gLevelsetGenCylinder(const gLevelsetGenCylinder &);
+ virtual gLevelset * clone() const { return new gLevelsetGenCylinder(*this); }
+ int type() const { return GENCYLINDER; }
};
class gLevelsetEllipsoid : public gLevelsetQuadric
{
-public:
- gLevelsetEllipsoid (const double *pt, const double *dir, const double &a,
- const double &b, const double &c, int tag=1);
- gLevelsetEllipsoid (const gLevelsetEllipsoid& );
- virtual gLevelset * clone() const{ return new gLevelsetEllipsoid(*this); }
+ public:
+ gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
+ const double &b, const double &c, int tag = 1);
+ gLevelsetEllipsoid(const gLevelsetEllipsoid &);
+ virtual gLevelset * clone() const { return new gLevelsetEllipsoid(*this); }
int type() const { return ELLIPS; }
};
class gLevelsetCone : public gLevelsetQuadric
{
-public:
- gLevelsetCone (const double *pt, const double *dir, const double &angle, int tag=1);
- gLevelsetCone (const gLevelsetCone& );
- virtual gLevelset * clone() const{ return new gLevelsetCone(*this); }
+ public:
+ gLevelsetCone(const double *pt, const double *dir, const double &angle, int tag = 1);
+ gLevelsetCone(const gLevelsetCone &);
+ virtual gLevelset * clone() const { return new gLevelsetCone(*this); }
int type() const { return CONE; }
};
class gLevelsetGeneralQuadric : public gLevelsetQuadric
{
-public:
- gLevelsetGeneralQuadric (const double *pt, const double *dir, const double &x2,
- const double &y2, const double &z2, const double &z,
- const double &c, int tag=1);
- gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& );
- virtual gLevelset * clone() const{ return new gLevelsetGeneralQuadric(*this); }
- int type() const {return QUADRIC;}
+ public:
+ gLevelsetGeneralQuadric(const double *pt, const double *dir,
+ const double &x2, const double &y2, const double &z2,
+ const double &z, const double &c, int tag = 1);
+ gLevelsetGeneralQuadric(const gLevelsetGeneralQuadric &);
+ virtual gLevelset * clone() const { return new gLevelsetGeneralQuadric(*this); }
+ int type() const { return QUADRIC; }
};
class gLevelsetPopcorn: public gLevelsetPrimitive
@@ -299,11 +297,11 @@ class gLevelsetPopcorn: public gLevelsetPrimitive
double sigma;
double r0;
double xc, yc, zc;
-public:
+ public:
gLevelsetPopcorn(double xc, double yc, double zc, double r0, double A,
- double sigma, int tag=1);
- ~gLevelsetPopcorn(){}
- double operator () (double x, double y, double z) const;
+ double sigma, int tag = 1);
+ ~gLevelsetPopcorn() {}
+ double operator()(double x, double y, double z) const;
int type() const { return UNKNOWN; }
};
@@ -316,51 +314,51 @@ class gLevelsetShamrock: public gLevelsetPrimitive
double xmid, a, b;
int c;
std::vector<double> iso_x, iso_y;
-public:
+ public:
gLevelsetShamrock(double xmid, double ymid, double zmid, double a, double b,
- int c=3, int tag=1);
- ~gLevelsetShamrock(){}
- double operator () (double x, double y, double z) const;
+ int c = 3, int tag = 1);
+ ~gLevelsetShamrock() {}
+ double operator()(double x, double y, double z) const;
int type() const { return UNKNOWN; }
};
class gLevelsetMathEval: public gLevelsetPrimitive
{
mathEvaluator *_expr;
-public:
- gLevelsetMathEval(std::string f, int tag=1);
- ~gLevelsetMathEval(){ if (_expr) delete _expr; }
- double operator () (double x, double y, double z) const;
+ public:
+ gLevelsetMathEval(std::string f, int tag = 1);
+ ~gLevelsetMathEval() { if(_expr) delete _expr; }
+ double operator()(double x, double y, double z) const;
int type() const { return UNKNOWN; }
};
class gLevelsetMathEvalAll: public gLevelsetPrimitive
{
mathEvaluator *_expr;
-public:
- gLevelsetMathEvalAll(std::vector<std::string> f, int tag=1);
- ~gLevelsetMathEvalAll(){ if (_expr) delete _expr; }
- double operator() (double x, double y, double z) const;
- void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const;
- void hessian (double x, double y, double z,
- double & dfdxx, double & dfdxy, double & dfdxz,
- double & dfdyx, double & dfdyy, double & dfdyz,
- double & dfdzx, double & dfdzy, double & dfdzz ) const;
+ public:
+ gLevelsetMathEvalAll(std::vector<std::string> f, int tag = 1);
+ ~gLevelsetMathEvalAll() { if(_expr) delete _expr; }
+ double operator()(double x, double y, double z) const;
+ void gradient(double x, double y, double z,
+ double &dfdx, double &dfdy, double &dfdz) const;
+ void hessian(double x, double y, double z,
+ double &dfdxx, double &dfdxy, double &dfdxz,
+ double &dfdyx, double &dfdyy, double &dfdyz,
+ double &dfdzx, double &dfdzy, double &dfdzz) const;
int type() const { return UNKNOWN; }
};
class gLevelsetSimpleFunction: public gLevelsetPrimitive
{
simpleFunction<double> *_f;
-public:
- gLevelsetSimpleFunction(simpleFunction<double> *f, int tag=1){
+ public:
+ gLevelsetSimpleFunction(simpleFunction<double> *f, int tag = 1) {
_f = f;
}
- ~gLevelsetSimpleFunction(){}
- double operator () (double x, double y, double z) const
+ ~gLevelsetSimpleFunction() {}
+ double operator()(double x, double y, double z) const
{
- return (*_f)(x,y,z);
+ return (*_f)(x, y, z);
}
int type() const { return UNKNOWN; }
};
@@ -371,13 +369,13 @@ class gLevelsetDistMesh: public gLevelsetPrimitive
const int _nbClose;
std::vector<GEntity*> _entities;
std::vector<MVertex*> _vertices;
- std::multimap<MVertex*,MElement*> _v2e;
+ std::multimap<MVertex*, MElement*> _v2e;
ANNkd_tree *_kdtree;
-public :
- gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5, int tag=1);
- double operator () (double x, double y, double z) const;
+ public :
+ gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5, int tag = 1);
+ double operator()(double x, double y, double z) const;
~gLevelsetDistMesh();
- int type() const { return UNKNOWN; }
+ int type() const { return LSMESH; }
};
#endif
@@ -386,14 +384,46 @@ class gLevelsetPostView : public gLevelsetPrimitive
{
int _viewIndex;
OctreePost *_octree;
-public:
- gLevelsetPostView(int index, int tag=1) ;
- ~gLevelsetPostView(){ if(_octree) delete _octree;}
- double operator () (double x, double y, double z) const;
+ public:
+ gLevelsetPostView(int index, int tag = 1);
+ ~gLevelsetPostView() { if(_octree) delete _octree; }
+ double operator()(double x, double y, double z) const;
int type() const { return UNKNOWN; }
};
#endif
+class gLevelsetNACA00 : public gLevelsetPrimitive
+{
+ double _x0, _y0, _c, _t;
+ public:
+ gLevelsetNACA00(double x0, double y0, double c, double t);
+ ~gLevelsetNACA00() {}
+ double operator()(double x, double y, double z) const;
+ void gradient(double x, double y, double z,
+ double &dfdx, double &dfdy, double &dfdz) const;
+ void hessian(double x, double y, double z,
+ double &dfdxx, double &dfdxy, double &dfdxz,
+ double &dfdyx, double &dfdyy, double &dfdyz,
+ double &dfdzx, double &dfdzy, double &dfdzz) const;
+ int type() const { return UNKNOWN; }
+ private:
+ void getClosestBndPoint(const double x, const double y, const double z,
+ double &xb, double &yb, double &curvRad,
+ bool &in) const;
+};
+
+class gLevelsetYarn : public gLevelsetPrimitive
+{
+ double minorAxis, majorAxis;
+ int typeLs;
+ std::vector<GEntity*> entities;
+ public:
+ gLevelsetYarn(int dim, int phys, double minA, double majA, int type, int tag = 1);
+ ~gLevelsetYarn() {}
+ double operator()(double x, double y, double z) const;
+ int type() const { return UNKNOWN; }
+};
+
// TOOLS
class gLevelsetTools : public gLevelset
@@ -401,32 +431,35 @@ class gLevelsetTools : public gLevelset
protected:
std::vector<gLevelset *> children;
bool _delChildren;//flag to delete only if called from gmsh Parser
-public:
- gLevelsetTools () {}
- gLevelsetTools (const std::vector<gLevelset *> &p, bool delC=false)
+ public:
+ gLevelsetTools() {}
+ gLevelsetTools(const std::vector<gLevelset *> &p, bool delC = false)
{
- children = p; _delChildren=delC;
+ children = p; _delChildren = delC;
}
- gLevelsetTools (const gLevelsetTools &);
- virtual ~gLevelsetTools () {
- if (_delChildren){
+ gLevelsetTools(const gLevelsetTools &);
+ virtual ~gLevelsetTools()
+ {
+ if(_delChildren){
for(int i = 0; i < (int)children.size(); i++)
- delete children[i];
+ delete children[i];
}
}
- double operator () (double x, double y, double z) const {
+ double operator()(double x, double y, double z) const
+ {
double d = (*children[0])(x, y, z);
- for (int i = 1; i < (int)children.size(); i++){
+ for(int i = 1; i < (int)children.size(); i++){
double dt = (*children[i])(x, y, z);
d = choose(d, dt);
}
return d;
}
- std::vector<gLevelset *> getChildren() const {
+ std::vector<gLevelset *> getChildren() const
+ {
if(children.size() != 1) return children;
return children[0]->getChildren();
}
- virtual double choose (double d1, double d2) const = 0;
+ virtual double choose(double d1, double d2) const = 0;
virtual int type2() const = 0;
virtual int type() const
{
@@ -447,17 +480,17 @@ public:
class gLevelsetReverse : public gLevelset
{
-protected:
+ protected:
gLevelset *ls;
-public:
- gLevelsetReverse (gLevelset *p) : ls(p){}
- double operator () (double x, double y, double z) const
+ public:
+ gLevelsetReverse(gLevelset *p) : ls(p) {}
+ double operator()(double x, double y, double z) const
{
return -(*ls)(x, y, z);
}
std::vector<gLevelset *> getChildren() const { return ls->getChildren(); }
- bool isPrimitive() const {return ls->isPrimitive();}
- virtual double choose (double d1, double d2) const { return -ls->choose(d1,d2); }
+ bool isPrimitive() const { return ls->isPrimitive(); }
+ virtual double choose(double d1, double d2) const { return -ls->choose(d1, d2); }
virtual int type() const { return ls->type(); }
int getTag() const { return ls->getTag(); }
};
@@ -466,14 +499,14 @@ public:
// others as tools that cut it
class gLevelsetCut : public gLevelsetTools
{
-public:
- gLevelsetCut (std::vector<gLevelset *> p, bool delC=false)
- : gLevelsetTools(p,delC) {}
- double choose (double d1, double d2) const
+ public:
+ gLevelsetCut(std::vector<gLevelset *> p, bool delC = false)
+ : gLevelsetTools(p, delC) {}
+ double choose(double d1, double d2) const
{
return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
}
- gLevelsetCut(const gLevelsetCut &lv) : gLevelsetTools(lv){}
+ gLevelsetCut(const gLevelsetCut &lv) : gLevelsetTools(lv) {}
virtual gLevelset * clone() const { return new gLevelsetCut(*this); }
int type2() const { return CUT; }
};
@@ -481,12 +514,13 @@ public:
// This levelset takes the minimum
class gLevelsetUnion : public gLevelsetTools
{
-public:
- gLevelsetUnion (std::vector<gLevelset *> p, bool delC=false)
- : gLevelsetTools(p,delC) { }
- gLevelsetUnion(const gLevelsetUnion &lv):gLevelsetTools(lv){}
- virtual gLevelset * clone() const{return new gLevelsetUnion(*this);}
- double choose (double d1, double d2) const
+ public:
+ gLevelsetUnion(std::vector<gLevelset *> p, bool delC = false)
+ : gLevelsetTools(p, delC) {}
+ gLevelsetUnion(const gLevelsetUnion &lv) : gLevelsetTools(lv) {}
+ virtual gLevelset * clone() const{ return new gLevelsetUnion(*this); }
+
+ double choose(double d1, double d2) const
{
return (d1 < d2) ? d1 : d2; // lesser of d1 and d2
}
@@ -496,13 +530,13 @@ public:
// This levelset takes the maximum
class gLevelsetIntersection : public gLevelsetTools
{
-public:
- gLevelsetIntersection (std::vector<gLevelset *> p, bool delC=false)
- : gLevelsetTools(p,delC) { }
- gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv) { }
+ public:
+ gLevelsetIntersection(std::vector<gLevelset *> p, bool delC = false)
+ : gLevelsetTools(p, delC) {}
+ gLevelsetIntersection(const gLevelsetIntersection &lv) : gLevelsetTools(lv) {}
virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
- double choose (double d1, double d2) const
- {
+
+ double choose(double d1, double d2) const {
return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
int type2() const { return INTER; }
@@ -511,21 +545,21 @@ public:
// Crack defined by a normal and a tangent levelset
class gLevelsetCrack : public gLevelsetTools
{
-public:
- gLevelsetCrack (std::vector<gLevelset *> p, bool delC=false)
+ public:
+ gLevelsetCrack(std::vector<gLevelset *> p, bool delC = false)
{
- if (p.size() != 2)
+ if(p.size() != 2)
printf("Error : gLevelsetCrack needs 2 levelsets\n");
children.push_back(p[0]);
children.push_back(new gLevelsetReverse(p[0]));
if(p[1]) children.push_back(p[1]);
_delChildren = delC;
}
- double choose (double d1, double d2) const
+ double choose(double d1, double d2) const
{
return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
- int type2() const {return CRACK;}
+ int type2() const { return CRACK; }
};
@@ -533,21 +567,21 @@ public:
class gLevelsetImproved : public gLevelset
{
-protected:
+ protected:
gLevelset *Ls;
-public:
- gLevelsetImproved(){}
+ public:
+ gLevelsetImproved() {}
gLevelsetImproved(const gLevelsetImproved &lv);
- double operator() (double x, double y, double z) const { return (*Ls)(x, y, z); }
+ double operator()(double x, double y, double z) const { return (*Ls)(x, y, z); }
std::vector<gLevelset *> getChildren() const { return Ls->getChildren(); }
- double choose (double d1, double d2) const { return Ls->choose(d1, d2); }
+ double choose(double d1, double d2) const { return Ls->choose(d1, d2); }
virtual int type() const = 0;
- bool isPrimitive() const {return Ls->isPrimitive();}
+ bool isPrimitive() const { return Ls->isPrimitive(); }
};
class gLevelsetBox : public gLevelsetImproved
{
-public:
+ public:
// create a box with parallel faces :
// pt is a corner of the box,
// dir1 is the direction of the first edge starting from pt,
@@ -564,7 +598,7 @@ public:
// face normal to dir1 and including pt : tag+5
gLevelsetBox(const double *pt, const double *dir1, const double *dir2,
const double *dir3, const double &a, const double &b,
- const double &c, int tag=1);
+ const double &c, int tag = 1);
// create a box with the 8 vertices (pt1,...,pt8).
// check if the faces are planar.
// tags of the faces are : face(pt5,pt6,pt7,pt8) : tag+0
@@ -575,15 +609,15 @@ public:
// face(pt1,pt5,pt8,pt4) : tag+5
gLevelsetBox(const double *pt1, const double *pt2, const double *pt3,
const double *pt4, const double *pt5, const double *pt6,
- const double *pt7, const double *pt8, int tag=1);
+ const double *pt7, const double *pt8, int tag = 1);
gLevelsetBox(const gLevelsetBox &);
- virtual gLevelset * clone() const{return new gLevelsetBox(*this);}
+ virtual gLevelset * clone() const { return new gLevelsetBox(*this); }
int type() const { return BOX; }
};
class gLevelsetCylinder : public gLevelsetImproved
{
-public:
+ public:
// create a cylinder : pt is the point in the middle of the cylinder base,
// dir is the direction of the cylinder axis,
// R is the outer radius of the cylinder,
@@ -591,11 +625,11 @@ public:
// tags of the faces are : exterior face : tag+0
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
- gLevelsetCylinder (const std::vector<double> &pt,
- const std::vector<double> &dir, const double &R,
- const double &H, int tag=1);
- gLevelsetCylinder (const double *pt, const double *dir, const double &R,
- const double &H, int tag=1);
+ gLevelsetCylinder(const std::vector<double> &pt,
+ const std::vector<double> &dir, const double &R,
+ const double &H, int tag = 1);
+ gLevelsetCylinder(const double *pt, const double *dir, const double &R,
+ const double &H, int tag = 1);
// create a cylinder : pt is the point in the middle of the cylinder base,
// dir is the direction of the cylinder axis,
// R is the outer radius of the cylinder,
@@ -605,16 +639,16 @@ public:
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
// interior face : tag+3
- gLevelsetCylinder (const double *pt, const double *dir, const double &R,
- const double &r, const double &H, int tag=1);
+ gLevelsetCylinder(const double *pt, const double *dir, const double &R,
+ const double &r, const double &H, int tag = 1);
gLevelsetCylinder(const gLevelsetCylinder &);
- virtual gLevelset * clone() const{return new gLevelsetCylinder(*this);}
+ virtual gLevelset * clone() const { return new gLevelsetCylinder(*this); }
int type() const { return CYLINDER; }
};
class gLevelsetConrod : public gLevelsetImproved
{
-public:
+ public:
// create a connecting rod :
// pt is the point in the middle of the first bore,
// dir1 is the direction of the rod,
@@ -643,24 +677,24 @@ public:
// top face of the second cylinder : tag+11
// interior face of the first cylinder : tag+12
// interior face of the second cylinder : tag+13
- gLevelsetConrod (const double *pt, const double *dir1, const double *dir2,
- const double &H1, const double &H2, const double &H3,
- const double &R1, const double &r1, const double &R2,
- const double &r2, const double &L1, const double &L2,
- const double &E, int tag=1);
+ gLevelsetConrod(const double *pt, const double *dir1, const double *dir2,
+ const double &H1, const double &H2, const double &H3,
+ const double &R1, const double &r1, const double &R2,
+ const double &r2, const double &L1, const double &L2,
+ const double &E, int tag = 1);
gLevelsetConrod(const gLevelsetConrod &);
- virtual gLevelset * clone() const{ return new gLevelsetConrod(*this); }
+ virtual gLevelset * clone() const { return new gLevelsetConrod(*this); }
int type() const { return CONROD; }
};
/*
class gLevelsetDisk : public gLevelsetTools
{
-public:
- // define the disk of given radius centered at a point pt and with outward
- // normal norm
- gLevelsetDisk (std::vector<double> pt, std::vector<double> dir,
- const double R, bool delC=false)
+ public:
+ // define the disk of given radius centered at a point pt and with outward
+ // normal norm
+ gLevelsetDisk(std::vector<double> pt, std::vector<double> dir,
+ const double R, bool delC = false)
{
double *ptP, *normP;
ptP[0] = pt[0]; ptP[1] = pt[1]; ptP[2] = pt[2];
@@ -679,32 +713,12 @@ public:
children.push_back(cyl);
_delChildren = delC;
}
- double choose (double d1, double d2) const
+ double choose(double d1, double d2) const
{
return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
- int type2() const {return DISK;}
+ int type2() const { return DISK; }
};
*/
-class gLevelsetNACA00: public gLevelsetPrimitive
-{
- double _x0, _y0, _c, _t;
-public:
- gLevelsetNACA00(double x0, double y0, double c, double t);
- ~gLevelsetNACA00() {}
- double operator () (double x, double y, double z) const;
- void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const;
- void hessian (double x, double y, double z,
- double & dfdxx, double & dfdxy, double & dfdxz,
- double & dfdyx, double & dfdyy, double & dfdyz,
- double & dfdzx, double & dfdzy, double & dfdzz ) const;
- int type() const { return UNKNOWN; }
-private:
- void getClosestBndPoint(const double x, const double y, const double z,
- double &xb, double &yb, double &curvRad,
- bool &in) const;
-};
-
#endif
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 4748ba728a..aba8d5f933 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -507,7 +507,6 @@ double ComputeScalarRep(int numComp, double *val)
void eigenvalue2x2(double mat[2][2], double v[2])
{
-
double a=1;
double b=-(mat[0][0]+mat[1][1]);
double c= (mat[0][0]*mat[1][1])-(mat[0][1]*mat[1][0]);
@@ -516,7 +515,6 @@ void eigenvalue2x2(double mat[2][2], double v[2])
v[0] = (-b+sqrt(det))/(2*a);
v[1] = (-b-sqrt(det))/(2*a);
-
}
void eigenvalue(double mat[3][3], double v[3])
@@ -694,7 +692,7 @@ bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
dx(0) = f(0) / J(0, 0);
else
if (!J.luSolve(f, dx))
- return false;
+ return false;
for (int i = 0; i < N; i++)
x(i) -= relax * dx(i);
@@ -706,9 +704,7 @@ bool newton_fd(bool (*func)(fullVector<double> &, fullVector<double> &, void *),
/*
min_a f(x+a*d);
-
f(x+a*d) = f(x) + f'(x) (
-
*/
void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
@@ -818,7 +814,6 @@ double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
// printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
if (check == 1) break;
}
-
return f;
}
@@ -851,7 +846,6 @@ void signedDistancePointTriangle(const SPoint3 &p1,const SPoint3 &p2, const SPoi
const SPoint3 &p, double &d, SPoint3 &closePt)
{
-
SVector3 t1 = p2 - p1;
SVector3 t2 = p3 - p1;
SVector3 t3 = p3 - p2;
@@ -908,9 +902,9 @@ void signedDistancePointTriangle(const SPoint3 &p1,const SPoint3 &p2, const SPoi
closePt = p3;
d = sign * std::min(fabs(d), p.distance(p3));
}
- }
-
+ }
}
+
void signedDistancesPointsTriangle(std::vector<double> &distances,
std::vector<SPoint3> &closePts,
const std::vector<SPoint3> &pts,
@@ -997,8 +991,6 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
}
}
-
-
void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p,
double &d, SPoint3 &closePt)
{
@@ -1154,9 +1146,9 @@ int computeDistanceRatio(const double &y, const double &yp, const double &x,
}
else{
if (x == xp){
- y1 = x1;
+ y1 = x1;
y2 = x2;
- x1 = -b;
+ x1 = -b;
x2 = -b;
}
else{
@@ -1169,32 +1161,32 @@ int computeDistanceRatio(const double &y, const double &yp, const double &x,
else{
y1 = -(b + x1) / a;
y2 = -(b + x2) / a;
- }
+ }
}
}
if (x1 == x2){
propdist = (y1 - y) / (yp - y);
if(propdist < 0.0){
- propdist = (y2 - y) / (yp - y);
+ propdist = (y2 - y) / (yp - y);
}
}
else{
if (xp != x){
propdist = (x1 - x) / (xp - x);
- if (propdist < 0.0){
- propdist = (x2 - x) / (xp - x);
- }
+ if (propdist < 0.0){
+ propdist = (x2 - x) / (xp - x);
+ }
}
else{
- if (yp != y){
- propdist = (y1 - y) / (yp - y);
- if(propdist < 0.0){
- propdist = (y2 - y) / (yp - y);
- }
- }
+ if (yp != y){
+ propdist = (y1 - y) / (yp - y);
+ if(propdist < 0.0){
+ propdist = (y2 - y) / (yp - y);
+ }
+ }
else{
- propdist = 0.01;
- }
+ propdist = 0.01;
+ }
}
}
*distance = propdist;
@@ -1202,13 +1194,53 @@ int computeDistanceRatio(const double &y, const double &yp, const double &x,
}
}
+void signedDistancesPointsEllipsePoint(std::vector<double> &distances,
+ std::vector<double> &distancesE,
+ std::vector<int> &isInYarn,
+ std::vector<SPoint3> &closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2,
+ const double radius)
+{
+ distances.clear();
+ distances.resize(pts.size());
+ distancesE.clear();
+ distancesE.resize(pts.size());
+ isInYarn.clear();
+ isInYarn.resize(pts.size());
+ closePts.clear();
+ closePts.resize(pts.size());
+ double d;
+ for (unsigned int i = 0; i < pts.size();i++){
+ SPoint3 closePt;
+ const SPoint3 &p = pts[i];
+ signedDistancePointLine(p1, p2, p, d, closePt);
+ closePts[i] = closePt;
+ distances[i] = d;
+ if (d <= radius){
+ isInYarn[i] = 1;
+ distancesE[i] = radius - d;
+ }
+ else{
+ isInYarn[i] = 0;
+ distancesE[i] = d - radius;
+ }
+ }
+}
+
void signedDistancesPointsEllipseLine(std::vector<double>&distances,
std::vector<double> &distancesE,
std::vector<int>&isInYarn,
std::vector<SPoint3>&closePts,
const std::vector<SPoint3> &pts,
const SPoint3 &p1,
- const SPoint3 &p2)
+ const SPoint3 &p2,
+ const double maxA,
+ const double minA,
+ const double maxB,
+ const double minB,
+ const int typeLevelSet)
{
distances.clear();
distances.resize(pts.size());
@@ -1229,37 +1261,127 @@ void signedDistancesPointsEllipseLine(std::vector<double>&distances,
int direction=0;
if (!(p.x()==closePt.x() && p.y()==closePt.y() && p.z()==closePt.z())){
double xp,yp,x,y,otherp,other,propdist;
- if (p1.x()==p2.x()){
- direction=1;
- if (fabs(closePt.x() - 0.0) < 0.00000001) isInYarn[i] = 1;
- if (fabs(closePt.x() - 2.2) < 0.00000001) isInYarn[i] = 4;
- if (fabs(closePt.x() - 4.4) < 0.00000001) isInYarn[i] = 2;
- if (fabs(closePt.x() - 6.6) < 0.00000001) isInYarn[i] = 5;
- if (fabs(closePt.x() - 8.8) < 0.00000001) isInYarn[i] = 3;
- if (fabs(closePt.x() - 11.0) < 0.00000001) isInYarn[i] = 1;
+ if (typeLevelSet==2){
+ if (p1.x()==p2.x()){
+ direction=1;
+ if (fabs(closePt.x()-0.0)<0.00000001) isInYarn[i]=1;
+ if (fabs(closePt.x()-2.2)<0.00000001) isInYarn[i]=1;
+ }
+ else{
+ if (p1.y()==p2.y()){
+ direction=2;
+ if (fabs(closePt.y()-0.0)<0.00000001) isInYarn[i]=6;
+ if (fabs(closePt.y()-2.2)<0.00000001) isInYarn[i]=6;
+ }
+ else{
+ printf("Error %lf %lf\n",closePt.x(),closePt.y());
+ }
+ }
}
- else{
- if (p1.y() == p2.y()){
- direction = 2;
- if (fabs(closePt.y() - 0.0) < 0.00000001) isInYarn[i] = 6;
- if (fabs(closePt.y() - 2.2) < 0.00000001) isInYarn[i] = 7;
- if (fabs(closePt.y() - 4.4) < 0.00000001) isInYarn[i] = 8;
- if (fabs(closePt.y() - 6.6) < 0.00000001) isInYarn[i] = 9;
- if (fabs(closePt.y() - 8.8) < 0.00000001) isInYarn[i] = 10;
- if (fabs(closePt.y() - 11.0) < 0.00000001) isInYarn[i] = 6;
+ if (typeLevelSet==1){
+ if (p1.x()==p2.x()){
+ direction=1;
+ //if (fabs(closePt.x() - 0.0) < 0.00000001) isInYarn[i] = 1;
+ if (fabs(closePt.x() - 2.2) < 0.00000001) isInYarn[i] = 4;
+ //if (fabs(closePt.x() - 4.4) < 0.00000001) isInYarn[i] = 2;
+ //if (fabs(closePt.x() - 6.6) < 0.00000001) isInYarn[i] = 5;
+ //if (fabs(closePt.x() - 8.8) < 0.00000001) isInYarn[i] = 3;
+ //if (fabs(closePt.x() - 11.0) < 0.00000001) isInYarn[i] = 1;
+ }
+ else{
+ if (p1.y() == p2.y()){
+ direction = 2;
+ //if (fabs(closePt.y() - 0.0) < 0.00000001) isInYarn[i] = 6;
+ if (fabs(closePt.y() - 2.2) < 0.00000001) isInYarn[i] = 7;
+ //if (fabs(closePt.y() - 4.4) < 0.00000001) isInYarn[i] = 8;
+ //if (fabs(closePt.y() - 6.6) < 0.00000001) isInYarn[i] = 9;
+ //if (fabs(closePt.y() - 8.8) < 0.00000001) isInYarn[i] = 10;
+ //if (fabs(closePt.y() - 11.0) < 0.00000001) isInYarn[i] = 6;
+ }
+ else{
+ printf("Error %lf %lf\n", closePt.x(), closePt.y());
+ }
+ }
+ }
+ if (typeLevelSet==4){
+ if (p1.x()==p2.x()){
+ direction=1;
+ if (fabs(closePt.x()-0.0)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=1;
+ if (fabs(closePt.x()-2.2)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=4;
+ if (fabs(closePt.x()-4.4)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=2;
+ if (fabs(closePt.x()-6.6)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=5;
+ if (fabs(closePt.x()-8.8)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=3;
+ if (fabs(closePt.x()-11.0)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=1;
+ if (fabs(closePt.x()-0.0)<0.00000001 && closePt.z()>0.35) isInYarn[i]=11;
+ if (fabs(closePt.x()-2.2)<0.00000001 && closePt.z()>0.35) isInYarn[i]=14;
+ if (fabs(closePt.x()-4.4)<0.00000001 && closePt.z()>0.35) isInYarn[i]=12;
+ if (fabs(closePt.x()-6.6)<0.00000001 && closePt.z()>0.35) isInYarn[i]=15;
+ if (fabs(closePt.x()-8.8)<0.00000001 && closePt.z()>0.35) isInYarn[i]=13;
+ if (fabs(closePt.x()-11.0)<0.00000001 && closePt.z()>0.35) isInYarn[i]=11;
+ }
+ else{
+ if (p1.y()==p2.y()){
+ direction=2;
+ if (fabs(closePt.y()-0.0)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=6;
+ if (fabs(closePt.y()-2.2)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=7;
+ if (fabs(closePt.y()-4.4)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=8;
+ if (fabs(closePt.y()-6.6)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=9;
+ if (fabs(closePt.y()-8.8)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=10;
+ if (fabs(closePt.y()-11.0)<0.00000001 && closePt.z()<=0.35) isInYarn[i]=6;
+ if (fabs(closePt.y()-0.0)<0.00000001 && closePt.z()>0.35) isInYarn[i]=16;
+ if (fabs(closePt.y()-2.2)<0.00000001 && closePt.z()>0.35) isInYarn[i]=17;
+ if (fabs(closePt.y()-4.4)<0.00000001 && closePt.z()>0.35) isInYarn[i]=18;
+ if (fabs(closePt.y()-6.6)<0.00000001 && closePt.z()>0.35) isInYarn[i]=19;
+ if (fabs(closePt.y()-8.8)<0.00000001 && closePt.z()>0.35) isInYarn[i]=20;
+ if (fabs(closePt.y()-11.0)<0.00000001 && closePt.z()>0.35) isInYarn[i]=16;
+ }
+ else{
+ printf("Error %lf %lf\n",closePt.x(),closePt.y());
+ }
+ }
+ }
+ if (typeLevelSet==3){
+ direction=3;
+ isInYarn[i]=1;
+ }
+ if (typeLevelSet==5){
+ if(p1.x()==p2.x()){
+ direction=1;
+ if (fabs(closePt.x() - 0.0) < 0.00000001) isInYarn[i] = 1;
+ if (fabs(closePt.x() - 3.225) < 0.00000001) isInYarn[i] = 2;
+ if (fabs(closePt.x() - 6.45) < 0.00000001) isInYarn[i] = 3;
+ if (fabs(closePt.x() - 9.675) < 0.00000001) isInYarn[i] = 4;
+ if (fabs(closePt.x() - 12.9) < 0.00000001) isInYarn[i] = 1;
}
else{
- printf("WTF %lf %lf\n", closePt.x(), closePt.y());
+ if (p1.y()==p2.y()){
+ direction=2;
+ if (fabs(closePt.y() - 0.0) < 0.00000001) isInYarn[i] = 5;
+ if (fabs(closePt.y() - 1.665) < 0.00000001) isInYarn[i] = 6;
+ if (fabs(closePt.y() - 3.33) < 0.00000001) isInYarn[i] = 7;
+ if (fabs(closePt.y() - 4.995) < 0.00000001) isInYarn[i] = 8;
+ if (fabs(closePt.y() - 6.66) < 0.00000001) isInYarn[i] = 5;
+ }
+ else{
+ printf("Error %lf %lf\n",closePt.x(),closePt.y());
+ }
}
}
changeReferential(direction, p, closePt, p1, p2, &xp, &yp,
&otherp, &x, &y, &other);
- int result;
+ int result = 1;
if (fabs(other-otherp) > 0.01){
- result = 1;
+ result = 1;
}
else{
- result = computeDistanceRatio(y, yp, x, xp, &propdist, 1.1, 0.0875);
+ if (direction==1){
+ result = computeDistanceRatio(y, yp, x, xp, &propdist, maxA, minA);
+ }
+ else{
+ if (direction==2){
+ result = computeDistanceRatio(y, yp, x, xp, &propdist, maxB, minB);
+ }
+ }
}
if (result == 1){
distancesE[i] = 1.e10;
@@ -1271,7 +1393,7 @@ void signedDistancesPointsEllipseLine(std::vector<double>&distances,
distancesE[i] = (1.0 / propdist) - 1.0;
}
else{
- distancesE[i] = (1.0 - (1.0 / propdist)) / 3.0;
+ distancesE[i] = (1.0 - (1.0 / propdist));
}
}
}
@@ -1308,7 +1430,7 @@ int intersection_segments(const SPoint3 &p1, const SPoint3 &p2,
double b[2] = {q1.x() - p1.x(), q1.y() - p1.y()};
sys2x2(A, b, x);
return (x[0] >= 0.0 && x[0] <= 1. &&
- x[1] >= 0.0 && x[1] <= 1.);
+ x[1] >= 0.0 && x[1] <= 1.);
}
}
@@ -1420,7 +1542,7 @@ void projectPointsToPlane(const std::vector<SPoint3> &pts, std::vector<SPoint3>
void transformPointsIntoOrthoBasis(const std::vector<SPoint3> &ptsProj,
std::vector<SPoint3> &pointsUV,
- const SPoint3 &ptCG, const mean_plane &meanPlane)
+ const SPoint3 &ptCG, const mean_plane &meanPlane)
{
pointsUV.resize(ptsProj.size());
SVector3 normal(meanPlane.a, meanPlane.b, meanPlane.c);
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index bb059023c3..00a9b80c65 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -135,12 +135,23 @@ int computeDistanceRatio(const double &y, const double &yp, const double &x,
const double &xp, double *distance, const double &r1,
const double &r2);
+void signedDistancesPointsEllipsePoint (std::vector<double>&distances,
+ std::vector<double>&distancesE,
+ std::vector<int>&isInYarn,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1, const SPoint3 &p2,
+ const double radius);
+
void signedDistancesPointsEllipseLine (std::vector<double>&distances,
std::vector<double>&distancesE,
std::vector<int>&isInYarn,
std::vector<SPoint3>&closePts,
const std::vector<SPoint3> &pts,
- const SPoint3 &p1, const SPoint3 &p2);
+ const SPoint3 &p1, const SPoint3 &p2,
+ const double maxA, const double minA,
+ const double maxB, const double minB,
+ const int typeLevelSet);
int intersection_segments (const SPoint3 &p1, const SPoint3 &p2,
const SPoint3 &q1, const SPoint3 &q2,
--
GitLab