73 double *ksi,
double *eta,
double *t,
79 double aa[2],bb[2],cc[2],dd[2];
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) );
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; }
98 if (fabs(d[kk]) < norm*zero)
_errorr(
"Length of line is 0.0");
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];
105 nr =
solv_2nle (norm*norm*zero*zero,aa,bb,cc,dd,ksi,eta);
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);
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);
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];
159 a1 = x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[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];
166 a8 = x[0] - x[1] + x[2] - x[3] - x[4] + x[5] - x[6] + x[7];
168 b1 = y[0] + y[1] + y[2] + y[3] + y[4] + y[5] + y[6] + y[7];
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];
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];
174 b7 = -y[0] + y[1] + y[2] - y[3] + y[4] - y[5] - y[6] + y[7];
175 b8 = y[0] - y[1] + y[2] - y[3] - y[4] + y[5] - y[6] + y[7];
177 c1 = z[0] + z[1] + z[2] + z[3] + z[4] + z[5] + z[6] + z[7];
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];
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];
184 c8 = z[0] - z[1] + z[2] - z[3] - z[4] + z[5] - z[6] + z[7];
193 fprintf (stderr,
"nc_hexa_3d: no convergence after 10 iterations");
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;
206 if ((r[0]*r[0]+r[1]*r[1]+r[2]*r[2]) < 1.e-20)
break;
208 p[0] = a2 + v * a5 + w * a6 + v * w * a8;
209 p[1] = a3 + u * a5 + w * a7 + u * w * a8;
210 p[2] = a4 + u * a6 + v * a7 + u * v * a8;
211 p[3] = b2 + v * b5 + w * b6 + v * w * b8;
212 p[4] = b3 + u * b5 + w * b7 + u * w * b8;
213 p[5] = b4 + u * b6 + v * b7 + u * v * b8;
214 p[6] = c2 + v * c5 + w * c6 + v * w * c8;
215 p[7] = c3 + u * c5 + w * c7 + u * w * c8;
216 p[8] = c4 + u * c6 + v * c7 + u * v * c8;
222 answer->
x -= delta[0];
223 answer->
y -= delta[1];
224 answer->
z -= delta[2];
230 if ((answer[0][i] < -(1.0+zero)) || (answer[0][i] > (1.0+zero)))
return 0;
446 double an, ax, ay, az;
451 ax = (N->
x>0 ? N->
x : -N->
x);
452 ay = (N->
y>0 ? N->
y : -N->
y);
453 az = (N->
z>0 ? N->
z : -N->
z);
456 if (ax > ay) {
if (ax > az) coord = 1; }
457 else {
if (ay > az) coord = 2; }
460 for (i=1, j=2, k=0; i<=n; i++, j++, k++)
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;
468 an = sqrt( ax*ax + ay*ay + az*az);
470 case 1: area *= (an / (2*ax));
break;
471 case 2: area *= (an / (2*ay));
break;
472 case 3: area *= (an / (2*az));
516 if (fabs(b) < zero) {
517 if (fabs(a) < zero)
return 2;
523 if (r < 0.0 || r > 1.0)
531 double uu, uv, vv, wu, wv, D;
540 D = uv * uv - uu * vv;
544 s = (uv * wv - vv * wu) / D;
545 if (s < 0.0-zero || s > 1.0+zero)
return 0;
547 t = (uv * wu - uu * wv) / D;
548 if (t < 0.0-zero || (s + t) > 1.0+zero)
return 0;
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;
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);
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);
621 double xi,
double eta,
double zeta)
624 double d1,d2,d3,d4,d5,d6,d7,d8,d9,*nx,*ny,*nz;
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__);
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;
687 jac = d1*d5*d9 + d4*d8*d3 + d7*d2*d6 - d7*d5*d3 - d8*d6*d1 - d9*d4*d2;
689 delete [] nx;
delete [] ny;
delete [] nz;
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...
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 + ...
int intersect_RayTriangle(double zero, const PoinT *P0, const PoinT *P1, const PoinT *V0, const PoinT *V1, const PoinT *V2, PoinT *I)
virtual long give_size(void) const
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...
void bePointAtAbscissa(const PoinT *p1, const PoinT *p2, double ksi)
receiver will be point at abscissa p1p2 with natural coord ksi
double give_length(void) const
Structs Elem3D, PoinT and VectoR; classes Array, Array1d, Xscal, Dscal, Xvctr, Lvctr, Dvctr, Xmtrx, Lmtrx and Dmtrx.
void solve_3(const double *m, const double *r, double *l)
Function solves the system of linear equations.
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.
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...
VectoR * beVectProduct(const VectoR *v1, const VectoR *v2)
vector product v1 x v2 (cross product)
Elem3D * add(const Elem3D *p)
VectoR * beP2P(const PoinT *p1, const PoinT *p2)
Receiver is set to vector from p1 to p2, i.e. 'this = p2 - p1'.
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...
double giveScalProduct(const Elem3D *v) const
scalar product this * e
PoinT * copy(const PoinT *p)
double area3D_Polygon(int n, const PoinT **V, const VectoR *N)
************************** WWW.SOFTSURFER.COM ALGORITHMS ************************** Copyright 2000...
void jac_3d(double &jac, Dvctr &x, Dvctr &y, Dvctr &z, double xi, double eta, double zeta)