MIDAS  0.75
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
geometrylib.cpp
Go to the documentation of this file.
1 #include "geometrylib.h"
2 
3 #include "mathlib.h"
4 #include "arrays.h"
5 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 
10 
11 namespace midaspace {
12 
16 
18 // Function computes edge points of common perpendicular of two line segments l1 and l2.
19 // The line segment l1 goes from point l1a to l1b.
20 // l2 (l2a l2b)
21 //
22 // Return 0 - both edge points of common perpendicular is lying inside of the line segments l1 and l2
23 // Return 1 - edge point of common perpendicular is lying outside of the line segment l1
24 // Return 2 - edge point of common perpendicular is lying outside of the line segment l2
25 // Return 3 - both edge points of common perpendicular is lying outside of the line segments l1 and l2
26 //*/
27 //int give_edge_points_commperpend_skewlines (const double *l1a, const double *l1b, const double *l2a, const double *l2b, double *p1, double *p2)
28 //{
29 // Dvctr x(3), y(3), z(3);
30 //
31 // x.be_point2point( l1a, l1b );
32 // y.be_point2point( l2a, l2b );
33 // z.be_vectproduct( x, y );
34 //
35 // double m[9], r[3], l[3];
36 //
37 // m[0] = -x(0); m[1] = y(0); m[2] = z(0); r[0] = l1a[0] - l2a[0];
38 // m[3] = -x(1); m[4] = y(1); m[5] = z(1); r[1] = l1a[1] - l2a[1];
39 // m[6] = -x(2); m[7] = y(2); m[8] = z(2); r[2] = l1a[2] - l2a[2];
40 //
41 // solve_3 (m, r, l);
42 //
43 // p1[0] = l1a[0] + l[0] * x(0); p2[0] = l2a[0] + l[1] * y(0);
44 // p1[1] = l1a[1] + l[0] * x(1); p2[1] = l2a[1] + l[1] * y(1);
45 // p1[2] = l1a[2] + l[0] * x(2); p2[2] = l2a[2] + l[1] * y(2);
46 //
47 // int answer = 0;
48 // if ( l[0] < 0.0 || 1.0 < l[0] ) answer += 1;
49 // if ( l[1] < 0.0 || 1.0 < l[1] ) answer += 2;
50 //
51 // return answer;
52 //}
53 
54 
70 long intersec_rectangle3d_line (double zero, double norm,
71  const PoinT *A, const PoinT *B, const PoinT *C, const PoinT *D,
72  const PoinT *U, const PoinT *V,
73  double *ksi, double *eta, double *t,
74  PoinT *I1, PoinT *I2)
75 {
76  long i,ii,jj,kk,nr;
77 
78  VectoR a,b,c,d,e;
79  double aa[2],bb[2],cc[2],dd[2];
80 
81  if (!norm) {
82  a.beP2P(C,A);
83  b.beP2P(D,B);
84  norm = 0.5 * ( sqrt(a.x*a.x+a.y*a.y+a.z*a.z) + sqrt(b.x*b.x+b.y*b.y+b.z*b.z) );
85  }
86 
87  a.beP2P(B,A); a.add(b.beP2P(D,C));
88  b.beP2P(A,B);
89  c.beP2P(A,D);
90  d.beP2P(V,U);
91  e.beP2P(U,A);
92 
93  if (fabs(d.x) > fabs(d.y)) if (fabs(d.x) > fabs(d.z)) { ii=1; jj=2; kk=0; }
94  else { ii=0; jj=1; kk=2; }
95  else if (fabs(d.y) > fabs(d.z)) { ii=2; jj=0; kk=1; }
96  else { ii=0; jj=1; kk=2; }
97 
98  if (fabs(d[kk]) < norm*zero) _errorr("Length of line is 0.0");
99 
100  aa[0] = a[ii]*d[kk] - a[kk]*d[ii]; aa[1] = a[jj]*d[kk] - a[kk]*d[jj];
101  bb[0] = b[ii]*d[kk] - b[kk]*d[ii]; bb[1] = b[jj]*d[kk] - b[kk]*d[jj];
102  cc[0] = c[ii]*d[kk] - c[kk]*d[ii]; cc[1] = c[jj]*d[kk] - c[kk]*d[jj];
103  dd[0] = e[ii]*d[kk] - e[kk]*d[ii]; dd[1] = e[jj]*d[kk] - e[kk]*d[jj];
104 
105  nr = solv_2nle (norm*norm*zero*zero,aa,bb,cc,dd,ksi,eta);
106 
107  // nefunguje solv_2nle resic
108  if (nr == -1) errol;
109  if (nr == -2) {
110  if (0 < intersect_RayTriangle (zero, U, V, A, B, C, I1)) errol;
111  if (0 < intersect_RayTriangle (zero, U, V, A, C, D, I1)) errol;
112  if (0 < intersect_RayTriangle (zero, U, V, A, B, D, I1)) errol;
113  if (0 < intersect_RayTriangle (zero, U, V, B, C, D, I1)) errol;
114  nr = 0;
115  }
116 
117  for (i=0; i<nr; i++) {
118  t[i] = (a[kk]*ksi[i]*eta[i] + b[kk]*ksi[i] + c[kk]*eta[i] + e[kk]) / -d[kk] ;
119  t[i] = 2.0*( t[i]-0.5);
120  ksi[i] = -2.0*(ksi[i]-0.5);
121  eta[i] = -2.0*(eta[i]-0.5);
122  }
123 
124  if (nr > 0) I1->bePointAtAbscissa (U,V,t[0]);
125  if (nr > 1) I2->bePointAtAbscissa (U,V,t[1]);
126 
127  if (0){
128  fprintf (stdout,"\n number of results %ld\n",nr);
129  if (nr > 0) fprintf (stdout,"\n ksi = %5.2f eta = %5.2f t = %5.2f x = %15.7f y = %15.7f z = %15.7f\n", ksi[0],eta[0],t[0],I1->x,I1->y,I1->z);
130  if (nr > 1) fprintf (stdout,"\n ksi = %5.2f eta = %5.2f t = %5.2f x = %15.7f y = %15.7f z = %15.7f\n", ksi[1],eta[1],t[1],I2->x,I2->y,I2->z);
131  }
132 
133  return nr;
134 }
135 
136 
148 long nc_brick_3d (double zero, double x[], double y[], double z[], const PoinT *point, PoinT *answer)
149 {
150  long i,nite;
151  double a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8,c1,c2,c3,c4,c5,c6,c7,c8;
152  double xp,yp,zp,u,v,w;
153  double r[3],p[9],delta[3];
154 
155  xp = point->x; // ksi
156  yp = point->y; // 7
157  zp = point->z; // /
158  // /
159  a1 = x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]; // 8----------------7----------------7
160  a2 = -x[0] - x[1] + x[2] + x[3] - x[4] - x[5] + x[6] + x[7]; // /| * /|
161  a3 = -x[0] + x[1] + x[2] - x[3] - x[4] + x[5] + x[6] - x[7]; // / | * / |
162  a4 = x[0] + x[1] + x[2] + x[3] - x[4] - x[5] - x[6] - x[7]; // / | * / |
163  a5 = x[0] - x[1] + x[2] - x[3] + x[4] - x[5] + x[6] - x[7]; // / | * / |
164  a6 = -x[0] - x[1] + x[2] + x[3] + x[4] + x[5] - x[6] - x[7]; // / | * / |
165  a7 = -x[0] + x[1] + x[2] - x[3] + x[4] - x[5] - x[6] + x[7]; // 8 | 2 * 6 |
166  a8 = x[0] - x[1] + x[2] - x[3] - x[4] + x[5] - x[6] + x[7]; // / | * / |
167  // / | * / |
168  b1 = y[0] + y[1] + y[2] + y[3] + y[4] + y[5] + y[6] + y[7]; // / 12 5 / 11
169  b2 = -y[0] - y[1] + y[2] + y[3] - y[4] - y[5] + y[6] + y[7]; // / | * / |
170  b3 = -y[0] + y[1] + y[2] - y[3] - y[4] + y[5] + y[6] - y[7]; // / | * / |
171  b4 = y[0] + y[1] + y[2] + y[3] - y[4] - y[5] - y[6] - y[7]; // 5----------------5----------------6 |
172  b5 = y[0] - y[1] + y[2] - y[3] + y[4] - y[5] + y[6] - y[7]; // | | * | |
173  b6 = -y[0] - y[1] + y[2] + y[3] + y[4] + y[5] - y[6] - y[7]; // | | * | | eta
174  b7 = -y[0] + y[1] + y[2] - y[3] + y[4] - y[5] - y[6] + y[7]; // | 6 | O * * * * *|* * *4----------->
175  b8 = y[0] - y[1] + y[2] - y[3] - y[4] + y[5] - y[6] + y[7]; // | | * | |
176  // | | * | |
177  c1 = z[0] + z[1] + z[2] + z[3] + z[4] + z[5] + z[6] + z[7]; // | 4----------*-----3----|-----------3
178  c2 = -z[0] - z[1] + z[2] + z[3] - z[4] - z[5] + z[6] + z[7]; // | / * | /
179  c3 = -z[0] + z[1] + z[2] - z[3] - z[4] + z[5] + z[6] - z[7]; // | / * | /
180  c4 = z[0] + z[1] + z[2] + z[3] - z[4] - z[5] - z[6] - z[7]; // 9 / 3 * 10 /
181  c5 = z[0] - z[1] + z[2] - z[3] + z[4] - z[5] + z[6] - z[7]; // | / * | /
182  c6 = -z[0] - z[1] + z[2] + z[3] + z[4] + z[5] - z[6] - z[7]; // | / * | /
183  c7 = -z[0] + z[1] + z[2] - z[3] + z[4] - z[5] - z[6] + z[7]; // | 4 1 | 2
184  c8 = z[0] - z[1] + z[2] - z[3] - z[4] + z[5] - z[6] + z[7]; // | / * | /
185  // | / * | /
186  // setup initial guess // | / * | /
187  answer->zero(); // | / * | /
188  // |/ * |/
189  // apply Newton-Raphson to solve the problem // 1----------------1----------------2
190  nite = 0; // |
191  do { // |
192  if ((++nite) > 10) { // |
193  fprintf (stderr,"nc_hexa_3d: no convergence after 10 iterations"); // V dzeta
194  return 0;
195  }
196  u = answer->x;
197  v = answer->y;
198  w = answer->z;
199 
200  // compute the residual
201  r[0] = a1 + u * a2 + v * a3 + w * a4 + u * v * a5 + u * w * a6 + v * w * a7 + u * v * w * a8 - 8.0 * xp;
202  r[1] = b1 + u * b2 + v * b3 + w * b4 + u * v * b5 + u * w * b6 + v * w * b7 + u * v * w * b8 - 8.0 * yp;
203  r[2] = c1 + u * c2 + v * c3 + w * c4 + u * v * c5 + u * w * c6 + v * w * c7 + u * v * w * c8 - 8.0 * zp;
204 
205  // check for convergence
206  if ((r[0]*r[0]+r[1]*r[1]+r[2]*r[2]) < 1.e-20) break; // sqrt(1.e-20) = 1.e-10
207 
208  p[0] = a2 + v * a5 + w * a6 + v * w * a8; // *** JE TO ODVOZENO PRO TYTO BAZOVE FCE ***
209  p[1] = a3 + u * a5 + w * a7 + u * w * a8; // N[1] = 0.125 * (1.0 - ksi) * (1.0 - eta) * (1.0 + dzeta);
210  p[2] = a4 + u * a6 + v * a7 + u * v * a8; // N[2] = 0.125 * (1.0 - ksi) * (1.0 + eta) * (1.0 + dzeta);
211  p[3] = b2 + v * b5 + w * b6 + v * w * b8; // N[3] = 0.125 * (1.0 + ksi) * (1.0 + eta) * (1.0 + dzeta);
212  p[4] = b3 + u * b5 + w * b7 + u * w * b8; // N[4] = 0.125 * (1.0 + ksi) * (1.0 - eta) * (1.0 + dzeta);
213  p[5] = b4 + u * b6 + v * b7 + u * v * b8; // N[5] = 0.125 * (1.0 - ksi) * (1.0 - eta) * (1.0 - dzeta);
214  p[6] = c2 + v * c5 + w * c6 + v * w * c8; // N[6] = 0.125 * (1.0 - ksi) * (1.0 + eta) * (1.0 - dzeta);
215  p[7] = c3 + u * c5 + w * c7 + u * w * c8; // N[7] = 0.125 * (1.0 + ksi) * (1.0 + eta) * (1.0 - dzeta);
216  p[8] = c4 + u * c6 + v * c7 + u * v * c8; // N[8] = 0.125 * (1.0 + ksi) * (1.0 - eta) * (1.0 - dzeta);
217 
218  // solve for corrections
219  solve_3 (p,r,delta);
220 
221  // update guess
222  answer->x -= delta[0];
223  answer->y -= delta[1];
224  answer->z -= delta[2];
225 
226  } while (1);
227 
228  // test if inside
229  for (i=0;i<3;i++)
230  if ((answer[0][i] < -(1.0+zero)) || (answer[0][i] > (1.0+zero))) return 0;
231 
232  return 1;
233 }
234 
235 
236 // long nc_brick_3d (double zero,const element *elem,const double *point,double *answer)
237 // {
238 // long i,nite;
239 // double x[8],y[8],z[8];
240 // double a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8,c1,c2,c3,c4,c5,c6,c7,c8;
241 // double xp,yp,zp,u,v,w;
242 // double r[3],p[9],delta[3];
243 //
244 // for (i=0;i<8;i++){
245 // x[i] = Nodes[elem->give_node(i]]->coord[0]; // *****
246 // y[i] = Nodes[elem->give_node(i]]->coord[1]; // *(1)*
247 // z[i] = Nodes[elem->give_node(i]]->coord[2]; // *****
248 // }
249 //
250 // xp = point[0];
251 // yp = point[1];
252 // zp = point[2];
253 //
254 // a1 = x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]; // 8----------------7----------------7
255 // a2 = x[0] - x[1] - x[2] + x[3] + x[4] - x[5] - x[6] + x[7]; // /| /|
256 // a3 = x[0] + x[1] - x[2] - x[3] + x[4] + x[5] - x[6] - x[7]; // / | / |
257 // a4 = x[0] + x[1] + x[2] + x[3] - x[4] - x[5] - x[6] - x[7]; // / | / |
258 // a5 = x[0] - x[1] + x[2] - x[3] + x[4] - x[5] + x[6] - x[7]; // / | / |
259 // a6 = x[0] - x[1] - x[2] + x[3] - x[4] + x[5] + x[6] - x[7]; // / | / |
260 // a7 = x[0] + x[1] - x[2] - x[3] - x[4] - x[5] + x[6] + x[7]; // 8 | 2 6 |
261 // a8 = x[0] - x[1] + x[2] - x[3] - x[4] + x[5] - x[6] + x[7]; // / | / |
262 // // / | / |
263 // b1 = y[0] + y[1] + y[2] + y[3] + y[4] + y[5] + y[6] + y[7]; // / 12 5 / 11
264 // b2 = y[0] - y[1] - y[2] + y[3] + y[4] - y[5] - y[6] + y[7]; // / | / |
265 // b3 = y[0] + y[1] - y[2] - y[3] + y[4] + y[5] - y[6] - y[7]; // / | / |
266 // b4 = y[0] + y[1] + y[2] + y[3] - y[4] - y[5] - y[6] - y[7]; // 5----------------5----------------6 |
267 // b5 = y[0] - y[1] + y[2] - y[3] + y[4] - y[5] + y[6] - y[7]; // | | | |
268 // b6 = y[0] - y[1] - y[2] + y[3] - y[4] + y[5] + y[6] - y[7]; // ksi | | | |
269 // b7 = y[0] + y[1] - y[2] - y[3] - y[4] - y[5] + y[6] + y[7]; // <---------|* * *6* * * * * * * * O | 4 |
270 // b8 = y[0] - y[1] + y[2] - y[3] - y[4] + y[5] - y[6] + y[7]; // | | ** | |
271 // // | | * * | |
272 // c1 = z[0] + z[1] + z[2] + z[3] + z[4] + z[5] + z[6] + z[7]; // | 4-------*--*-----3----|-----------3
273 // c2 = z[0] - z[1] - z[2] + z[3] + z[4] - z[5] - z[6] + z[7]; // | / * * | /
274 // c3 = z[0] + z[1] - z[2] - z[3] + z[4] + z[5] - z[6] - z[7]; // | / * * | /
275 // c4 = z[0] + z[1] + z[2] + z[3] - z[4] - z[5] - z[6] - z[7]; // 9 / 3 * 10 /
276 // c5 = z[0] - z[1] + z[2] - z[3] + z[4] - z[5] + z[6] - z[7]; // | / / * | /
277 // c6 = z[0] - z[1] - z[2] + z[3] - z[4] + z[5] + z[6] - z[7]; // | / / * | /
278 // c7 = z[0] + z[1] - z[2] - z[3] - z[4] - z[5] + z[6] + z[7]; // | 4 / 1 | 2
279 // c8 = z[0] - z[1] + z[2] - z[3] - z[4] + z[5] - z[6] + z[7]; // | / / * | /
280 // // | / / * | /
281 // // setup initial guess // | / / * | /
282 // answer[0] = answer[1] = answer[2] = 0.0; // | / / * | /
283 // // |/ / * |/
284 // // apply Newton-Raphson to solve the problem // 1-------/--------1----------------2
285 // nite = 0; // / |
286 // do { // / |
287 // if ((++nite) > 10) { // eta / |
288 // fprintf (stderr,"nc_hexa_3d: no convergence after 10 iterations"); // L V dzeta
289 // return 0;
290 // }
291 // u = answer[0];
292 // v = answer[1];
293 // w = answer[2];
294 //
295 // // compute the residual
296 // r[0] = a1 + u * a2 + v * a3 + w * a4 + u * v * a5 + u * w * a6 + v * w * a7 + u * v * w * a8 - 8.0 * xp;
297 // r[1] = b1 + u * b2 + v * b3 + w * b4 + u * v * b5 + u * w * b6 + v * w * b7 + u * v * w * b8 - 8.0 * yp;
298 // r[2] = c1 + u * c2 + v * c3 + w * c4 + u * v * c5 + u * w * c6 + v * w * c7 + u * v * w * c8 - 8.0 * zp;
299 //
300 // // check for convergence
301 // if ((r[0]*r[0]+r[1]*r[1]+r[2]*r[2]) < 1.e-20) break; // sqrt(1.e-20) = 1.e-10
302 //
303 // p[0] = a2 + v * a5 + w * a6 + v * w * a8; // *** JE TO ODVOZENO PRO TYTO BAZOVE FCE ***
304 // p[1] = a3 + u * a5 + w * a7 + u * w * a8; // N[1] = 0.125 * (1.0 + ksi) * (1.0 + eta) * (1.0 + dzeta);
305 // p[2] = a4 + u * a6 + v * a7 + u * v * a8; // N[2] = 0.125 * (1.0 - ksi) * (1.0 + eta) * (1.0 + dzeta);
306 // p[3] = b2 + v * b5 + w * b6 + v * w * b8; // N[3] = 0.125 * (1.0 - ksi) * (1.0 - eta) * (1.0 + dzeta);
307 // p[4] = b3 + u * b5 + w * b7 + u * w * b8; // N[4] = 0.125 * (1.0 + ksi) * (1.0 - eta) * (1.0 + dzeta);
308 // p[5] = b4 + u * b6 + v * b7 + u * v * b8; // N[5] = 0.125 * (1.0 + ksi) * (1.0 + eta) * (1.0 - dzeta);
309 // p[6] = c2 + v * c5 + w * c6 + v * w * c8; // N[6] = 0.125 * (1.0 - ksi) * (1.0 + eta) * (1.0 - dzeta);
310 // p[7] = c3 + u * c5 + w * c7 + u * w * c8; // N[7] = 0.125 * (1.0 - ksi) * (1.0 - eta) * (1.0 - dzeta);
311 // p[8] = c4 + u * c6 + v * c7 + u * v * c8; // N[8] = 0.125 * (1.0 + ksi) * (1.0 - eta) * (1.0 - dzeta);
312 //
313 // // solve for corrections
314 // solve_3 (p,r,delta);
315 //
316 // // update guess
317 // answer[0] -= delta[0];
318 // answer[1] -= delta[1];
319 // answer[2] -= delta[2];
320 //
321 // } while (1);
322 //
323 // // test if inside
324 // for (i=0;i<3;i++)
325 // if ((answer[i] < -(1.0+zero)) || (answer[i] > (1.0+zero))) return 0;
326 //
327 // return 1;
328 // }
329 
338 
339 // a Triangle is given by three points: Point V0, V1, V2
340 // a Polygon is given by:
341 // int n = number of vertex points
342 // Point* V[] = an array of points with V[n]=V[0], V[n+1]=V[1]
343 
344 // Note: for efficiency low-level functions are declared to be inline.
345 
346 
347 // // isLeft(): tests if a point is Left|On|Right of an infinite line.
348 // // Input: three points P0, P1, and P2
349 // // Return: >0 for P2 left of the line through P0 and P1
350 // // =0 for P2 on the line
351 // // <0 for P2 right of the line
352 // inline int
353 // isLeft( Point P0, Point P1, Point P2 )
354 // {
355 // return ( (P1.x - P0.x) * (P2.y - P0.y)
356 // - (P2.x - P0.x) * (P1.y - P0.y) );
357 // }
358 // //===================================================================
359 //
360 // // orientation2D_Triangle(): test the orientation of a triangle
361 // // Input: three vertex points V0, V1, V2
362 // // Return: >0 for counterclockwise
363 // // =0 for none (degenerate)
364 // // <0 for clockwise
365 // inline int
366 // orientation2D_Triangle( Point V0, Point V1, Point V2 )
367 // {
368 // return isLeft(V0, V1, V2);
369 // }
370 // //===================================================================
371 //
372 // // area2D_Triangle(): compute the area of a triangle
373 // // Input: three vertex points V0, V1, V2
374 // // Return: the (float) area of T
375 // inline float
376 // area2D_Triangle( Point V0, Point V1, Point V2 )
377 // {
378 // return (float)isLeft(V0, V1, V2) / 2.0;
379 // }
380 // //===================================================================
381 //
382 // // orientation2D_Polygon(): tests the orientation of a simple polygon
383 // // Input: int n = the number of vertices in the polygon
384 // // Point* V = an array of n+1 vertices with V[n]=V[0]
385 // // Return: >0 for counterclockwise
386 // // =0 for none (degenerate)
387 // // <0 for clockwise
388 // // Note: this algorithm is faster than computing the signed area.
389 // int
390 // orientation2D_Polygon( int n, Point* V )
391 // {
392 // // first find rightmost lowest vertex of the polygon
393 // int rmin = 0;
394 // int xmin = V[0].x;
395 // int ymin = V[0].y;
396 //
397 // for (int i=1; i<n; i++) {
398 // if (V[i].y > ymin)
399 // continue;
400 // if (V[i].y == ymin) { // just as low
401 // if (V[i].x < xmin) // and to left
402 // continue;
403 // }
404 // rmin = i; // a new rightmost lowest vertex
405 // xmin = V[i].x;
406 // ymin = V[i].y;
407 // }
408 //
409 // // test orientation at this rmin vertex
410 // // ccw <=> the edge leaving is left of the entering edge
411 // if (rmin == 0)
412 // return isLeft( V[n-1], V[0], V[1] );
413 // else
414 // return isLeft( V[rmin-1], V[rmin], V[rmin+1] );
415 // }
416 // //===================================================================
417 //
418 // // area2D_Polygon(): computes the area of a 2D polygon
419 // // Input: int n = the number of vertices in the polygon
420 // // Point* V = an array of n+2 vertices
421 // // with V[n]=V[0] and V[n+1]=V[1]
422 // // Return: the (float) area of the polygon
423 // float
424 // area2D_Polygon( int n, Point* V )
425 // {
426 // float area = 0;
427 // int i, j, k; // indices
428 //
429 // for (i=1, j=2, k=0; i<=n; i++, j++, k++) {
430 // area += V[i].x * (V[j].y - V[k].y);
431 // }
432 // return area / 2.0;
433 // }
434 // //===================================================================
435 
436 // area3D_Polygon(): computes the area of a 3D planar polygon
437 // Input: int n = the number of vertices in the polygon
438 // Point* V = an array of n+2 vertices in a plane
439 // with V[n]=V[0] and V[n+1]=V[1]
440 // Point N = unit normal vector of the polygon's plane
441 // Return: the (float) area of the polygon
442  //
443 double area3D_Polygon (int n, const PoinT **V, const VectoR *N)
444 {
445  double area = 0;
446  double an, ax, ay, az; // abs value of normal and its coords
447  int coord; // coord to ignore: 1=x, 2=y, 3=z
448  int i, j, k; // loop indices
449 
450  // select largest abs coordinate to ignore for projection
451  ax = (N->x>0 ? N->x : -N->x); // abs x-coord
452  ay = (N->y>0 ? N->y : -N->y); // abs y-coord
453  az = (N->z>0 ? N->z : -N->z); // abs z-coord
454 
455  coord = 3; // ignore z-coord
456  if (ax > ay) { if (ax > az) coord = 1; } // ignore x-coord
457  else { if (ay > az) coord = 2; } // ignore y-coord
458 
459  // compute area of the 2D projection
460  for (i=1, j=2, k=0; i<=n; i++, j++, k++)
461  switch (coord) {
462  case 1: area += (V[i]->y * (V[j]->z - V[k]->z)); continue;
463  case 2: area += (V[i]->x * (V[j]->z - V[k]->z)); continue;
464  case 3: area += (V[i]->x * (V[j]->y - V[k]->y)); continue;
465  }
466 
467  // scale to get area before projection
468  an = sqrt( ax*ax + ay*ay + az*az); // length of normal vector
469  switch (coord) {
470  case 1: area *= (an / (2*ax)); break;
471  case 2: area *= (an / (2*ay)); break;
472  case 3: area *= (an / (2*az));
473  }
474 
475  return area;
476 }
477 
478 //===================================================================
479 // http://www.softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm
480 // Line and Ray and Segment with defining points {Point P0, P1;}
481 // (a Line is infinite, Rays and Segments start at P0)
482 // (a Ray extends beyond P1, but a Segment ends at P1)
483 // Plane with a point and a normal {Point V0; Vector n;}
484 // Triangle with defining vertices {Point V0, V1, V2;}
485 // Polyline and Polygon with n vertices {int n; Point *V;}
486 // (a Polygon has V[n]=V[0])
487 //===================================================================
488 // intersect_RayTriangle(): intersect a ray with a 3D triangle
489 // Input: a ray R, and a triangle T
490 // Output: *I = intersection point (when it exists)
491 // Return: -1 = triangle is degenerate (a segment or point)
492 // 0 = disjoint (no intersect)
493 // 1 = intersect in unique point I1
494 // 2 = are in the same plane
495 int intersect_RayTriangle (double zero, const PoinT *P0, const PoinT *P1, const PoinT *V0, const PoinT *V1, const PoinT *V2, PoinT* I)
496 {
497  VectoR u, v, n; // triangle vectors
498  VectoR dir, w0, w; // ray vectors
499  double r, a, b; // params to calc ray-plane intersect
500 
501  // get triangle edge vectors and plane normal
502  u.beP2P (V0, V1);
503  v.beP2P (V0, V2);
504  n.beVectProduct (&u, &v); // cross product
505 
506  // triangle is degenerate
507  if (u.give_length() < zero || v.give_length() < zero || n.give_length() < zero)
508  return -1;
509 
510  dir.beP2P (P0, P1); // ray direction vector
511  w0.beP2P (V0, P0);
512 
513  a = - n.giveScalProduct(&w0);
514  b = n.giveScalProduct(&dir);
515  // ray is parallel to triangle plane
516  if (fabs(b) < zero) {
517  if (fabs(a) < zero) return 2; // ray lies in triangle plane
518  else return 0; // ray disjoint from plane
519  }
520 
521  // get intersect point of ray with triangle plane
522  r = a / b;
523  if (r < 0.0 || r > 1.0) // ray goes away from triangle plane
524  return 0; // => no intersect
525  // for a ray, do not test if(r > 1.0)
526 
527 
528  I->copy(P0)->add(dir.tms(r)); // intersect point of ray and plane
529 
530  // is I inside T?
531  double uu, uv, vv, wu, wv, D;
532  uu = u.giveScalProduct(&u);
533  uv = u.giveScalProduct(&v);
534  vv = v.giveScalProduct(&v);
535 
536  w.beP2P (V0, I);
537  wu = w.giveScalProduct(&u);
538  wv = w.giveScalProduct(&v);
539 
540  D = uv * uv - uu * vv;
541 
542  // get and test parametric coords
543  double s, t;
544  s = (uv * wv - vv * wu) / D;
545  if (s < 0.0-zero || s > 1.0+zero) return 0; // I is outside T
546 
547  t = (uv * wu - uu * wv) / D;
548  if (t < 0.0-zero || (s + t) > 1.0+zero) return 0; // I is outside T
549 
550  return 1; // I is in T
551 }
552 
553 
554 
555 
556 
557 
558 //===================================================================
559 
560 
567 void dx_bf_lin_hex_3d (double *n,double y,double z)
568 {
569  n[0]=(1.0+y)*(1.0+z)/8.0;
570  n[1]=(1.0+y)*(1.0+z)/(-8.0);
571  n[2]=(1.0-y)*(1.0+z)/(-8.0);
572  n[3]=(1.0-y)*(1.0+z)/8.0;
573  n[4]=(1.0+y)*(1.0-z)/8.0;
574  n[5]=(1.0+y)*(1.0-z)/(-8.0);
575  n[6]=(1.0-y)*(1.0-z)/(-8.0);
576  n[7]=(1.0-y)*(1.0-z)/8.0;
577 }
578 
579 
580 
587 void dy_bf_lin_hex_3d (double *n,double x,double z)
588 {
589  n[0]=(1.0+x)*(1.0+z)/8.0;
590  n[1]=(1.0-x)*(1.0+z)/8.0;
591  n[2]=(1.0-x)*(1.0+z)/(-8.0);
592  n[3]=(1.0+x)*(1.0+z)/(-8.0);
593  n[4]=(1.0+x)*(1.0-z)/8.0;
594  n[5]=(1.0-x)*(1.0-z)/8.0;
595  n[6]=(1.0-x)*(1.0-z)/(-8.0);
596  n[7]=(1.0+x)*(1.0-z)/(-8.0);
597 }
598 
599 
600 
607 void dz_bf_lin_hex_3d (double *n,double x,double y)
608 {
609  n[0]=(1.0+x)*(1.0+y)/8.0;
610  n[1]=(1.0-x)*(1.0+y)/8.0;
611  n[2]=(1.0-x)*(1.0-y)/8.0;
612  n[3]=(1.0+x)*(1.0-y)/8.0;
613  n[4]=(1.0+x)*(1.0+y)/(-8.0);
614  n[5]=(1.0-x)*(1.0+y)/(-8.0);
615  n[6]=(1.0-x)*(1.0-y)/(-8.0);
616  n[7]=(1.0+x)*(1.0-y)/(-8.0);
617 }
618 
619 
620 void jac_3d (double &jac, Dvctr &x,Dvctr &y,Dvctr &z,
621  double xi,double eta,double zeta)
622 {
623  long i;
624  double d1,d2,d3,d4,d5,d6,d7,d8,d9,*nx,*ny,*nz;
625 
626  nx = new double [x.give_size()];
627  ny = new double [x.give_size()];
628  nz = new double [x.give_size()];
629 
630  switch (x.give_size()){
631  //case 4:{
632  // dx_bf_lin_tet (nx);
633  // dy_bf_lin_tet (ny);
634  // dz_bf_lin_tet (nz);
635  // break;
636  //}
637  //case 6:{
638  // dx_bf_lin_wed_3d (nx,zeta);
639  // dy_bf_lin_wed_3d (ny,zeta);
640  // dz_bf_lin_wed_3d (nz,xi,eta);
641  // break;
642  //}
643  case 8:{
644  dx_bf_lin_hex_3d (nx,eta,zeta);
645  dy_bf_lin_hex_3d (ny,xi,zeta);
646  dz_bf_lin_hex_3d (nz,xi,eta);
647  break;
648  }
649  //case 10:{
650  // dx_bf_quad_tet (nx,xi,eta,zeta);
651  // dy_bf_quad_tet (ny,xi,eta,zeta);
652  // dz_bf_quad_tet (nz,xi,eta,zeta);
653  // break;
654  //}
655  //case 15:{
656  // dx_bf_quad_wed_3d (nx,xi,eta,zeta);
657  // dy_bf_quad_wed_3d (ny,xi,eta,zeta);
658  // dz_bf_quad_wed_3d (nz,xi,eta,zeta);
659  // break;
660  //}
661  //case 20:{
662  // dx_bf_quad_hex_3d (nx,xi,eta,zeta);
663  // dy_bf_quad_hex_3d (ny,xi,eta,zeta);
664  // dz_bf_quad_hex_3d (nz,xi,eta,zeta);
665  // break;
666  //}
667  default:{
668  fprintf (stderr,"\n\n wrong number of nodes on 3D element in the function");
669  fprintf (stderr,"\n jac_3d (%s, line %d).\n",__FILE__,__LINE__);
670  }
671  }
672 
673  d1=0.0; d2=0.0; d3=0.0; d4=0.0;
674  d5=0.0; d6=0.0; d7=0.0; d8=0.0; d9=0.0;
675  for (i=0;i<x.give_size();i++) {
676  d1+=nx[i]*x[i];
677  d2+=ny[i]*x[i];
678  d3+=nz[i]*x[i];
679  d4+=nx[i]*y[i];
680  d5+=ny[i]*y[i];
681  d6+=nz[i]*y[i];
682  d7+=nx[i]*z[i];
683  d8+=ny[i]*z[i];
684  d9+=nz[i]*z[i];
685  }
686 
687  jac = d1*d5*d9 + d4*d8*d3 + d7*d2*d6 - d7*d5*d3 - d8*d6*d1 - d9*d4*d2;
688 
689  delete [] nx; delete [] ny; delete [] nz;
690 }
691 
692 
693 } // namespace midaspace
midaspace::dy_bf_lin_hex_3d
void dy_bf_lin_hex_3d(double *n, double x, double z)
function computes derivatives of tri-linear base functions on hexahedron with respect of natural coor...
Definition: geometrylib.cpp:587
midaspace::solv_2nle
long solv_2nle(double zero, const double *a, const double *b, const double *c, const double *d, double *x, double *y)
function solves system of two non-linear equations: a[0]*x*y + b[0]*x + c[0]*y + d[0] = 0 a[1]*x*y + ...
Definition: mathlib.cpp:169
midaspace::intersect_RayTriangle
int intersect_RayTriangle(double zero, const PoinT *P0, const PoinT *P1, const PoinT *V0, const PoinT *V1, const PoinT *V2, PoinT *I)
Definition: geometrylib.cpp:495
midaspace::Xvctr::give_size
virtual long give_size(void) const
Definition: arrays.h:387
midaspace::nc_brick_3d
long nc_brick_3d(double zero, double x[], double y[], double z[], const PoinT *point, PoinT *answer)
Function computes natural coordinates of 'point' on 'element', 'point' is entered in cartesian coordi...
Definition: geometrylib.cpp:148
geometrylib.h
Geometry functions.
midaspace::PoinT::bePointAtAbscissa
void bePointAtAbscissa(const PoinT *p1, const PoinT *p2, double ksi)
receiver will be point at abscissa p1p2 with natural coord ksi
Definition: arrays.h:107
errol
#define errol
Definition: gelib.h:142
midaspace::VectoR::give_length
double give_length(void) const
Definition: arrays.h:188
midaspace::Elem3D::zero
Elem3D * zero(void)
Definition: arrays.h:64
arrays.h
Structs Elem3D, PoinT and VectoR; classes Array, Array1d, Xscal, Dscal, Xvctr, Lvctr, Dvctr, Xmtrx, Lmtrx and Dmtrx.
midaspace::solve_3
void solve_3(const double *m, const double *r, double *l)
Function solves the system of linear equations.
Definition: mathlib.cpp:103
midaspace::intersec_rectangle3d_line
long intersec_rectangle3d_line(double zero, double norm, const PoinT *A, const PoinT *B, const PoinT *C, const PoinT *D, const PoinT *U, const PoinT *V, double *ksi, double *eta, double *t, PoinT *I1, PoinT *I2)
ZDROJE http://softsurfer.com/algorithm_archive.htm.
Definition: geometrylib.cpp:70
midaspace::Elem3D::tms
Elem3D * tms(double val)
Definition: arrays.h:59
midaspace::dx_bf_lin_hex_3d
void dx_bf_lin_hex_3d(double *n, double y, double z)
function computes derivatives of tri-linear base functions on hexahedron with respect of natural coor...
Definition: geometrylib.cpp:567
midaspace::VectoR::beVectProduct
VectoR * beVectProduct(const VectoR *v1, const VectoR *v2)
vector product v1 x v2 (cross product)
Definition: arrays.h:160
_errorr
#define _errorr(_1)
Definition: gelib.h:151
mathlib.h
Mathematic functions.
midaspace::Elem3D::add
Elem3D * add(const Elem3D *p)
Definition: arrays.h:61
midaspace::VectoR::beP2P
VectoR * beP2P(const PoinT *p1, const PoinT *p2)
Receiver is set to vector from p1 to p2, i.e. 'this = p2 - p1'.
Definition: arrays.h:153
midaspace::dz_bf_lin_hex_3d
void dz_bf_lin_hex_3d(double *n, double x, double y)
function computes derivatives of tri-linear base functions on hexahedron with respect of natural coor...
Definition: geometrylib.cpp:607
midaspace::Elem3D::giveScalProduct
double giveScalProduct(const Elem3D *v) const
scalar product this * e
Definition: arrays.h:80
midaspace::PoinT::copy
PoinT * copy(const PoinT *p)
Definition: arrays.h:99
midaspace::area3D_Polygon
double area3D_Polygon(int n, const PoinT **V, const VectoR *N)
************************** WWW.SOFTSURFER.COM ALGORITHMS ************************** Copyright 2000...
Definition: geometrylib.cpp:443
midaspace::jac_3d
void jac_3d(double &jac, Dvctr &x, Dvctr &y, Dvctr &z, double xi, double eta, double zeta)
Definition: geometrylib.cpp:620
This reference manual is part of the MIDAS documentation. Copyright (c) 2005-2013 Ladislav Termit Svoboda. Generated on Fri Mar 21 2014 11:39:42 for MIDAS by doxygen 1.8.6.