1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributor(s):
7 //   Gaetan Bricteux
8 
9 #ifndef GMSH_LEVELSET_H
10 #define GMSH_LEVELSET_H
11 
12 #include <string>
13 #include <cmath>
14 #include <stdio.h>
15 #include <stdlib.h> // for abs()
16 #include <vector>
17 #include "GmshMessage.h"
18 #include "fullMatrix.h"
19 #include "GModel.h"
20 #include "MVertex.h"
21 #include "GmshConfig.h"
22 #include "cartesian.h"
23 #include "simpleFunction.h"
24 
25 #if defined(HAVE_ANN)
26 class ANNkd_tree;
27 #endif
28 
29 #if defined(HAVE_POST)
30 #include "PView.h"
31 #include "OctreePost.h"
32 #endif
33 
34 class mathEvaluator;
35 
36 // PRIMITIVE LEVELSET
37 #define UNKNOWN 0
38 #define SPHERE 1
39 #define PLANE 2
40 #define GENCYLINDER 3
41 #define ELLIPS 4
42 #define CONE 5
43 #define QUADRIC 6
44 #define BOX 7
45 #define CYLINDER 8
46 #define CONROD 9
47 #define LSMESH 10
48 #define LSPOINTS 11 // don't define 'POINTS' as it's reserved by win32
49 
50 // TOOLS
51 #define CUT 12
52 #define UNION 13
53 #define INTER 14
54 #define CRACK 15
55 #define DISK 16
56 
57 class gLevelset;
58 
59 class gLevelsetLessThan {
60 public:
61   bool operator()(const gLevelset *l1, const gLevelset *l2) const;
62 };
63 
64 class gLevelset : public simpleFunction<double> {
65 protected:
66   // negative values of the levelset are inside the domain.
67   static const short insideDomain = -1;
68   // unique levelset id, must be greater than 0
69   int tag_;
70   // max tag in all levelsets
71   static int maxTag_;
72   // all levelsets
73   static std::set<gLevelset *, gLevelsetLessThan> all_;
74 
75 public:
76   gLevelset(int tag = 0)
77   {
78     if(tag <= 0)
79       tag_ = gLevelset::maxTag_++;
80     else
81       tag_ = tag;
82   }
83   gLevelset(const gLevelset &);
~gLevelset()84   virtual ~gLevelset() {}
85   static gLevelset *find(int tag);
86   static void add(gLevelset *l);
clone()87   virtual gLevelset *clone() const { return nullptr; }
operator()88   virtual double operator()(double x, double y, double z) const { return 0.; }
isInsideDomain(const double & x,const double & y,const double & z)89   bool isInsideDomain(const double &x, const double &y, const double &z) const
90   {
91     return this->operator()(x, y, z) * insideDomain > 0.;
92   }
isOutsideDomain(const double & x,const double & y,const double & z)93   bool isOutsideDomain(const double &x, const double &y, const double &z) const
94   {
95     return this->operator()(x, y, z) * insideDomain < 0.;
96   }
isOnBorder(const double & x,const double & y,const double & z)97   bool isOnBorder(const double &x, const double &y, const double &z) const
98   {
99     return this->operator()(x, y, z) == 0.;
100   }
getChildren()101   virtual std::vector<gLevelset *> getChildren() const
102   {
103     return std::vector<gLevelset *>();
104   }
choose(double d1,double d2)105   virtual double choose(double d1, double d2) const { return 0.; }
type()106   virtual int type() const { return 0; }
isPrimitive()107   virtual bool isPrimitive() const { return false; }
setTag(int t)108   void setTag(int t) { tag_ = t; }
getTag()109   virtual int getTag() const { return tag_; }
110   void getPrimitives(std::vector<gLevelset *> &primitives);
111   void getPrimitivesPO(std::vector<gLevelset *> &primitives);
112   void getRPN(std::vector<gLevelset *> &gLsRPN);
H(const double & x,const double & y,const double & z)113   double H(const double &x, const double &y, const double &z) const
114   {
115     if(isInsideDomain(x, y, z) || isOnBorder(x, y, z)) return 1.0;
116     return 0.0;
117   }
print()118   void print() const
119   {
120     printf("LS : ");
121     switch(type()) {
122     case SPHERE: printf("SPHERE"); break;
123     case PLANE: printf("PLANE"); break;
124     case GENCYLINDER: printf("GENCYLINDER"); break;
125     case ELLIPS: printf("ELLIPS"); break;
126     case CONE: printf("CONE"); break;
127     case QUADRIC: printf("QUADRIC"); break;
128     case BOX: printf("BOX"); break;
129     case CYLINDER: printf("CYLINDER"); break;
130     case CONROD: printf("CONROD"); break;
131     case CUT: printf("CUT"); break;
132     case UNION: printf("UNION"); break;
133     case INTER: printf("INTER"); break;
134     case LSMESH: printf("LSMESH"); break;
135     case LSPOINTS: printf("LSPOINTS"); break;
136     }
137     printf(" Tag=%d\n", getTag());
138   }
139 };
140 
141 // PRIMITIVES
142 
143 class gLevelsetPrimitive : public gLevelset {
144 public:
gLevelsetPrimitive()145   gLevelsetPrimitive() : gLevelset() {}
gLevelsetPrimitive(const gLevelsetPrimitive & lv)146   gLevelsetPrimitive(const gLevelsetPrimitive &lv) : gLevelset(lv) {}
gLevelsetPrimitive(int tag)147   gLevelsetPrimitive(int tag) : gLevelset(tag) {}
148   virtual double operator()(double x, double y, double z) const = 0;
getChildren()149   std::vector<gLevelset *> getChildren() const
150   {
151     std::vector<gLevelset *> p;
152     return p;
153   }
choose(double d1,double d2)154   double choose(double d1, double d2) const
155   {
156     Msg::Error("Cannot use function \"choose\" with a primitive!\n");
157     return d1;
158   }
159   virtual int type() const = 0;
isPrimitive()160   virtual bool isPrimitive() const { return true; }
161 };
162 
163 class gLevelsetSphere : public gLevelsetPrimitive {
164 protected:
165   double xc, yc, zc, r;
166 
167 public:
168   gLevelsetSphere(const double &x, const double &y, const double &z,
169                   const double &R, int tag = 0);
operator()170   virtual double operator()(double x, double y, double z) const
171   {
172     if(r >= 0.)
173       return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) +
174                   (zc - z) * (zc - z)) -
175              r;
176     return (-r - sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) +
177                       (zc - z) * (zc - z)));
178   }
179   void gradient(double x, double y, double z, double &dfdx, double &dfdy,
180                 double &dfdz) const;
181   void hessian(double x, double y, double z, double &dfdxx, double &dfdxy,
182                double &dfdxz, double &dfdyx, double &dfdyy, double &dfdyz,
183                double &dfdzx, double &dfdzy, double &dfdzz) const;
type()184   int type() const { return SPHERE; }
185 };
186 
187 class gLevelsetPlane : public gLevelsetPrimitive {
188 protected:
189   double a, b, c, d;
190 
191 public:
192   // define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
193   gLevelsetPlane(const double _a, const double _b, const double _c,
194                  const double _d, int tag = 0)
gLevelsetPrimitive(tag)195     : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d)
196   {
197   }
198   // define the plane passing through the point pt and with outward normal norm
199   gLevelsetPlane(const std::vector<double> &pt, const std::vector<double> &norm,
200                  int tag = 0);
201   gLevelsetPlane(const double *pt, const double *norm, int tag = 0);
202   // define the plane passing through the 3 points pt1,pt2,pt3 and with outward
203   // normal (pt1,pt2)x(pt1,pt3)
204   gLevelsetPlane(const double *pt1, const double *pt2, const double *pt3,
205                  int tag = 0);
206   // copy constructor
207   gLevelsetPlane(const gLevelsetPlane &lv);
clone()208   virtual gLevelset *clone() const { return new gLevelsetPlane(*this); }
209   // return negative value inward and positive value outward
operator()210   virtual double operator()(double x, double y, double z) const
211   {
212     return a * x + b * y + c * z + d;
213   }
type()214   int type() const { return PLANE; }
215 };
216 
217 class gLevelsetPoints : public gLevelsetPrimitive {
218 protected:
219   fullMatrix<double> points;
220   fullMatrix<double> surf;
221   fullMatrix<double> matAInv;
222   double delta;
223   std::map<SPoint3, double> mapP;
224   fullMatrix<double> generateRbfMat(int p, int index,
225                                     const fullMatrix<double> &nodes1,
226                                     const fullMatrix<double> &nodes2) const;
227   void RbfOp(int p, int index, const fullMatrix<double> &cntrs,
228              const fullMatrix<double> &nodes, fullMatrix<double> &D,
229              bool isLocal = false) const;
230   void evalRbfDer(int p, int index, const fullMatrix<double> &cntrs,
231                   const fullMatrix<double> &nodes,
232                   const fullMatrix<double> &fValues,
233                   fullMatrix<double> &fApprox, bool isLocal = false) const;
234   void setup_level_set(const fullMatrix<double> &cntrs,
235                        fullMatrix<double> &level_set_nodes,
236                        fullMatrix<double> &level_set_funvals);
237 
238 public:
239   // define the data points
240   gLevelsetPoints(fullMatrix<double> &_centers, int tag = 0);
241   // copy constructor
242   gLevelsetPoints(const gLevelsetPoints &lv);
clone()243   virtual gLevelset *clone() const { return new gLevelsetPoints(*this); }
244   // return negative value inward and positive value outward
245   virtual double operator()(double x, double y, double z) const;
246   void computeLS(std::vector<MVertex *> &vert);
type()247   int type() const { return LSPOINTS; }
248 };
249 
250 class gLevelsetQuadric : public gLevelsetPrimitive {
251 protected:
252   double A[3][3], B[3], C;
253   void translate(const double transl[3]);
254   void rotate(const double rotate[3][3]);
255   void computeRotationMatrix(const double dir[3], double t[3][3]);
256   void computeRotationMatrix(const double dir1[3], const double dir2[3],
257                              double t[3][3]);
258   void Ax(const double x[3], double res[3], double fact = 1.0);
259   void xAx(const double x[3], double &res, double fact = 1.0);
260   void init();
261 
262 public:
gLevelsetPrimitive(tag)263   gLevelsetQuadric(int tag = 0) : gLevelsetPrimitive(tag) { init(); }
264   gLevelsetQuadric(const gLevelsetQuadric &);
~gLevelsetQuadric()265   virtual ~gLevelsetQuadric() {}
266   double operator()(double x, double y, double z) const;
267   virtual int type() const = 0;
268 };
269 
270 class gLevelsetGenCylinder : public gLevelsetQuadric {
271 public:
272   gLevelsetGenCylinder(const double *pt, const double *dir, const double &R,
273                        int tag = 0);
274   gLevelsetGenCylinder(const gLevelsetGenCylinder &);
clone()275   virtual gLevelset *clone() const { return new gLevelsetGenCylinder(*this); }
type()276   int type() const { return GENCYLINDER; }
277 };
278 
279 class gLevelsetEllipsoid : public gLevelsetQuadric {
280 public:
281   gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
282                      const double &b, const double &c, int tag = 0);
283   gLevelsetEllipsoid(const gLevelsetEllipsoid &);
clone()284   virtual gLevelset *clone() const { return new gLevelsetEllipsoid(*this); }
type()285   int type() const { return ELLIPS; }
286 };
287 
288 class gLevelsetCone : public gLevelsetQuadric {
289 public:
290   gLevelsetCone(const double *pt, const double *dir, const double &angle,
291                 int tag = 0);
292   gLevelsetCone(const gLevelsetCone &);
clone()293   virtual gLevelset *clone() const { return new gLevelsetCone(*this); }
type()294   int type() const { return CONE; }
295 };
296 
297 class gLevelsetGeneralQuadric : public gLevelsetQuadric {
298 public:
299   gLevelsetGeneralQuadric(const double *pt, const double *dir, const double &x2,
300                           const double &y2, const double &z2, const double &z,
301                           const double &c, int tag = 0);
302   gLevelsetGeneralQuadric(const gLevelsetGeneralQuadric &);
clone()303   virtual gLevelset *clone() const
304   {
305     return new gLevelsetGeneralQuadric(*this);
306   }
type()307   int type() const { return QUADRIC; }
308 };
309 
310 class gLevelsetPopcorn : public gLevelsetPrimitive {
311   double A;
312   double sigma;
313   double r0;
314   double xc, yc, zc;
315 
316 public:
317   gLevelsetPopcorn(double xc, double yc, double zc, double r0, double A,
318                    double sigma, int tag = 0);
~gLevelsetPopcorn()319   ~gLevelsetPopcorn() {}
320   double operator()(double x, double y, double z) const;
type()321   int type() const { return UNKNOWN; }
322 };
323 
324 // creates the 2D (-approximate- signed distance !) level set corresponding to
325 // the "shamrock-like" iso-zero from Dobrzynski and Frey, "Anisotropic delaunay
326 // mesh adaptation for unsteady simulations", 17th International Meshing
327 // Rountable (2008)(177–194)
328 class gLevelsetShamrock : public gLevelsetPrimitive {
329   double xmid, a, b;
330   int c;
331   std::vector<double> iso_x, iso_y;
332 
333 public:
334   gLevelsetShamrock(double xmid, double ymid, double zmid, double a, double b,
335                     int c = 3, int tag = 0);
~gLevelsetShamrock()336   ~gLevelsetShamrock() {}
337   double operator()(double x, double y, double z) const;
type()338   int type() const { return UNKNOWN; }
339 };
340 
341 class gLevelsetMathEval : public gLevelsetPrimitive {
342   mathEvaluator *_expr;
343 
344 public:
345   gLevelsetMathEval(const std::string &f, int tag = 0);
346   ~gLevelsetMathEval();
347   double operator()(double x, double y, double z) const;
type()348   int type() const { return UNKNOWN; }
349 };
350 
351 class gLevelsetMathEvalAll : public gLevelsetPrimitive {
352   mathEvaluator *_expr;
353 
354 public:
355   gLevelsetMathEvalAll(std::vector<std::string> f, int tag = 0);
356   ~gLevelsetMathEvalAll();
357   double operator()(double x, double y, double z) const;
358   void gradient(double x, double y, double z, double &dfdx, double &dfdy,
359                 double &dfdz) const;
360   void hessian(double x, double y, double z, double &dfdxx, double &dfdxy,
361                double &dfdxz, double &dfdyx, double &dfdyy, double &dfdyz,
362                double &dfdzx, double &dfdzy, double &dfdzz) const;
type()363   int type() const { return UNKNOWN; }
364 };
365 
366 class gLevelsetSimpleFunction : public gLevelsetPrimitive {
367   simpleFunction<double> *_f;
368 
369 public:
370   gLevelsetSimpleFunction(simpleFunction<double> *f, int tag = 0) { _f = f; }
~gLevelsetSimpleFunction()371   ~gLevelsetSimpleFunction() {}
operator()372   double operator()(double x, double y, double z) const
373   {
374     return (*_f)(x, y, z);
375   }
type()376   int type() const { return UNKNOWN; }
377 };
378 
379 #if defined(HAVE_ANN)
380 class gLevelsetDistMesh : public gLevelsetPrimitive {
381   const int _nbClose;
382   std::vector<GEntity *> _entities;
383   std::vector<MVertex *> _vertices;
384   std::multimap<MVertex *, MElement *> _v2e;
385   ANNkd_tree *_kdtree;
386 
387 public:
388   gLevelsetDistMesh(GModel *gm, const std::string &physical, int nbClose = 5,
389                     int tag = 0);
390   double operator()(double x, double y, double z) const;
391   ~gLevelsetDistMesh();
type()392   int type() const { return LSMESH; }
393 };
394 #endif
395 
396 #if defined(HAVE_POST)
397 class gLevelsetPostView : public gLevelsetPrimitive {
398   int _viewIndex;
399   OctreePost *_octree;
400 
401 public:
402   gLevelsetPostView(int index, int tag = 0);
~gLevelsetPostView()403   ~gLevelsetPostView()
404   {
405     if(_octree) delete _octree;
406   }
407   double operator()(double x, double y, double z) const;
type()408   int type() const { return UNKNOWN; }
409 };
410 #endif
411 
412 class gLevelsetNACA00 : public gLevelsetPrimitive {
413   double _x0, _y0, _c, _t;
414 
415 public:
416   gLevelsetNACA00(double x0, double y0, double c, double t);
~gLevelsetNACA00()417   ~gLevelsetNACA00() {}
418   double operator()(double x, double y, double z) const;
419   void gradient(double x, double y, double z, double &dfdx, double &dfdy,
420                 double &dfdz) const;
421   void hessian(double x, double y, double z, double &dfdxx, double &dfdxy,
422                double &dfdxz, double &dfdyx, double &dfdyy, double &dfdyz,
423                double &dfdzx, double &dfdzy, double &dfdzz) const;
type()424   int type() const { return UNKNOWN; }
425 
426 private:
427   void getClosestBndPoint(const double x, const double y, const double z,
428                           double &xb, double &yb, double &curvRad,
429                           bool &in) const;
430 };
431 
432 class gLevelsetYarn : public gLevelsetPrimitive {
433   // double minorAxis, majorAxis;
434   // int typeLs;
435   std::vector<GEntity *> entities;
436 
437 public:
438   gLevelsetYarn(int dim, int phys, double minA, double majA, int type,
439                 int tag = 0);
~gLevelsetYarn()440   ~gLevelsetYarn() {}
441   double operator()(double x, double y, double z) const;
type()442   int type() const { return UNKNOWN; }
443 };
444 
445 // TOOLS
446 
447 class gLevelsetTools : public gLevelset {
448 protected:
449   std::vector<gLevelset *> children;
450   bool _delChildren; // flag to delete only if called from gmsh Parser
451 public:
gLevelset(tag)452   gLevelsetTools(int tag = 0) : gLevelset(tag) {}
453   gLevelsetTools(const std::vector<gLevelset *> &p, bool delC = false,
454                  int tag = 0)
gLevelset(tag)455     : gLevelset(tag)
456   {
457     children = p;
458     _delChildren = delC;
459   }
460   gLevelsetTools(const gLevelsetTools &);
~gLevelsetTools()461   virtual ~gLevelsetTools()
462   {
463     if(_delChildren) {
464       for(int i = 0; i < (int)children.size(); i++) delete children[i];
465     }
466   }
operator()467   double operator()(double x, double y, double z) const
468   {
469     double d = (*children[0])(x, y, z);
470     for(int i = 1; i < (int)children.size(); i++) {
471       double dt = (*children[i])(x, y, z);
472       d = choose(d, dt);
473     }
474     return d;
475   }
getChildren()476   std::vector<gLevelset *> getChildren() const
477   {
478     if(children.size() != 1) return children;
479     return children[0]->getChildren();
480   }
481   virtual double choose(double d1, double d2) const = 0;
482   virtual int type2() const = 0;
type()483   virtual int type() const
484   {
485     if(children.size() != 1) return type2();
486     return children[0]->type();
487   }
isPrimitive()488   virtual bool isPrimitive() const
489   {
490     if(children.size() != 1) return false;
491     return children[0]->isPrimitive();
492   }
getTag()493   int getTag() const
494   {
495     if(children.size() != 1) return tag_;
496     return children[0]->getTag();
497   }
498 };
499 
500 class gLevelsetReverse : public gLevelset {
501 protected:
502   gLevelset *ls;
503 
504 public:
gLevelset(tag)505   gLevelsetReverse(gLevelset *p, int tag = 0) : gLevelset(tag), ls(p) {}
operator()506   double operator()(double x, double y, double z) const
507   {
508     return -(*ls)(x, y, z);
509   }
getChildren()510   std::vector<gLevelset *> getChildren() const { return ls->getChildren(); }
isPrimitive()511   virtual bool isPrimitive() const { return ls->isPrimitive(); }
choose(double d1,double d2)512   virtual double choose(double d1, double d2) const
513   {
514     return -ls->choose(d1, d2);
515   }
type()516   virtual int type() const { return ls->type(); }
getTag()517   int getTag() const { return ls->getTag(); }
518 };
519 
520 // This levelset takes the first levelset in the list as the object and the
521 // others as tools that cut it
522 class gLevelsetCut : public gLevelsetTools {
523 public:
524   gLevelsetCut(const std::vector<gLevelset *> &p, bool delC = false,
525                int tag = 0)
gLevelsetTools(p,delC,tag)526     : gLevelsetTools(p, delC, tag)
527   {
528   }
choose(double d1,double d2)529   double choose(double d1, double d2) const
530   {
531     return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
532   }
gLevelsetCut(const gLevelsetCut & lv)533   gLevelsetCut(const gLevelsetCut &lv) : gLevelsetTools(lv) {}
clone()534   virtual gLevelset *clone() const { return new gLevelsetCut(*this); }
type2()535   int type2() const { return CUT; }
536 };
537 
538 // This levelset takes the minimum
539 class gLevelsetUnion : public gLevelsetTools {
540 public:
541   gLevelsetUnion(const std::vector<gLevelset *> &p, bool delC = false,
542                  int tag = 0)
gLevelsetTools(p,delC,tag)543     : gLevelsetTools(p, delC, tag)
544   {
545   }
gLevelsetUnion(const gLevelsetUnion & lv)546   gLevelsetUnion(const gLevelsetUnion &lv) : gLevelsetTools(lv) {}
clone()547   virtual gLevelset *clone() const { return new gLevelsetUnion(*this); }
548 
choose(double d1,double d2)549   double choose(double d1, double d2) const
550   {
551     return (d1 < d2) ? d1 : d2; // lesser of d1 and d2
552   }
type2()553   int type2() const { return UNION; }
554 };
555 
556 // This levelset takes the maximum
557 class gLevelsetIntersection : public gLevelsetTools {
558 public:
559   gLevelsetIntersection(const std::vector<gLevelset *> &p, bool delC = false,
560                         int tag = 0)
gLevelsetTools(p,delC,tag)561     : gLevelsetTools(p, delC, tag)
562   {
563   }
gLevelsetIntersection(const gLevelsetIntersection & lv)564   gLevelsetIntersection(const gLevelsetIntersection &lv) : gLevelsetTools(lv) {}
clone()565   virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
566 
choose(double d1,double d2)567   double choose(double d1, double d2) const
568   {
569     return (d1 > d2) ? d1 : d2; // greater of d1 and d2
570   }
type2()571   int type2() const { return INTER; }
572 };
573 
574 // Crack defined by a normal and a tangent levelset
575 class gLevelsetCrack : public gLevelsetTools {
576 public:
577   gLevelsetCrack(std::vector<gLevelset *> p, bool delC = false, int tag = 0)
gLevelsetTools(tag)578     : gLevelsetTools(tag)
579   {
580     if(p.size() != 2) printf("Error : gLevelsetCrack needs 2 levelsets\n");
581     children.push_back(p[0]);
582     children.push_back(new gLevelsetReverse(p[0]));
583     if(p[1]) children.push_back(p[1]);
584     _delChildren = delC;
585   }
choose(double d1,double d2)586   double choose(double d1, double d2) const
587   {
588     return (d1 > d2) ? d1 : d2; // greater of d1 and d2
589   }
type2()590   int type2() const { return CRACK; }
591 };
592 
593 // IMPROVED LEVELSET
594 
595 class gLevelsetImproved : public gLevelset {
596 protected:
597   gLevelset *Ls;
598 
599 public:
gLevelset(tag)600   gLevelsetImproved(int tag = 0) : gLevelset(tag) {}
601   gLevelsetImproved(const gLevelsetImproved &lv);
operator()602   double operator()(double x, double y, double z) const
603   {
604     return (*Ls)(x, y, z);
605   }
getChildren()606   std::vector<gLevelset *> getChildren() const { return Ls->getChildren(); }
choose(double d1,double d2)607   double choose(double d1, double d2) const { return Ls->choose(d1, d2); }
608   virtual int type() const = 0;
isPrimitive()609   virtual bool isPrimitive() const { return Ls->isPrimitive(); }
610 };
611 
612 class gLevelsetBox : public gLevelsetImproved {
613 public:
614   // create a box with parallel faces :
615   //    pt is a corner of the box,
616   //    dir1 is the direction of the first edge starting from pt,
617   //    dir2 is the direction of the second edge starting from pt,
618   //    dir3 is the direction of the third edge starting from pt,
619   //    a is the length of the first edge starting from pt,
620   //    b is the length of the second edge starting from pt,
621   //    c is the length of the third edge starting from pt.
622   // tags of the faces are : face normal to dir3 and not including pt : tag+0
623   //                         face normal to dir3 and     including pt : tag+1
624   //                         face normal to dir2 and     including pt : tag+2
625   //                         face normal to dir2 and not including pt : tag+3
626   //                         face normal to dir1 and not including pt : tag+4
627   //                         face normal to dir1 and     including pt : tag+5
628   gLevelsetBox(const double *pt, const double *dir1, const double *dir2,
629                const double *dir3, const double &a, const double &b,
630                const double &c, int tag = 0);
631   // create a box with the 8 vertices (pt1,...,pt8).
632   // check if the faces are planar.
633   // tags of the faces are : face(pt5,pt6,pt7,pt8) : tag+0
634   //                         face(pt1,pt4,pt3,pt2) : tag+1
635   //                         face(pt1,pt2,pt6,pt5) : tag+2
636   //                         face(pt3,pt4,pt8,pt7) : tag+3
637   //                         face(pt2,pt3,pt7,pt6) : tag+4
638   //                         face(pt1,pt5,pt8,pt4) : tag+5
639   gLevelsetBox(const double *pt1, const double *pt2, const double *pt3,
640                const double *pt4, const double *pt5, const double *pt6,
641                const double *pt7, const double *pt8, int tag = 0);
642   gLevelsetBox(const gLevelsetBox &);
clone()643   virtual gLevelset *clone() const { return new gLevelsetBox(*this); }
type()644   int type() const { return BOX; }
645 };
646 
647 class gLevelsetCylinder : public gLevelsetImproved {
648 public:
649   // create a cylinder : pt is the point in the middle of the cylinder base,
650   //                     dir is the direction of the cylinder axis,
651   //                     R is the outer radius of the cylinder,
652   //                     H is the height of the cylinder.
653   // tags of the faces are : exterior face :             tag+0
654   //                         plane face including pt :   tag+1
655   //                         plane face opposite to pt : tag+2
656   gLevelsetCylinder(const std::vector<double> &pt,
657                     const std::vector<double> &dir, const double &R,
658                     const double &H, int tag = 0);
659   gLevelsetCylinder(const double *pt, const double *dir, const double &R,
660                     const double &H, int tag = 0);
661   // create a cylinder : pt is the point in the middle of the cylinder base,
662   //                     dir is the direction of the cylinder axis,
663   //                     R is the outer radius of the cylinder,
664   //                     r is the inner radius of the cylinder,
665   //                     H is the height of the cylinder.
666   // tags of the faces are : exterior face :             tag+0
667   //                         plane face including pt :   tag+1
668   //                         plane face opposite to pt : tag+2
669   //                         interior face :             tag+3
670   gLevelsetCylinder(const double *pt, const double *dir, const double &R,
671                     const double &r, const double &H, int tag = 0);
672   gLevelsetCylinder(const gLevelsetCylinder &);
clone()673   virtual gLevelset *clone() const { return new gLevelsetCylinder(*this); }
type()674   int type() const { return CYLINDER; }
675 };
676 
677 class gLevelsetConrod : public gLevelsetImproved {
678 public:
679   // create a connecting rod :
680   //    pt is the point in the middle of the first bore,
681   //    dir1 is the direction of the rod,
682   //    dir2 is the direction of the axis of the bore,
683   //    H1 is height of the first cylinder,
684   //    H2 is the height of the second cylinder,
685   //    H3 is the height of the rod,
686   //    R1 is the outer radius of the first cylinder,
687   //    r1 is the inner radius of the first cylinder,
688   //    R2 is the outer radius of the second cylinder,
689   //    r2 is the inner radius of the second cylinder,
690   //    L1 is the width of the rod in the plane passing through the middle
691   //       of the first bore,
692   //    L2 is the width of the rod in the plane passing through the middle
693   //       of the second bore,
694   //    E is the distance between the axis of the cylinders.
695   // tags of the faces are : bottom face (+dir2) of the bore :      tag+2
696   //                         top    face (-dir2) of the bore :      tag+3
697   //                         rear   face (-dir1xdir2) of the bore : tag+4
698   //                         front  face (+dir1xdir2) of the bore : tag+5
699   //                         exterior face of the first cylinder :  tag+6
700   //                         bottom   face of the first cylinder :  tag+7
701   //                         top      face of the first cylinder :  tag+8
702   //                         exterior face of the second cylinder : tag+9
703   //                         bottom   face of the second cylinder : tag+10
704   //                         top      face of the second cylinder : tag+11
705   //                         interior face of the first  cylinder : tag+12
706   //                         interior face of the second cylinder : tag+13
707   gLevelsetConrod(const double *pt, const double *dir1, const double *dir2,
708                   const double &H1, const double &H2, const double &H3,
709                   const double &R1, const double &r1, const double &R2,
710                   const double &r2, const double &L1, const double &L2,
711                   const double &E, int tag = 0);
712   gLevelsetConrod(const gLevelsetConrod &);
clone()713   virtual gLevelset *clone() const { return new gLevelsetConrod(*this); }
type()714   int type() const { return CONROD; }
715 };
716 
717 #endif
718