1 #include <mystdlib.h>
2
3 #include <linalg.hpp>
4 #include <csg.hpp>
5
6 namespace netgen
7 {
8
Parallelogram3d(Point<3> ap1,Point<3> ap2,Point<3> ap3)9 Parallelogram3d :: Parallelogram3d (Point<3> ap1, Point<3> ap2, Point<3> ap3)
10 {
11 p1 = ap1;
12 p2 = ap2;
13 p3 = ap3;
14
15 CalcData();
16 }
17
~Parallelogram3d()18 Parallelogram3d ::~Parallelogram3d ()
19 {
20 ;
21 }
22
SetPoints(Point<3> ap1,Point<3> ap2,Point<3> ap3)23 void Parallelogram3d :: SetPoints (Point<3> ap1,
24 Point<3> ap2,
25 Point<3> ap3)
26 {
27 p1 = ap1;
28 p2 = ap2;
29 p3 = ap3;
30
31 CalcData();
32 }
33
CalcData()34 void Parallelogram3d :: CalcData()
35 {
36 v12 = p2 - p1;
37 v13 = p3 - p1;
38 p4 = p2 + v13;
39
40 n = Cross (v12, v13);
41 n.Normalize();
42 }
43
44 int Parallelogram3d ::
IsIdentic(const Surface & s2,int & inv,double eps) const45 IsIdentic (const Surface & s2, int & inv, double eps) const
46 {
47 int id =
48 (fabs (s2.CalcFunctionValue (p1)) <= eps) &&
49 (fabs (s2.CalcFunctionValue (p2)) <= eps) &&
50 (fabs (s2.CalcFunctionValue (p3)) <= eps);
51
52 if (id)
53 {
54 Vec<3> n2;
55 n2 = s2.GetNormalVector(p1);
56 inv = (n * n2) < 0;
57 }
58 return id;
59 }
60
61
CalcFunctionValue(const Point<3> & point) const62 double Parallelogram3d :: CalcFunctionValue (const Point<3> & point) const
63 {
64 return n * (point - p1);
65 }
66
CalcGradient(const Point<3> & point,Vec<3> & grad) const67 void Parallelogram3d :: CalcGradient (const Point<3> & point,
68 Vec<3> & grad) const
69 {
70 grad = n;
71 }
72
CalcHesse(const Point<3> & point,Mat<3> & hesse) const73 void Parallelogram3d :: CalcHesse (const Point<3> & point, Mat<3> & hesse) const
74 {
75 hesse = 0;
76 }
77
HesseNorm() const78 double Parallelogram3d :: HesseNorm () const
79 {
80 return 0;
81 }
82
GetSurfacePoint() const83 Point<3> Parallelogram3d :: GetSurfacePoint () const
84 {
85 return p1;
86 }
87
Print(ostream & str) const88 void Parallelogram3d :: Print (ostream & str) const
89 {
90 str << "Parallelogram3d " << p1 << " - " << p2 << " - " << p3 << endl;
91 }
92
93
94 void Parallelogram3d ::
GetTriangleApproximation(TriangleApproximation & tas,const Box<3> & bbox,double facets) const95 GetTriangleApproximation (TriangleApproximation & tas,
96 const Box<3> & bbox,
97 double facets) const
98 {
99 tas.AddPoint (p1);
100 tas.AddPoint (p2);
101 tas.AddPoint (p3);
102 tas.AddPoint (p4);
103 tas.AddTriangle (TATriangle (0, 0, 1, 2));
104 tas.AddTriangle (TATriangle (0, 2, 1, 3));
105 }
106
107
108
109
110
111
112
113
114
115
Brick(Point<3> ap1,Point<3> ap2,Point<3> ap3,Point<3> ap4)116 Brick :: Brick (Point<3> ap1, Point<3> ap2,
117 Point<3> ap3, Point<3> ap4)
118 {
119 faces.SetSize (6);
120 surfaceids.SetSize (6);
121 surfaceactive.SetSize(6);
122
123 p1 = ap1; p2 = ap2;
124 p3 = ap3; p4 = ap4;
125
126 for (int i = 0; i < 6; i++)
127 {
128 faces[i] = new Plane (Point<3>(0,0,0), Vec<3> (0,0,1));
129 surfaceactive[i] = 1;
130 }
131
132 CalcData();
133 }
134
~Brick()135 Brick :: ~Brick ()
136 {
137 for (int i = 0; i < 6; i++)
138 delete faces[i];
139 }
140
CreateDefault()141 Primitive * Brick :: CreateDefault ()
142 {
143 return new Brick (Point<3> (0,0,0),
144 Point<3> (1,0,0),
145 Point<3> (0,1,0),
146 Point<3> (0,0,1));
147 }
148
149
150
Copy() const151 Primitive * Brick :: Copy () const
152 {
153 return new Brick (p1, p2, p3, p4);
154 }
155
Transform(Transformation<3> & trans)156 void Brick :: Transform (Transformation<3> & trans)
157 {
158 trans.Transform (p1);
159 trans.Transform (p2);
160 trans.Transform (p3);
161 trans.Transform (p4);
162
163 CalcData();
164 }
165
166
167
168
169
170
171
172
173
BoxInSolid(const BoxSphere<3> & box) const174 INSOLID_TYPE Brick :: BoxInSolid (const BoxSphere<3> & box) const
175 {
176 /*
177 int i;
178 double maxval;
179 for (i = 1; i <= 6; i++)
180 {
181 double val = faces.Get(i)->CalcFunctionValue (box.Center());
182 if (i == 1 || val > maxval)
183 maxval = val;
184 }
185
186 if (maxval > box.Diam()) return IS_OUTSIDE;
187 if (maxval < -box.Diam()) return IS_INSIDE;
188 return DOES_INTERSECT;
189 */
190
191 bool inside = 1;
192 bool outside = 0;
193
194 Point<3> p[8];
195 for (int j = 0; j < 8; j++)
196 p[j] = box.GetPointNr(j);
197
198 for (int i = 0; i < 6; i++)
199 {
200 bool outsidei = 1;
201 for (int j = 0; j < 8; j++)
202 {
203 // Point<3> p = box.GetPointNr (j);
204 double val = faces[i]->Plane::CalcFunctionValue (p[j]);
205
206 if (val > 0) inside = 0;
207 if (val < 0) outsidei = 0;
208 }
209 if (outsidei) outside = 1;
210 }
211
212 if (outside) return IS_OUTSIDE;
213 if (inside) return IS_INSIDE;
214 return DOES_INTERSECT;
215 }
216
PointInSolid(const Point<3> & p,double eps) const217 INSOLID_TYPE Brick :: PointInSolid (const Point<3> & p,
218 double eps) const
219 {
220 double maxval = faces[0] -> Plane::CalcFunctionValue (p);
221 for (int i = 1; i < 6; i++)
222 {
223 double val = faces[i] -> Plane::CalcFunctionValue (p);
224 if (val > maxval) maxval = val;
225 }
226
227 if (maxval > eps) return IS_OUTSIDE;
228 if (maxval < -eps) return IS_INSIDE;
229 return DOES_INTERSECT;
230 }
231
232
VecInSolid(const Point<3> & p,const Vec<3> & v,double eps) const233 INSOLID_TYPE Brick :: VecInSolid (const Point<3> & p,
234 const Vec<3> & v,
235 double eps) const
236 {
237 INSOLID_TYPE result = IS_INSIDE;
238 for (int i = 0; i < faces.Size(); i++)
239 {
240 INSOLID_TYPE hres = faces[i]->VecInSolid(p, v, eps);
241 if (hres == IS_OUTSIDE || result == IS_OUTSIDE) result = IS_OUTSIDE;
242 else if (hres == DOES_INTERSECT || result == DOES_INTERSECT) result = DOES_INTERSECT;
243 else result = IS_INSIDE;
244 }
245 return result;
246
247 /*
248 INSOLID_TYPE is = IS_INSIDE;
249 Vec<3> grad;
250 double scal;
251
252 for (int i = 0; i < faces.Size(); i++)
253 {
254 if (faces[i] -> PointOnSurface (p, eps))
255 {
256 GetSurface(i).CalcGradient (p, grad);
257 scal = v * grad;
258
259 if (scal >= eps)
260 is = IS_OUTSIDE;
261 if (scal >= -eps && is == IS_INSIDE)
262 is = DOES_INTERSECT;
263 }
264 }
265 return is;
266 */
267
268 /*
269 Point<3> p2 = p + 1e-2 * v;
270 return PointInSolid (p2, eps);
271 */
272 }
273
274
275
276
277
VecInSolid2(const Point<3> & p,const Vec<3> & v1,const Vec<3> & v2,double eps) const278 INSOLID_TYPE Brick :: VecInSolid2 (const Point<3> & p,
279 const Vec<3> & v1,
280 const Vec<3> & v2,
281 double eps) const
282 {
283 INSOLID_TYPE result = IS_INSIDE;
284 for (int i = 0; i < faces.Size(); i++)
285 {
286 INSOLID_TYPE hres = faces[i]->VecInSolid2(p, v1, v2, eps);
287 if (hres == IS_OUTSIDE || result == IS_OUTSIDE) result = IS_OUTSIDE;
288 else if (hres == DOES_INTERSECT || result == DOES_INTERSECT) result = DOES_INTERSECT;
289 else result = IS_INSIDE;
290 }
291 return result;
292 }
293
VecInSolid3(const Point<3> & p,const Vec<3> & v1,const Vec<3> & v2,double eps) const294 INSOLID_TYPE Brick :: VecInSolid3 (const Point<3> & p,
295 const Vec<3> & v1,
296 const Vec<3> & v2,
297 double eps) const
298 {
299 INSOLID_TYPE result = IS_INSIDE;
300 for (int i = 0; i < faces.Size(); i++)
301 {
302 INSOLID_TYPE hres = faces[i]->VecInSolid3(p, v1, v2, eps);
303 if (hres == IS_OUTSIDE || result == IS_OUTSIDE) result = IS_OUTSIDE;
304 else if (hres == DOES_INTERSECT || result == DOES_INTERSECT) result = DOES_INTERSECT;
305 else result = IS_INSIDE;
306 }
307 return result;
308 }
309
VecInSolid4(const Point<3> & p,const Vec<3> & v,const Vec<3> & v2,const Vec<3> & m,double eps) const310 INSOLID_TYPE Brick :: VecInSolid4 (const Point<3> & p,
311 const Vec<3> & v,
312 const Vec<3> & v2,
313 const Vec<3> & m,
314 double eps) const
315 {
316 INSOLID_TYPE result = IS_INSIDE;
317 for (int i = 0; i < faces.Size(); i++)
318 {
319 INSOLID_TYPE hres = faces[i]->VecInSolid4(p, v, v2, m, eps);
320 if (hres == IS_OUTSIDE || result == IS_OUTSIDE) result = IS_OUTSIDE;
321 else if (hres == DOES_INTERSECT || result == DOES_INTERSECT) result = DOES_INTERSECT;
322 else result = IS_INSIDE;
323 }
324 return result;
325 }
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345 void Brick ::
GetPrimitiveData(const char * & classname,ARRAY<double> & coeffs) const346 GetPrimitiveData (const char *& classname, ARRAY<double> & coeffs) const
347 {
348 classname = "brick";
349 coeffs.SetSize(12);
350 coeffs.Elem(1) = p1(0);
351 coeffs.Elem(2) = p1(1);
352 coeffs.Elem(3) = p1(2);
353
354 coeffs.Elem(4) = p2(0);
355 coeffs.Elem(5) = p2(1);
356 coeffs.Elem(6) = p2(2);
357
358 coeffs.Elem(7) = p3(0);
359 coeffs.Elem(8) = p3(1);
360 coeffs.Elem(9) = p3(2);
361
362 coeffs.Elem(10) = p4(0);
363 coeffs.Elem(11) = p4(1);
364 coeffs.Elem(12) = p4(2);
365 }
366
SetPrimitiveData(ARRAY<double> & coeffs)367 void Brick :: SetPrimitiveData (ARRAY<double> & coeffs)
368 {
369 p1(0) = coeffs.Elem(1);
370 p1(1) = coeffs.Elem(2);
371 p1(2) = coeffs.Elem(3);
372
373 p2(0) = coeffs.Elem(4);
374 p2(1) = coeffs.Elem(5);
375 p2(2) = coeffs.Elem(6);
376
377 p3(0) = coeffs.Elem(7);
378 p3(1) = coeffs.Elem(8);
379 p3(2) = coeffs.Elem(9);
380
381 p4(0) = coeffs.Elem(10);
382 p4(1) = coeffs.Elem(11);
383 p4(2) = coeffs.Elem(12);
384
385 CalcData();
386 }
387
388
389
CalcData()390 void Brick :: CalcData()
391 {
392 v12 = p2 - p1;
393 v13 = p3 - p1;
394 v14 = p4 - p1;
395
396 Point<3> pi[8];
397 int i1, i2, i3;
398 int i, j;
399
400 i = 0;
401 for (i3 = 0; i3 <= 1; i3++)
402 for (i2 = 0; i2 <= 1; i2++)
403 for (i1 = 0; i1 <= 1; i1++)
404 {
405 pi[i] = p1 + i1 * v12 + i2 * v13 + i3 * v14;
406 i++;
407 }
408
409 static int lface[6][4] =
410 { { 1, 3, 2, 4 },
411 { 5, 6, 7, 8 },
412 { 1, 2, 5, 6 },
413 { 3, 7, 4, 8 },
414 { 1, 5, 3, 7 },
415 { 2, 4, 6, 8 } };
416
417 ARRAY<double> data(6);
418 for (i = 0; i < 6; i++)
419 {
420 const Point<3> lp1 = pi[lface[i][0]-1];
421 const Point<3> lp2 = pi[lface[i][1]-1];
422 const Point<3> lp3 = pi[lface[i][2]-1];
423
424 Vec<3> n = Cross ((lp2-lp1), (lp3-lp1));
425 n.Normalize();
426
427 for (j = 0; j < 3; j++)
428 {
429 data[j] = lp1(j);
430 data[j+3] = n(j);
431 }
432 faces[i] -> SetPrimitiveData (data);
433 /*
434 {
435 faces.Elem(i+1) -> SetPoints
436 (pi[lface[i][0]-1],
437 pi[lface[i][1]-1],
438 pi[lface[i][2]-1]);
439 }
440 */
441 }
442 }
443
444
Reduce(const BoxSphere<3> & box)445 void Brick :: Reduce (const BoxSphere<3> & box)
446 {
447 double val;
448 // Point<3> p;
449 Point<3> p[8];
450 for(int j=0;j<8;j++)
451 p[j]=box.GetPointNr(j);
452
453 for (int i = 0; i < 6; i++)
454 {
455 bool hasout = 0;
456 bool hasin = 0;
457 for (int j = 0; j < 8; j++)
458 {
459 // p = box.GetPointNr (j);
460 val = faces[i]->Plane::CalcFunctionValue (p[j]);
461 if (val > 0) hasout = 1;
462 else if (val < 0) hasin = 1;
463 if (hasout && hasin) break;
464 }
465 surfaceactive[i] = hasout && hasin;
466 }
467 }
468
UnReduce()469 void Brick :: UnReduce ()
470 {
471 for (int i = 0; i < 6; i++)
472 surfaceactive[i] = 1;
473 }
474
475
476
OrthoBrick(const Point<3> & ap1,const Point<3> & ap2)477 OrthoBrick :: OrthoBrick (const Point<3> & ap1, const Point<3> & ap2)
478 : Brick (ap1,
479 Point<3> (ap2(0), ap1(1), ap1(2)),
480 Point<3> (ap1(0), ap2(1), ap1(2)),
481 Point<3> (ap1(0), ap1(1), ap2(2)))
482 {
483 pmin = ap1;
484 pmax = ap2;
485 }
486
BoxInSolid(const BoxSphere<3> & box) const487 INSOLID_TYPE OrthoBrick :: BoxInSolid (const BoxSphere<3> & box) const
488 {
489 if (pmin(0) > box.PMax()(0) ||
490 pmin(1) > box.PMax()(1) ||
491 pmin(2) > box.PMax()(2) ||
492 pmax(0) < box.PMin()(0) ||
493 pmax(1) < box.PMin()(1) ||
494 pmax(2) < box.PMin()(2))
495 return IS_OUTSIDE;
496
497 if (pmin(0) < box.PMin()(0) &&
498 pmin(1) < box.PMin()(1) &&
499 pmin(2) < box.PMin()(2) &&
500 pmax(0) > box.PMax()(0) &&
501 pmax(1) > box.PMax()(1) &&
502 pmax(2) > box.PMax()(2))
503 return IS_INSIDE;
504
505 return DOES_INTERSECT;
506 }
507
508
Reduce(const BoxSphere<3> & box)509 void OrthoBrick :: Reduce (const BoxSphere<3> & box)
510 {
511 surfaceactive.Elem(1) =
512 (box.PMin()(2) < pmin(2)) && (pmin(2) < box.PMax()(2));
513 surfaceactive.Elem(2) =
514 (box.PMin()(2) < pmax(2)) && (pmax(2) < box.PMax()(2));
515
516 surfaceactive.Elem(3) =
517 (box.PMin()(1) < pmin(1)) && (pmin(1) < box.PMax()(1));
518 surfaceactive.Elem(4) =
519 (box.PMin()(1) < pmax(1)) && (pmax(1) < box.PMax()(1));
520
521 surfaceactive.Elem(5) =
522 (box.PMin()(0) < pmin(0)) && (pmin(0) < box.PMax()(0));
523 surfaceactive.Elem(6) =
524 (box.PMin()(0) < pmax(0)) && (pmax(0) < box.PMax()(0));
525 }
526 }
527