1 #ifndef FILE_SURFACE 2 #define FILE_SURFACE 3 4 /**************************************************************************/ 5 /* File: surface.hh */ 6 /* Author: Joachim Schoeberl */ 7 /* Date: 1. Dez. 95 */ 8 /**************************************************************************/ 9 10 11 namespace netgen 12 { 13 14 class TriangleApproximation; 15 16 17 /** 18 Basis class for implicit surface geometry. 19 This class is used for generation of surface meshes 20 in NETGEN 21 */ 22 class Surface 23 { 24 protected: 25 /// invert normal vector 26 bool inverse; 27 /// maximal h in surface 28 double maxh; 29 /// name of surface 30 char * name; 31 /// boundary condition nr 32 int bcprop; 33 /// boundary condition label 34 string bcname; 35 36 public: 37 Surface (); 38 /** @name Tangential plane. 39 The tangential plane is used for surface mesh generation. 40 */ 41 42 virtual ~Surface(); 43 44 protected: 45 /** @name Points in the surface defining tangential plane. 46 Tangential plane is taken in p1, the local x-axis 47 is directed to p2. 48 */ 49 //@{ 50 /// 51 Point<3> p1; 52 /// 53 Point<3> p2; 54 //@} 55 /** @name Base-vectos for local coordinate system. */ 56 //@{ 57 /// in plane, directed p1->p2 58 Vec<3> ex; 59 /// in plane 60 Vec<3> ey; 61 /// outer normal direction 62 Vec<3> ez; 63 //@} 64 public: 65 66 void SetName (const char * aname); Name() const67 const char * Name () const { return name; } 68 69 //@{ 70 /** 71 Defines tangential plane in ap1. 72 The local x-coordinate axis points to the direction of ap2 */ 73 virtual void DefineTangentialPlane (const Point<3> & ap1, 74 const Point<3> & ap2); 75 76 /// Transforms 3d point p3d to local coordinates pplane 77 virtual void ToPlane (const Point<3> & p3d, Point<2> & pplane, 78 double h, int & zone) const; 79 80 /// Transforms point pplane in local coordinates to 3d point 81 virtual void FromPlane (const Point<2> & pplane, 82 Point<3> & p3d, double h) const; 83 //@} 84 85 86 /// Project point p onto surface (closest point) 87 virtual void Project (Point<3> & p) const; 88 89 /// Project along direction 90 virtual void SkewProject(Point<3> & p, const Vec<3> & direction) const; 91 92 /// Is current surface identic to surface 2 ? IsIdentic(const Surface &,int &,double) const93 virtual int IsIdentic (const Surface & /* s2 */, int & /* inv */, 94 double /* eps */) const 95 { return 0; } 96 97 /// 98 virtual int PointOnSurface (const Point<3> & p, 99 double eps = 1e-6) const; 100 101 102 /** @name Implicit function. 103 Calculate function value and derivatives. 104 */ 105 //@{ 106 /// Calculate implicit function value in point point 107 virtual double CalcFunctionValue (const Point<3> & point) const = 0; 108 109 /** 110 Calc gradient of implicit function. 111 gradient should be O(1) at surface 112 */ 113 virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const = 0; 114 115 /** 116 Calculate second derivatives of implicit function. 117 */ 118 virtual void CalcHesse (const Point<3> & point, Mat<3> & hesse) const; 119 120 /** 121 Returns outer normal vector. 122 */ 123 // virtual void GetNormalVector (const Point<3> & p, Vec<3> & n) const; 124 virtual Vec<3> GetNormalVector (const Point<3> & p) const; 125 126 /** 127 Upper bound for spectral norm of Hesse-matrix 128 */ 129 virtual double HesseNorm () const = 0; 130 131 /** 132 Upper bound for spectral norm of Hesse-matrix in the 133 rad - environment of point c. 134 */ HesseNormLoc(const Point<3> &,double) const135 virtual double HesseNormLoc (const Point<3> & /* c */, 136 double /* rad */) const 137 { return HesseNorm (); } 138 //@} 139 140 141 /// 142 virtual double MaxCurvature () const; 143 /// 144 virtual double MaxCurvatureLoc (const Point<3> & /* c */ , 145 double /* rad */) const; 146 147 /** Returns any point in the surface. 148 Needed to start surface mesh generation e.g. on sphere */ 149 virtual Point<3> GetSurfacePoint () const = 0; 150 151 /// Inverse() const152 bool Inverse () const { return inverse; } 153 /// SetInverse(bool ainverse)154 void SetInverse (bool ainverse) { inverse = ainverse; } 155 /// 156 virtual void Print (ostream & str) const = 0; 157 158 /// Reduce(const BoxSphere<3> &)159 virtual void Reduce (const BoxSphere<3> & /* box */) { }; 160 /// UnReduce()161 virtual void UnReduce () { }; 162 163 /// set max h in surface SetMaxH(double amaxh)164 void SetMaxH (double amaxh) { maxh = amaxh; } 165 /// GetMaxH() const166 double GetMaxH () const { return maxh; } 167 /// GetBCProperty() const168 int GetBCProperty () const { return bcprop; } 169 /// SetBCProperty(int abc)170 void SetBCProperty (int abc) { bcprop = abc; } 171 172 /** Determine local mesh-size. 173 Find 174 \[ h \leq hmax, \] 175 such that 176 \[ h \times \kappa (x) \leq c \qquad \mbox{in} B(x, h), \] 177 where kappa(x) is the curvature in x. */ 178 virtual double LocH (const Point<3> & p, double x, 179 double c, double hmax) const; 180 181 /** 182 Gets Approximation by triangles, 183 where qual is about the number of triangles per radius 184 */ GetTriangleApproximation(TriangleApproximation &,const Box<3> &,double) const185 virtual void GetTriangleApproximation (TriangleApproximation & /* tas */, 186 const Box<3> & /* boundingbox */, 187 double /* facets */ ) const { }; 188 189 GetBCName() const190 string GetBCName() const { return bcname; } 191 SetBCName(string abc)192 void SetBCName( string abc ) { bcname = abc; } 193 }; 194 195 operator <<(ostream & ost,const Surface & surf)196 inline ostream & operator<< (ostream & ost, const Surface & surf) 197 { 198 surf.Print(ost); 199 return ost; 200 } 201 202 203 204 typedef enum { IS_OUTSIDE = 0, IS_INSIDE = 1, DOES_INTERSECT = 2} 205 INSOLID_TYPE; 206 207 208 209 210 class DummySurface : public Surface 211 { CalcFunctionValue(const Point<3> &) const212 virtual double CalcFunctionValue (const Point<3> & /* point */) const 213 { return 0; } 214 CalcGradient(const Point<3> &,Vec<3> & grad) const215 virtual void CalcGradient (const Point<3> & /* point */, Vec<3> & grad) const 216 { grad = Vec<3> (0,0,0); } 217 GetSurfacePoint() const218 virtual Point<3> GetSurfacePoint () const 219 { return Point<3> (0,0,0); } 220 HesseNorm() const221 virtual double HesseNorm () const 222 { return 0; } 223 Project(Point<3> &) const224 virtual void Project (Point<3> & /* p */) const 225 { ; } 226 Print(ostream & ost) const227 virtual void Print (ostream & ost) const 228 { ost << "dummy surface"; } 229 }; 230 231 232 233 class Primitive 234 { 235 236 public: 237 238 Primitive (); 239 240 virtual ~Primitive(); 241 242 243 /* 244 Check, whether box intersects solid defined by surface. 245 246 return values: 247 0 .. box outside solid \\ 248 1 .. box in solid \\ 249 2 .. can't decide (allowed, iff box is close to solid) 250 */ 251 virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const = 0; 252 virtual INSOLID_TYPE PointInSolid (const Point<3> & p, 253 double eps) const = 0; 254 255 virtual void GetTangentialSurfaceIndices (const Point<3> & p, 256 Array<int> & surfind, double eps) const; 257 258 virtual INSOLID_TYPE VecInSolid (const Point<3> & p, 259 const Vec<3> & v, 260 double eps) const = 0; 261 262 // checks if lim s->0 lim t->0 p + t(v1 + s v2) in solid 263 virtual INSOLID_TYPE VecInSolid2 (const Point<3> & p, 264 const Vec<3> & v1, 265 const Vec<3> & v2, 266 double eps) const; 267 268 // checks if p + s v1 + s*s/2 v2 is inside 269 virtual INSOLID_TYPE VecInSolid3 (const Point<3> & p, 270 const Vec<3> & v1, 271 const Vec<3> & v2, 272 double eps) const; 273 274 // like VecInSolid2, but second order approximation 275 virtual INSOLID_TYPE VecInSolid4 (const Point<3> & p, 276 const Vec<3> & v, 277 const Vec<3> & v2, 278 const Vec<3> & m, 279 double eps) const; 280 281 virtual void GetTangentialVecSurfaceIndices (const Point<3> & p, const Vec<3> & v, 282 Array<int> & surfind, double eps) const; 283 284 virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2, 285 Array<int> & surfind, double eps) const; 286 287 CalcSpecialPoints(Array<Point<3>> &) const288 virtual void CalcSpecialPoints (Array<Point<3> > & /* pts */) const { ; } AnalyzeSpecialPoint(const Point<3> &,Array<Point<3>> &) const289 virtual void AnalyzeSpecialPoint (const Point<3> & /* pt */, 290 Array<Point<3> > & /* specpts */) const { ; } SpecialPointTangentialVector(const Point<3> &,int,int) const291 virtual Vec<3> SpecialPointTangentialVector (const Point<3> & /* p */, 292 int /* s1 */, int /* s2 */) const 293 { return Vec<3> (0,0,0); } 294 295 296 virtual int GetNSurfaces() const = 0; 297 virtual Surface & GetSurface (int i = 0) = 0; 298 virtual const Surface & GetSurface (int i = 0) const = 0; 299 300 Array<int> surfaceids; 301 Array<int> surfaceactive; 302 303 int GetSurfaceId (int i = 0) const; 304 void SetSurfaceId (int i, int id); SurfaceActive(int i) const305 int SurfaceActive (int i) const { return surfaceactive[i]; } SurfaceInverted(int=0) const306 virtual int SurfaceInverted (int /* i */ = 0) const { return 0; } 307 308 virtual void GetPrimitiveData (const char *& classname, 309 Array<double> & coeffs) const; 310 virtual void SetPrimitiveData (Array<double> & coeffs); 311 static Primitive * CreatePrimitive (const char * classname); 312 313 Reduce(const BoxSphere<3> &)314 virtual void Reduce (const BoxSphere<3> & /* box */) { }; UnReduce()315 virtual void UnReduce () { }; 316 317 virtual Primitive * Copy () const; 318 virtual void Transform (Transformation<3> & trans); 319 }; 320 321 322 323 324 class OneSurfacePrimitive : public Surface, public Primitive 325 { 326 public: 327 OneSurfacePrimitive(); 328 ~OneSurfacePrimitive(); 329 330 virtual INSOLID_TYPE PointInSolid (const Point<3> & p, 331 double eps) const; 332 virtual INSOLID_TYPE VecInSolid (const Point<3> & p, 333 const Vec<3> & v, 334 double eps) const; 335 virtual INSOLID_TYPE VecInSolid2 (const Point<3> & p, 336 const Vec<3> & v1, 337 const Vec<3> & v2, 338 double eps) const; 339 340 virtual INSOLID_TYPE VecInSolid3 (const Point<3> & p, 341 const Vec<3> & v1, 342 const Vec<3> & v2, 343 double eps) const; 344 345 virtual INSOLID_TYPE VecInSolid4 (const Point<3> & p, 346 const Vec<3> & v, 347 const Vec<3> & v2, 348 const Vec<3> & m, 349 double eps) const; 350 351 virtual int GetNSurfaces() const; 352 virtual Surface & GetSurface (int i = 0); 353 virtual const Surface & GetSurface (int i = 0) const; 354 }; 355 356 357 358 359 360 361 /** 362 Projects point to edge. 363 The point hp is projected to the edge descibed by f1 and f2. 364 It is assumed that the edge is non-degenerated, and the 365 (generalized) Newton method converges. 366 */ 367 extern void ProjectToEdge (const Surface * f1, 368 const Surface * f2, 369 Point<3> & hp); 370 371 372 } 373 374 #endif 375