gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
libol1.c
Go to the documentation of this file.
1 // All libOL code is Copyright 2012-2018 - by Loïc Maréchal / INRIA.
2 // This program is a free software.
3 // You can redistribute it and/or modify it under the terms of the MIT License
4 // as published by the Open Source Initiative.
5 //
6 // This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
7 // without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
8 // See the MIT License for more details.
9 //
10 // You should have received a copy of the MIT License along with this program
11 // as the file LICENSE.txt; if not, please see:
12 // https://opensource.org/licenses/MIT
13 //
14 //
15 // MIT License
16 //
17 // Copyright (c) 2012-2018 Loïc Maréchal / INRIA
18 //
19 // Permission is hereby granted, free of charge, to any person obtaining a copy
20 // of this software and associated documentation files (the "Software"), to deal
21 // in the Software without restriction, including without limitation the rights
22 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
23 // copies of the Software, and to permit persons to whom the Software is
24 // furnished to do so, subject to the following conditions:
25 //
26 // The above copyright notice and this permission notice shall be included in all
27 // copies or substantial portions of the Software.
28 //
29 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
30 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
31 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
32 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
33 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
34 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
35 // SOFTWARE.
36 
37 
38 /*----------------------------------------------------------------------------*/
39 /* */
40 /* LIB OCTREE LOCALISATION V1.70 */
41 /* */
42 /*----------------------------------------------------------------------------*/
43 /* */
44 /* Description: Octree for mesh localization */
45 /* Author: Loic MARECHAL */
46 /* Creation date: mar 16 2012 */
47 /* Last modification: feb 03 2021 */
48 /* */
49 /*----------------------------------------------------------------------------*/
50 
51 
52 /*----------------------------------------------------------------------------*/
53 /* ANSI C headers */
54 /*----------------------------------------------------------------------------*/
55 
56 #include <assert.h>
57 #include <fcntl.h>
58 #include <limits.h>
59 #include <float.h>
60 #include <math.h>
61 #include <stdio.h>
62 #include <stdlib.h>
63 #include <string.h>
64 #include <stdint.h>
65 #include "libol1.h"
66 
67 
68 /*----------------------------------------------------------------------------*/
69 /* Defines and macros */
70 /*----------------------------------------------------------------------------*/
71 
72 #define MaxItmOct 20
73 #define MaxOctLvl 10
74 #define MinGrdLvl 5
75 #define ItmPerBuc 100
76 #define MemBlkSiz 100000
77 #define TngFlg 1
78 #define AniFlg 2
79 #define MaxThr 256
80 #define MIN(a,b) ((a) < (b) ? (a) : (b))
81 #define MAX(a,b) ((a) > (b) ? (a) : (b))
82 #define POW(a) ((a)*(a))
83 #define CUB(a) ((a)*(a)*(a))
84 
85 
86 /*----------------------------------------------------------------------------*/
87 /* Local structures */
88 /*----------------------------------------------------------------------------*/
89 
90 typedef struct
91 {
92  itg idx;
93  fpn crd[3];
94 }VerSct;
95 
96 typedef struct
97 {
98  VerSct *ver[2];
99  fpn tng[3], siz;
100  itg idx;
101 }EdgSct;
102 
103 typedef struct
104 {
105  EdgSct edg[3];
106  VerSct *ver[3];
107  fpn nrm[3];
108  itg idx;
109  float ani;
110 }TriSct;
111 
112 typedef struct
113 {
114  VerSct *ver[4];
115  EdgSct edg[4];
116  TriSct tri[2];
117  fpn nrm[3];
118  itg idx;
119 }QadSct;
120 
121 typedef struct
122 {
123  VerSct *ver[4];
124  EdgSct edg[6];
125  TriSct tri[4];
126  itg idx;
127  float ani;
128 }TetSct;
129 
130 typedef struct
131 {
132  VerSct *ver[5];
133  EdgSct edg[8];
134  TriSct tri[5][2];
135  itg idx;
136 }PyrSct;
137 
138 typedef struct
139 {
140  VerSct *ver[6];
141  EdgSct edg[9];
142  TriSct tri[5][2];
143  itg idx;
144 }PriSct;
145 
146 typedef struct
147 {
148  VerSct *ver[8];
149  EdgSct edg[12];
150  QadSct qad[6];
151  itg idx;
152 }HexSct;
153 
154 typedef struct LnkSctPtr
155 {
156  itg typ, idx;
157  struct LnkSctPtr *nex;
158 }LnkSct;
159 
160 typedef struct MemSctPtr
161 {
162  size_t siz;
163  void *adr;
164  struct MemSctPtr *nex;
165 }MemSct;
166 
167 typedef struct
168 {
169  VerSct ver[8], BakVer[8];
170  EdgSct edg, BakEdg;
171  TriSct tri, BakTri;
172  QadSct qad, BakQad;
173  TetSct tet, BakTet;
174  PyrSct pyr, BakPyr;
175  PriSct pri, BakPri;
176  HexSct hex, BakHex;
177  char *FlgTab;
178  itg *TagTab, tag;
179 }MshThrSct;
180 
181 typedef struct
182 {
183  MshThrSct thr[ MaxThr ];
184  size_t UsrSiz[ LolNmbTyp ], NmbItm[ LolNmbTyp ];
185  fpn aniso, eps;
186  itg BasIdx;
187  char *UsrPtr[ LolNmbTyp ];
188 }MshSct;
189 
190 typedef struct OctSctPtr
191 {
192  union
193  {
194  LnkSct *lnk;
195  struct OctSctPtr *son;
196  };
197  unsigned char NmbVer, NmbEdg, NmbFac, NmbVol, lvl, sub, ani, MaxItm;
198 }OctSct;
199 
200 typedef struct
201 {
202  OctSct *oct;
203  itg idx;
204  char pos[3];
205 }BucSct;
206 
207 typedef struct
208 {
209  VerSct ver[8];
210  HexSct hex;
211  itg tag, *ThrTag;
212  BucSct **ThrStk;
213 }OctThrSct;
214 
215 typedef struct
216 {
217  OctThrSct thr[ MaxThr ];
218  itg MaxLvl, NmbFreOct, NmbOct, GrdLvl, NmbBuc, NmbThr;
219  size_t MemUse;
220  fpn eps, MaxSiz, MinSiz, BucSiz, bnd[2][3];
221  OctSct oct, *CurOctBlk;
222  BucSct *grd;
223  LnkSct *NexFreLnk;
224  MemSct *NexMem;
225  MshSct *msh;
226 }TreSct;
227 
228 
229 /*----------------------------------------------------------------------------*/
230 /* Prototypes of octree procedures */
231 /*----------------------------------------------------------------------------*/
232 
233 static void SetMshBox (TreSct *, MshSct *);
234 static void AddVer (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
235 static void AddEdg (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
236 static void AddTri (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
237 static void AddQad (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
238 static void AddTet (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
239 static void SubOct (MshSct *, TreSct *, OctSct *, fpn *, fpn *);
240 static void LnkItm (TreSct *, OctSct *, itg, itg, char);
241 static OctSct *GetCrd (OctSct *, itg, fpn *, fpn *, fpn *);
242 static void GetBox (TreSct *, OctSct *, itg, itg *, itg, itg *,
243  char *, fpn [2][3], fpn, fpn *, fpn *, itg );
244 static itg BoxIntBox (fpn [2][3], fpn [2][3], fpn);
245 static void SetItm (MshSct *, itg, itg, itg, itg);
246 static void AniTri (MshSct *, itg);
247 static void SetSonCrd (itg, fpn *, fpn *, fpn *, fpn *);
248 static void GetOctLnk (MshSct *, itg, fpn *, itg *, fpn *, OctSct *,
249  fpn *, fpn *, itg (void *, itg), void *, itg);
250 static void IntRayOct (TreSct *, MshSct *, fpn *, fpn *, itg *, fpn *,
251  OctSct *, fpn *, fpn *, itg (void *, itg), void *,
252  itg);
253 static void GetBucBox (TreSct *, BucSct *, fpn *, fpn *);
254 static BucSct *GetBucNgb (TreSct *, BucSct *, itg);
255 static fpn DisVerOct (fpn *, fpn *, fpn *);
256 static itg VerInsOct (fpn *, fpn *, fpn *);
257 static char *GetPtrItm (MshSct *, itg, itg);
258 static void BakMshItm (MshSct *);
259 static void RstMshItm (MshSct *);
260 
261 
262 /*----------------------------------------------------------------------------*/
263 /* Prototypes of meshing procedures */
264 /*----------------------------------------------------------------------------*/
265 
266 static itg EdgIntEdg (EdgSct *, EdgSct *, VerSct *, fpn);
267 static fpn DisVerTri (MshSct *, fpn *, TriSct *);
268 static fpn DisVerQad (MshSct *, fpn *, QadSct *);
269 static fpn DisVerTet (MshSct *, fpn *, TetSct *);
270 static fpn GetTriSrf (TriSct *);
271 static fpn GetVolTet (TetSct *);
272 static fpn DisVerEdg (fpn *, EdgSct *);
273 static void GetTriVec (TriSct *, fpn *);
274 static void SetTriNrm (TriSct *);
275 static void SetTmpHex (HexSct *, fpn *, fpn *);
276 static itg VerInsTet (VerSct *, TetSct *, fpn);
277 static itg VerInsHex (VerSct *, HexSct *);
278 static itg EdgIntHex (EdgSct *, HexSct *, fpn);
279 static itg TriIntHex (TriSct *, HexSct *, fpn);
280 static itg QadIntHex (QadSct *, HexSct *, fpn);
281 static itg TetIntHex (TetSct *, HexSct *, fpn);
282 static itg EdgIntQad (HexSct *, itg, EdgSct *, VerSct *, fpn);
283 static itg EdgIntTri (TriSct *, EdgSct *, VerSct *, fpn);
284 static itg VerInsTri (TriSct *, VerSct *, fpn);
285 static itg VerInsEdg (EdgSct *, VerSct *, fpn);
286 static void SetEdgTng (EdgSct *);
287 static fpn GetTriAni (TriSct *);
288 
289 
290 /*----------------------------------------------------------------------------*/
291 /* Prototypes of geometric procedures */
292 /*----------------------------------------------------------------------------*/
293 
294 static void PrjVerLin (fpn *, fpn *, fpn *, fpn *);
295 static fpn PrjVerPla (fpn *, fpn *, fpn *, fpn *);
296 static void LinCmbVec3 (fpn, fpn *, fpn, fpn *, fpn *);
297 static void ClrVec (fpn *);
298 static void CpyVec (fpn *, fpn *);
299 static void AddVec2 (fpn *, fpn *);
300 static void SubVec2 (fpn *, fpn *);
301 static void SubVec3 (fpn *, fpn *, fpn *);
302 static void AddScaVec1 (fpn, fpn *);
303 static void AddScaVec2 (fpn, fpn *, fpn *);
304 static void MulVec1 (fpn, fpn *);
305 static void MulVec2 (fpn, fpn *, fpn *);
306 static void NrmVec (fpn *);
307 static void CrsPrd (fpn *, fpn *, fpn *);
308 static fpn DotPrd (fpn *, fpn *);
309 static fpn dis (fpn *, fpn *);
310 static fpn DisPow (fpn *, fpn *);
311 static fpn DisVerPla (fpn *, fpn *, fpn *);
312 static fpn GetNrmVec (fpn *);
313 static fpn VerInsBox (fpn *, fpn *, fpn *, fpn);
314 static itg LinIntBox (fpn *, fpn *, fpn *, fpn *, fpn);
315 static void LinIntPla (fpn *, fpn *, fpn *, fpn *, fpn *);
316 
317 
318 /*----------------------------------------------------------------------------*/
319 /* Prototypes of memory handling procedures */
320 /*----------------------------------------------------------------------------*/
321 
322 static void *NewMem (TreSct *, size_t);
323 static void FreAllMem (TreSct *);
324 
325 
326 /*----------------------------------------------------------------------------*/
327 /* Local mesh connectivity tables */
328 /*----------------------------------------------------------------------------*/
329 
330 static const itg TetEdg[6][2] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
331 static const itg TetFac[4][3] = { {1,2,3}, {2,0,3}, {3,0,1}, {0,2,1} };
332 static const itg TetFacEdg[4][3] = { {5,4,3}, {2,5,1}, {0,4,2}, {3,1,0} };
333 static const itg tvpe[12][2] = { {3,2}, {0,1}, {4,5}, {7,6}, {3,7}, {2,6},
334  {1,5}, {0,4}, {3,0}, {7,4}, {6,5}, {2,1} };
335 static const itg tvpf[6][4] = { {3,0,4,7}, {5,1,2,6}, {3,2,1,0},
336  {5,6,7,4},{3,7,6,2}, {5,4,0,1} };
337 
338 
339 /*----------------------------------------------------------------------------*/
340 /* Allocate and build a new octree from user's data */
341 /*----------------------------------------------------------------------------*/
342 
343 int64_t LolNewOctree(itg NmbVer, fpn *PtrCrd1, fpn *PtrCrd2,
344  itg NmbEdg, itg *PtrEdg1, itg *PtrEdg2,
345  itg NmbTri, itg *PtrTri1, itg *PtrTri2,
346  itg NmbQad, itg *PtrQad1, itg *PtrQad2,
347  itg NmbTet, itg *PtrTet1, itg *PtrTet2,
348  itg NmbPyr, itg *PtrPyr1, itg *PtrPyr2,
349  itg NmbPri, itg *PtrPri1, itg *PtrPri2,
350  itg NmbHex, itg *PtrHex1, itg *PtrHex2,
351  itg BasIdx, itg NmbThr)
352 {
353  itg i, j, k, t, EdgIdx, TotItmCnt = 0, MaxItmCnt, idx = 0;
354  fpn crd[3];
355  BucSct *buc;
356  TreSct *tre = NULL;
357  MshSct *msh;
358  OctThrSct *OctThr;
359  MshThrSct *MshThr;
360 
361  // Make sure we have a vertex table and a valid base index (0 or 1)
362  if(!NmbVer || !PtrCrd1 || !PtrCrd2 || (BasIdx < 0) || (BasIdx > 1) )
363  return(0);
364 
365  NmbThr = MAX(NmbThr, 1);
366  NmbThr = MIN(NmbThr, MaxThr);
367  MaxItmCnt = NmbVer;
368 
369  // Setup a single octant octree
370  tre = calloc(1, sizeof(TreSct));
371  assert(tre);
372 
373  // Setup the mesh structure
374  msh = NewMem(tre, sizeof(MshSct));
375  msh->BasIdx = BasIdx;
376 
377  msh->NmbItm[ LolTypVer ] = NmbVer;
378  TotItmCnt += NmbVer;
379  msh->UsrPtr[ LolTypVer ] = (char *)PtrCrd1;
380  msh->UsrSiz[ LolTypVer ] = (char *)PtrCrd2 - (char *)PtrCrd1;
381 
382  msh->NmbItm[ LolTypEdg ] = NmbEdg;
383  MaxItmCnt = MAX(MaxItmCnt, NmbEdg);
384  TotItmCnt += NmbEdg;
385  msh->UsrPtr[ LolTypEdg ] = (char *)PtrEdg1;
386  msh->UsrSiz[ LolTypEdg ] = (char *)PtrEdg2 - (char *)PtrEdg1;
387 
388  msh->NmbItm[ LolTypTri ] = NmbTri;
389  MaxItmCnt = MAX(MaxItmCnt, NmbTri);
390  TotItmCnt += NmbTri;
391  msh->UsrPtr[ LolTypTri ] = (char *)PtrTri1;
392  msh->UsrSiz[ LolTypTri ] = (char *)PtrTri2 - (char *)PtrTri1;
393 
394  msh->NmbItm[ LolTypQad ] = NmbQad;
395  MaxItmCnt = MAX(MaxItmCnt, NmbQad);
396  TotItmCnt += NmbQad;
397  msh->UsrPtr[ LolTypQad ] = (char *)PtrQad1;
398  msh->UsrSiz[ LolTypQad ] = (char *)PtrQad2 - (char *)PtrQad1;
399 
400  msh->NmbItm[ LolTypTet ] = NmbTet;
401  MaxItmCnt = MAX(MaxItmCnt, NmbTet);
402  TotItmCnt += NmbTet;
403  msh->UsrPtr[ LolTypTet ] = (char *)PtrTet1;
404  msh->UsrSiz[ LolTypTet ] = (char *)PtrTet2 - (char *)PtrTet1;
405 
406  for(t=0;t<NmbThr;t++)
407  {
408  MshThr = &msh->thr[t];
409 
410  MshThr->FlgTab = NewMem(tre, (MaxItmCnt + 1) * sizeof(char) );
411  memset(MshThr->FlgTab, 0, (MaxItmCnt + 1) * sizeof(char));
412 
413  MshThr->TagTab = NewMem(tre, (MaxItmCnt + 1) * sizeof(itg) );
414  memset(MshThr->TagTab, 0, (MaxItmCnt + 1) * sizeof(itg));
415 
416  MshThr->tag = 0;
417  }
418 
419  // Setup the octree structure
420  SetMshBox(tre, msh);
421  tre->msh = msh;
422  msh->eps = tre->eps;
423 
424  // Set the grid size depending on the number of entities in the mesh
425  tre->GrdLvl = MAX(MinGrdLvl, log((TotItmCnt) / ItmPerBuc) / (3 * log(2)));
426  tre->NmbBuc = 1<<tre->GrdLvl;
427  tre->grd = NewMem(tre, CUB(tre->NmbBuc) * sizeof(BucSct));
428  tre->NmbThr = NmbThr;
429 
430  for(t=0;t<NmbThr;t++)
431  {
432  MshThr = &msh->thr[t];
433  OctThr = &tre->thr[t];
434 
435  OctThr->ThrStk = NewMem(tre, CUB(tre->NmbBuc) * sizeof(void *));
436  OctThr->ThrTag = NewMem(tre, CUB(tre->NmbBuc) * sizeof(itg));
437 
438  // Setup the temporary edge for local geometric calculations
439  for(i=0;i<2;i++)
440  MshThr->edg.ver[i] = &MshThr->ver[i];
441 
442  // Setup the temporary triangle for local geometric calculations
443  for(i=0;i<3;i++)
444  {
445  MshThr->tri.ver[i] = &MshThr->ver[i];
446  MshThr->tri.edg[i].ver[0] = &MshThr->ver[ (i+1)%3 ];
447  MshThr->tri.edg[i].ver[1] = &MshThr->ver[ (i+2)%3 ];
448  }
449 
450  // Setup the temporary quad for local geometric calculations
451  for(i=0;i<4;i++)
452  {
453  MshThr->qad.ver[i] = &MshThr->ver[i];
454  MshThr->qad.edg[i].ver[0] = &MshThr->ver[ i ];
455  MshThr->qad.edg[i].ver[1] = &MshThr->ver[ (i+1)%4 ];
456  }
457 
458  // The quad is made of two triangles
459  MshThr->qad.tri[0].ver[0] = &MshThr->ver[0];
460  MshThr->qad.tri[0].ver[1] = &MshThr->ver[1];
461  MshThr->qad.tri[0].ver[2] = &MshThr->ver[2];
462 
463  MshThr->qad.tri[1].ver[0] = &MshThr->ver[3];
464  MshThr->qad.tri[1].ver[1] = &MshThr->ver[2];
465  MshThr->qad.tri[1].ver[2] = &MshThr->ver[1];
466 
467  for(i=0;i<2;i++)
468  for(j=0;j<3;j++)
469  {
470  MshThr->qad.tri[i].edg[j].ver[0] = MshThr->qad.tri[i].ver[ (j+1)%3 ];
471  MshThr->qad.tri[i].edg[j].ver[1] = MshThr->qad.tri[i].ver[ (j+2)%3 ];
472  }
473 
474  // Setup the temporary tetrahedron for local geometric calculations
475  for(i=0;i<4;i++)
476  MshThr->tet.ver[i] = &MshThr->ver[i];
477 
478  for(i=0;i<4;i++)
479  {
480  for(j=0;j<3;j++)
481  MshThr->tet.tri[i].ver[j] = MshThr->tet.ver[ TetFac[i][j] ];
482 
483  for(j=0;j<3;j++)
484  {
485  EdgIdx = TetFacEdg[i][j];
486  MshThr->tet.tri[i].edg[j].ver[0] = MshThr->tet.ver[ TetEdg[ EdgIdx ][0] ];
487  MshThr->tet.tri[i].edg[j].ver[1] = MshThr->tet.ver[ TetEdg[ EdgIdx ][1] ];
488  }
489  }
490 
491  for(i=0;i<6;i++)
492  for(j=0;j<2;j++)
493  MshThr->tet.edg[i].ver[j] = MshThr->tet.ver[ TetEdg[i][j] ];
494 
495  // Setup the temporary hex for local geometric calculations
496  for(i=0;i<8;i++)
497  OctThr->hex.ver[i] = &OctThr->ver[i];
498 
499  for(i=0;i<12;i++)
500  {
501  OctThr->hex.edg[i].ver[0] = &OctThr->ver[ tvpe[i][0] ];
502  OctThr->hex.edg[i].ver[1] = &OctThr->ver[ tvpe[i][1] ];
503  }
504 
505  for(i=0;i<6;i++)
506  for(j=0;j<4;j++)
507  OctThr->hex.qad[i].ver[j] = &OctThr->ver[ tvpf[i][j] ];
508 
509  for(i=0;i<4;i++)
510  {
511  OctThr->hex.edg[i].tng[0] = 1;
512  OctThr->hex.edg[i].tng[1] = 0;
513  OctThr->hex.edg[i].tng[2] = 0;
514 
515  OctThr->hex.edg[i+4].tng[0] = 0;
516  OctThr->hex.edg[i+4].tng[1] = 1;
517  OctThr->hex.edg[i+4].tng[2] = 0;
518 
519  OctThr->hex.edg[i+8].tng[0] = 0;
520  OctThr->hex.edg[i+8].tng[1] = 0;
521  OctThr->hex.edg[i+8].tng[2] = 1;
522  }
523 
524  OctThr->hex.qad[0].nrm[0] = 1;
525  OctThr->hex.qad[0].nrm[1] = 0;
526  OctThr->hex.qad[0].nrm[2] = 0;
527 
528  OctThr->hex.qad[1].nrm[0] = -1;
529  OctThr->hex.qad[1].nrm[1] = 0;
530  OctThr->hex.qad[1].nrm[2] = 0;
531 
532  OctThr->hex.qad[2].nrm[0] = 0;
533  OctThr->hex.qad[2].nrm[1] = 1;
534  OctThr->hex.qad[2].nrm[2] = 0;
535 
536  OctThr->hex.qad[3].nrm[0] = 0;
537  OctThr->hex.qad[3].nrm[1] = -1;
538  OctThr->hex.qad[3].nrm[2] = 0;
539 
540  OctThr->hex.qad[4].nrm[0] = 0;
541  OctThr->hex.qad[4].nrm[1] = 0;
542  OctThr->hex.qad[4].nrm[2] = 1;
543 
544  OctThr->hex.qad[5].nrm[0] = 0;
545  OctThr->hex.qad[5].nrm[1] = 0;
546  OctThr->hex.qad[5].nrm[2] = -1;
547  }
548 
549  // Insert each vertex in the octree
550  for(i=0;i<msh->NmbItm[ LolTypVer ];i++)
551  {
552  SetItm(msh, LolTypVer, i + BasIdx, 0, 0);
553  AddVer(msh, tre, &tre->oct, tre->bnd[0], tre->bnd[1]);
554  }
555 
556  // Insert each edge in the octree
557  for(i=0;i<msh->NmbItm[ LolTypEdg ];i++)
558  {
559  SetItm(msh, LolTypEdg, i + BasIdx, 0, 0);
560  AddEdg(msh, tre, &tre->oct, tre->bnd[0], tre->bnd[1]);
561  }
562 
563  // Insert each triangle in the octree
564  for(i=0;i<msh->NmbItm[ LolTypTri ];i++)
565  {
566  SetItm(msh, LolTypTri, i + BasIdx, TngFlg | AniFlg, 0);
567  AddTri(msh, tre, &tre->oct, tre->bnd[0], tre->bnd[1]);
568  }
569 
570  // Insert each quad in the octree
571  for(i=0;i<msh->NmbItm[ LolTypQad ];i++)
572  {
573  SetItm(msh, LolTypQad, i + BasIdx, TngFlg | AniFlg, 0);
574  AddQad(msh, tre, &tre->oct, tre->bnd[0], tre->bnd[1]);
575  }
576 
577  // Insert each tetrahedron in the octree
578  for(i=0;i<msh->NmbItm[ LolTypTet ];i++)
579  {
580  SetItm(msh, LolTypTet, i + BasIdx, TngFlg, 0);
581  AddTet(msh, tre, &tre->oct, tre->bnd[0], tre->bnd[1]);
582  }
583 
584  // Setup an acceleration grid whose buckets point to an octant
585  tre->BucSiz = (tre->bnd[1][0] - tre->bnd[0][0]) / (fpn)tre->NmbBuc;
586  crd[0] = tre->bnd[0][0] + tre->BucSiz / 2.;
587 
588  for(i=0;i<tre->NmbBuc;i++)
589  {
590  crd[1] = tre->bnd[0][1] + tre->BucSiz / 2.;
591 
592  for(j=0;j<tre->NmbBuc;j++)
593  {
594  crd[2] = tre->bnd[0][2] + tre->BucSiz / 2.;
595 
596  for(k=0;k<tre->NmbBuc;k++)
597  {
598  buc = &tre->grd[ i * POW(tre->NmbBuc) + j * tre->NmbBuc + k ];
599  buc->oct = GetCrd(&tre->oct, tre->GrdLvl,
600  crd, tre->bnd[0], tre->bnd[1]);
601  buc->pos[0] = i;
602  buc->pos[1] = j;
603  buc->pos[2] = k;
604  buc->idx = idx++;
605  crd[2] += tre->BucSiz;
606  }
607 
608  crd[1] += tre->BucSiz;
609  }
610 
611  crd[0] += tre->BucSiz;
612  }
613 
614  // Return the octree's unique ID
615  return((int64_t)tre);
616 }
617 
618 
619 /*----------------------------------------------------------------------------*/
620 /* Free the octants and the links */
621 /*----------------------------------------------------------------------------*/
622 
623 size_t LolFreeOctree(int64_t OctIdx)
624 {
625  TreSct *tre = (TreSct *)OctIdx;
626  size_t MemUse = tre->MemUse;
627 
628  FreAllMem(tre);
629  memset(tre, 0, sizeof(TreSct));
630 
631  return(MemUse);
632 }
633 
634 
635 /*----------------------------------------------------------------------------*/
636 /* Search the octree for mesh entities included in this box */
637 /*----------------------------------------------------------------------------*/
638 
639 itg LolGetBoundingBox( int64_t OctIdx, itg typ, itg MaxItm, itg *ItmTab,
640  fpn MinCrd[3], fpn MaxCrd[3], itg ThrIdx )
641 {
642  itg i, NmbItm = 0;
643  fpn box[2][3] = { {MinCrd[0], MinCrd[1], MinCrd[2]},
644  {MaxCrd[0], MaxCrd[1], MaxCrd[2]} };
645  TreSct *tre = (TreSct *)OctIdx;
646 
647  GetBox( tre, &tre->oct, typ, &NmbItm, MaxItm, ItmTab,
648  tre->msh->thr[ ThrIdx ].FlgTab, box, tre->eps,
649  tre->bnd[0], tre->bnd[1], ThrIdx );
650 
651  for(i=0;i<NmbItm;i++)
652  tre->msh->thr[ ThrIdx ].FlgTab[ ItmTab[i] ] = 0;
653 
654  return(NmbItm);
655 }
656 
657 
658 /*----------------------------------------------------------------------------*/
659 /* Recusrsive coordinates search */
660 /*----------------------------------------------------------------------------*/
661 
662 static OctSct *GetCrd( OctSct *oct, itg MaxLvl, fpn VerCrd[3],
663  fpn MinCrd[3], fpn MaxCrd[3] )
664 {
665  itg SonIdx;
666  fpn MidCrd[3], OctMin[3], OctMax[3], SonMin[3], SonMax[3];
667 
668  CpyVec(MinCrd, OctMin);
669  CpyVec(MaxCrd, OctMax);
670 
671  // Search for the smallest octant containing
672  // this vertex but stop at the grid level
673  while(oct->sub && (oct->lvl < MaxLvl))
674  {
675  LinCmbVec3(.5, OctMin, .5, OctMax, MidCrd);
676 
677  SonIdx = ((VerCrd[0] < MidCrd[0]) ? 0 : 1)
678  | ((VerCrd[1] < MidCrd[1]) ? 0 : 2)
679  | ((VerCrd[2] < MidCrd[2]) ? 0 : 4);
680 
681  SetSonCrd(SonIdx, SonMin, SonMax, OctMin, OctMax);
682  CpyVec(SonMin, OctMin);
683  CpyVec(SonMax, OctMax);
684  oct = oct->son + SonIdx;
685  }
686 
687  return(oct);
688 }
689 
690 
691 /*----------------------------------------------------------------------------*/
692 /* Recusrsive box search */
693 /*----------------------------------------------------------------------------*/
694 
695 static void GetBox( TreSct *tre, OctSct *oct, itg typ, itg *NmbItm,
696  itg MaxItm, itg *ItmTab, char *FlgTab, fpn box[2][3],
697  fpn eps, fpn MinCrd[3], fpn MaxCrd[3], itg ThrIdx )
698 {
699  itg i;
700  LnkSct *lnk;
701  HexSct hex;
702  OctThrSct *ThrOct = &tre->thr[ ThrIdx ];
703  MshThrSct *ThrMsh = &tre->msh->thr[ ThrIdx ];
704  fpn xmid = (MinCrd[0] + MaxCrd[0])/2.;
705  fpn ymid = (MinCrd[1] + MaxCrd[1])/2.;
706  fpn zmid = (MinCrd[2] + MaxCrd[2])/2.;
707  fpn son[8][2][3] = {
708  { {MinCrd[0], MinCrd[1], MinCrd[2]}, {xmid, ymid, zmid} },
709  { {xmid, MinCrd[1], MinCrd[2]}, {MaxCrd[0], ymid, zmid} },
710  { {MinCrd[0], ymid, MinCrd[2]}, {xmid, MaxCrd[1], zmid} },
711  { {xmid, ymid, MinCrd[2]}, {MaxCrd[0], MaxCrd[1], zmid} },
712  { {MinCrd[0], MinCrd[1], zmid}, {xmid, ymid, MaxCrd[2]} },
713  { {xmid, MinCrd[1], zmid}, {MaxCrd[0], ymid, MaxCrd[2]} },
714  { {MinCrd[0], ymid, zmid}, {xmid, MaxCrd[1], MaxCrd[2]} },
715  { {xmid, ymid, zmid}, {MaxCrd[0], MaxCrd[1], MaxCrd[2]} } };
716 
717  if(oct->sub)
718  {
719  // Recursively intersect the box with the octree
720  for(i=0;i<8;i++)
721  if(BoxIntBox(box, son[i], eps))
722  GetBox( tre, oct->son+i, typ, NmbItm, MaxItm, ItmTab,
723  FlgTab, box, eps, son[i][0], son[i][1], ThrIdx );
724  }
725  else if((lnk = oct->lnk) && (*NmbItm < MaxItm) )
726  {
727  SetTmpHex(&ThrOct->hex, box[0], box[1]);
728 
729  // When a terminal octant is reached, add its linked entities
730  // to the table and flag them to avoid duplicates
731  do
732  {
733  if(lnk->typ != typ)
734  continue;
735 
736  if(lnk->typ == LolTypVer)
737  {
738  SetItm(tre->msh, LolTypVer, lnk->idx, 0, ThrIdx);
739 
740  if(!VerInsHex(&ThrMsh->ver[0], &ThrOct->hex))
741  continue;
742  }
743  else if(lnk->typ == LolTypEdg)
744  {
745  SetItm(tre->msh, LolTypEdg, lnk->idx, 0, ThrIdx);
746 
747  if(!EdgIntHex(&ThrMsh->edg, &ThrOct->hex, tre->eps))
748  continue;
749  }
750  else if(lnk->typ == LolTypTri)
751  {
752  SetItm(tre->msh, LolTypTri, lnk->idx, TngFlg, ThrIdx);
753 
754  if(!TriIntHex(&ThrMsh->tri, &ThrOct->hex, tre->eps))
755  continue;
756  }
757  else if(lnk->typ == LolTypQad)
758  {
759  SetItm(tre->msh, LolTypQad, lnk->idx, TngFlg, ThrIdx);
760 
761  if(!QadIntHex(&ThrMsh->qad, &ThrOct->hex, tre->eps))
762  continue;
763  }
764  else if(lnk->typ == LolTypTet)
765  {
766  SetItm(tre->msh, LolTypTet, lnk->idx, TngFlg, ThrIdx);
767 
768  if(!TetIntHex(&ThrMsh->tet, &ThrOct->hex, tre->eps))
769  continue;
770  }
771 
772  if(!FlgTab[ lnk->idx ])
773  {
774  ItmTab[ (*NmbItm)++ ] = lnk->idx;
775  FlgTab[ lnk->idx ] = 1;
776  }
777  }while( (lnk = lnk->nex) && (*NmbItm < MaxItm) );
778  }
779 }
780 
781 
782 /*----------------------------------------------------------------------------*/
783 /* Search for the nearest item from a vertex in the octree */
784 /*----------------------------------------------------------------------------*/
785 
786 itg LolGetNearest(int64_t OctIdx, itg typ, fpn *VerCrd, fpn *MinDis, fpn MaxDis,
787  itg (UsrPrc)(void *, itg), void *UsrDat, itg ThrIdx)
788 {
789  TreSct *tre = (TreSct *)OctIdx;
790  itg i, ins = 0, out = 0, MinItm = 0, ini[3], *tag, len;
791  fpn MinCrd[3], MaxCrd[3], vec[3];
792  MshSct *msh = tre->msh;
793  BucSct *IniBuc, *buc, *ngb, **stk;
794  OctThrSct *ThrOct = &tre->thr[ ThrIdx ];
795  MshThrSct *ThrMsh = &tre->msh->thr[ ThrIdx ];
796 
797  if( (ThrIdx < 0) || (ThrIdx >= tre->NmbThr) )
798  return(0);
799 
800  ThrOct->tag++;
801  ThrMsh->tag = ThrOct->tag;
802  tag = ThrOct->ThrTag;
803  stk = ThrOct->ThrStk;
804  len = tre->NmbBuc;
805  *MinDis = (MaxDis > 0.) ? POW(MaxDis) : DBL_MAX;
806 
807  // Get the vertex's integer coordinates in the grid
808  // and clip it if it stands outside the bounding box
809  for(i=0;i<3;i++)
810  {
811  ini[i] = (VerCrd[i] - tre->bnd[0][i]) / tre->BucSiz;
812  ini[i] = MAX(0, ini[i]);
813  ini[i] = MIN(len-1, ini[i]);
814  }
815 
816  IniBuc = &tre->grd[ ini[0] * POW(len) + ini[1] * len + ini[2] ];
817 
818  // Push the octant containing the starting point on the lifo stack
819  stk[ ins++ ] = IniBuc;
820  tag[ IniBuc->idx ] = ThrOct->tag;
821 
822  // Flood fill processing of the grid :
823  // check octant's contents distance against the closest item
824  while(ins > out)
825  {
826  buc = stk[ out++ ];
827  GetBucBox( tre, buc, MinCrd, MaxCrd );
828  GetOctLnk( msh, typ, VerCrd, &MinItm, MinDis, buc->oct,
829  MinCrd, MaxCrd, UsrPrc, UsrDat, ThrIdx );
830 
831  // Push unprocessed neighbours on the stack as long as they are not too far
832  for(i=0;i<6;i++)
833  {
834  if( !(ngb = GetBucNgb(tre, buc, i)) || (tag[ ngb->idx ] == ThrOct->tag) )
835  continue;
836 
837  GetBucBox(tre, ngb, MinCrd, MaxCrd);
838 
839  if(DisVerOct(VerCrd, MinCrd, MaxCrd) < *MinDis)
840  {
841  stk[ ins++ ] = ngb;
842  tag[ ngb->idx ] = ThrOct->tag;
843  }
844  }
845  }
846 
847  *MinDis = sqrt(*MinDis);
848 
849  return(MinItm);
850 }
851 
852 
853 /*----------------------------------------------------------------------------*/
854 /* Find the closest intersected triangle along a line */
855 /*----------------------------------------------------------------------------*/
856 
857 itg LolIntersectSurface(int64_t OctIdx, fpn *VerCrd, fpn *VerTng, fpn *MinDis,
858  fpn MaxDis, itg (UsrPrc)(void *, itg), void *UsrDat,
859  itg ThrIdx )
860 {
861  TreSct *tre = (TreSct *)OctIdx;
862  OctThrSct *ThrOct = &tre->thr[ ThrIdx ];
863  MshThrSct *ThrMsh = &tre->msh->thr[ ThrIdx ];
864  itg i, ins=0, out=0, MinItm = 0, ini[3], *tag, len;
865  fpn MinCrd[3], MaxCrd[3], vec[3];
866  MshSct *msh = tre->msh;
867  BucSct *IniBuc, *buc, *ngb, **stk;
868 
869  ThrOct->tag++;
870  ThrMsh->tag = ThrOct->tag;
871  tag = ThrOct->ThrTag;
872  stk = ThrOct->ThrStk;
873  len = tre->NmbBuc;
874  *MinDis = (MaxDis > 0.) ? POW(MaxDis) : DBL_MAX;
875 
876  // Get the vertex's integer coordinates in the grid
877  // and clip it if it stands outside the bounding box
878  for(i=0;i<3;i++)
879  {
880  ini[i] = (VerCrd[i] - tre->bnd[0][i]) / tre->BucSiz;
881  ini[i] = MAX(0, ini[i]);
882  ini[i] = MIN(tre->NmbBuc-1, ini[i]);
883  }
884 
885  IniBuc = &tre->grd[ ini[0] * POW(len) + ini[1] * len + ini[2] ];
886 
887  // Push the octant containing the starting point on the lifo stack
888  stk[ ins++ ] = IniBuc;
889  tag[ IniBuc->idx ] = ThrOct->tag;
890 
891  // Flood fill processing of the grid :
892  // check octant's contents distance against the closest item
893  while(ins > out)
894  {
895  buc = stk[ out++ ];
896  GetBucBox( tre, buc, MinCrd, MaxCrd);
897  IntRayOct( tre, msh, VerCrd, VerTng, &MinItm, MinDis, buc->oct,
898  MinCrd, MaxCrd, UsrPrc, UsrDat, ThrIdx );
899 
900  // Push unprocessed neighbours intersected by the line
901  // on the stack as long as they are not too far
902  for(i=0;i<6;i++)
903  {
904  if( !(ngb = GetBucNgb(tre, buc, i)) || (tag[ ngb->idx ] == ThrOct->tag) )
905  continue;
906 
907  GetBucBox(tre, ngb, MinCrd, MaxCrd);
908 
909  if(!LinIntBox(VerCrd, VerTng, MinCrd, MaxCrd, tre->eps))
910  continue;
911 
912  if(DisVerOct(VerCrd, MinCrd, MaxCrd) < *MinDis)
913  {
914  stk[ ins++ ] = ngb;
915  tag[ ngb->idx ] = ThrOct->tag;
916  }
917  }
918  }
919 
920  *MinDis = sqrt(*MinDis);
921 
922  return(MinItm);
923 }
924 
925 
926 /*----------------------------------------------------------------------------*/
927 /* Project a vertex on a given enitity: vertex, edge, triangle or quad */
928 /*----------------------------------------------------------------------------*/
929 
930 itg LolProjectVertex(int64_t OctIdx, fpn *VerCrd, itg typ,
931  itg MinItm, fpn *MinCrd, itg ThrIdx)
932 {
933  TreSct *tre = (TreSct *)OctIdx;
934  OctThrSct *ThrOct = &tre->thr[ ThrIdx ];
935  MshThrSct *ThrMsh = &tre->msh->thr[ ThrIdx ];
936  MshSct *msh = tre->msh;
937  VerSct TmpVer;
938  TriSct *tri;
939  itg i, EdgFlg = 0;
940  fpn CurDis, MinDis = DBL_MAX;
941 
942  if(typ == LolTypVer)
943  {
944  // Vertex case, there is only one possible projection:
945  // the target vertex itself
946  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, MinItm), MinCrd);
947  return(1);
948  }
949  else if(typ == LolTypEdg)
950  {
951  // Edge case, the closest position may be on the edge itself
952  SetItm(msh, LolTypEdg, MinItm, 0, ThrIdx);
953  PrjVerLin(VerCrd, ThrMsh->edg.ver[0]->crd, ThrMsh->edg.tng, TmpVer.crd);
954 
955  if(VerInsEdg(&ThrMsh->edg, &TmpVer, tre->eps))
956  {
957  CpyVec(TmpVer.crd, MinCrd);
958  return(2);
959  }
960 
961  // Or one of its two vertices
962  if(dis(VerCrd, ThrMsh->edg.ver[0]->crd) < dis(VerCrd, ThrMsh->edg.ver[1]->crd))
963  CpyVec(ThrMsh->edg.ver[0]->crd, MinCrd);
964  else
965  CpyVec(ThrMsh->edg.ver[1]->crd, MinCrd);
966 
967  return(1);
968  }
969  else if(typ == LolTypTri)
970  {
971  // Triangle case, the closest position may be on the triangle itself
972  SetItm(msh, LolTypTri, MinItm, TngFlg, ThrIdx);
973  PrjVerPla(VerCrd, ThrMsh->tri.ver[0]->crd, ThrMsh->tri.nrm, TmpVer.crd);
974 
975  if(VerInsTri(&ThrMsh->tri, &TmpVer, tre->eps))
976  {
977  CpyVec(TmpVer.crd, MinCrd);
978  return(3); // the closest projection is inside the triangle
979  }
980 
981  // Or fall inside one of its three edges
982  for(i=0;i<3;i++)
983  {
984  PrjVerLin( VerCrd, ThrMsh->tri.edg[i].ver[0]->crd,
985  ThrMsh->tri.edg[i].tng, TmpVer.crd );
986 
987  if(VerInsEdg(&ThrMsh->tri.edg[i], &TmpVer, tre->eps)
988  && (dis(VerCrd, TmpVer.crd) < MinDis) )
989  {
990  MinDis = dis(VerCrd, TmpVer.crd);
991  CpyVec(TmpVer.crd, MinCrd);
992  EdgFlg = 2;
993  }
994  }
995 
996  // Or one of the three vertices
997  for(i=0;i<3;i++)
998  {
999  CurDis = dis(VerCrd, ThrMsh->tri.ver[i]->crd);
1000 
1001  if(CurDis < MinDis)
1002  {
1003  MinDis = CurDis;
1004  CpyVec(ThrMsh->tri.ver[i]->crd, MinCrd);
1005  EdgFlg = 0;
1006  }
1007  }
1008 
1009  if(EdgFlg)
1010  return(2); // the closest projection is on an edge
1011  else
1012  return(1); // the closest projection is on a vertex
1013  }
1014  else if(typ == LolTypQad)
1015  {
1016  SetItm(msh, LolTypQad, MinItm, TngFlg, ThrIdx);
1017 
1018  // Compute the projection on the two triangles
1019  for(i=0;i<2;i++)
1020  {
1021  tri = &ThrMsh->qad.tri[i];
1022  PrjVerPla(VerCrd, tri->ver[0]->crd, tri->nrm, TmpVer.crd);
1023 
1024  if(VerInsTri(tri, &TmpVer, tre->eps))
1025  {
1026  CpyVec(TmpVer.crd, MinCrd);
1027  return(3); // the closest projection is inside the quad
1028  }
1029  }
1030 
1031  // Check the projections on the four edges
1032  for(i=0;i<4;i++)
1033  {
1034  PrjVerLin( VerCrd, ThrMsh->qad.edg[i].ver[0]->crd,
1035  ThrMsh->qad.edg[i].tng, TmpVer.crd );
1036 
1037  if(VerInsEdg(&ThrMsh->qad.edg[i], &TmpVer, tre->eps)
1038  && (dis(VerCrd, TmpVer.crd) < MinDis) )
1039  {
1040  MinDis = dis(VerCrd, TmpVer.crd);
1041  CpyVec(TmpVer.crd, MinCrd);
1042  EdgFlg = 2;
1043  }
1044  }
1045 
1046  // Or one of the four vertices
1047  for(i=0;i<4;i++)
1048  {
1049  CurDis = dis(VerCrd, ThrMsh->qad.ver[i]->crd);
1050 
1051  if(CurDis < MinDis)
1052  {
1053  MinDis = CurDis;
1054  CpyVec(ThrMsh->qad.ver[i]->crd, MinCrd);
1055  EdgFlg = 0;
1056  }
1057  }
1058 
1059  if(EdgFlg)
1060  return(2); // the closest projection is on an edge
1061  else
1062  return(1); // the closest projection is on a vertex
1063  }
1064  else
1065  return(0);
1066 }
1067 
1068 
1069 /*----------------------------------------------------------------------------*/
1070 /* Extract the bounding box from a grid's bucket */
1071 /*----------------------------------------------------------------------------*/
1072 
1073 static void GetBucBox( TreSct *tre, BucSct *buc,
1074  fpn MinCrd[3], fpn MaxCrd[3] )
1075 {
1076  itg i;
1077 
1078  for(i=0;i<3;i++)
1079  {
1080  MinCrd[i] = tre->bnd[0][i] + buc->pos[i] * tre->BucSiz;
1081  MaxCrd[i] = tre->bnd[0][i] + (buc->pos[i]+1) * tre->BucSiz;
1082  }
1083 
1084 }
1085 
1086 
1087 /*----------------------------------------------------------------------------*/
1088 /* Get a bucket's neighbour from the grid */
1089 /*----------------------------------------------------------------------------*/
1090 
1091 static BucSct *GetBucNgb(TreSct *tre, BucSct *buc, itg dir)
1092 {
1093  if( (dir == 0) && (buc->pos[0] > 0) )
1094  return(&tre->grd[ (buc->pos[0]-1) * POW(tre->NmbBuc)
1095  + buc->pos[1] * tre->NmbBuc + buc->pos[2] ]);
1096 
1097  if( (dir == 1) && (buc->pos[0] < tre->NmbBuc-1) )
1098  return(&tre->grd[ (buc->pos[0]+1) * POW(tre->NmbBuc)
1099  + buc->pos[1] * tre->NmbBuc + buc->pos[2] ]);
1100 
1101  if( (dir == 2) && (buc->pos[1] > 0) )
1102  return(&tre->grd[ buc->pos[0] * POW(tre->NmbBuc)
1103  + (buc->pos[1]-1) * tre->NmbBuc + buc->pos[2] ]);
1104 
1105  if( (dir == 3) && (buc->pos[1] < tre->NmbBuc-1) )
1106  return(&tre->grd[ buc->pos[0] * POW(tre->NmbBuc)
1107  + (buc->pos[1]+1) * tre->NmbBuc + buc->pos[2] ]);
1108 
1109  if( (dir == 4) && (buc->pos[2] > 0) )
1110  return(&tre->grd[ buc->pos[0] * POW(tre->NmbBuc)
1111  + buc->pos[1] * tre->NmbBuc + buc->pos[2]-1 ]);
1112 
1113  if( (dir == 5) && (buc->pos[2] < tre->NmbBuc-1) )
1114  return(&tre->grd[ buc->pos[0] * POW(tre->NmbBuc)
1115  + buc->pos[1] * tre->NmbBuc + buc->pos[2]+1 ]);
1116 
1117  return(NULL);
1118 }
1119 
1120 
1121 /*----------------------------------------------------------------------------*/
1122 /* Compute the distance between a point and an octant */
1123 /*----------------------------------------------------------------------------*/
1124 
1125 static fpn DisVerOct(fpn VerCrd[3], fpn MinCrd[3], fpn MaxCrd[3])
1126 {
1127  itg i;
1128  fpn ClpCrd[3];
1129 
1130  // Project the vertex on the octant's surface
1131  for(i=0;i<3;i++)
1132  {
1133  ClpCrd[i] = MAX(VerCrd[i], MinCrd[i]);
1134  ClpCrd[i] = MIN(ClpCrd[i], MaxCrd[i]);
1135  }
1136 
1137  // The distance between the projection and the vertex is the shortest
1138  // if this latter stands OUTSIDE the octant
1139  return(DisPow(ClpCrd, VerCrd));
1140 }
1141 
1142 
1143 /*----------------------------------------------------------------------------*/
1144 /* Search for the nearest item from a vertex from an octant */
1145 /*----------------------------------------------------------------------------*/
1146 
1147 static void GetOctLnk( MshSct *msh, itg typ, fpn VerCrd[3], itg *MinItm,
1148  fpn *MinDis, OctSct *oct, fpn MinCrd[3], fpn MaxCrd[3],
1149  itg (UsrPrc)(void *, itg), void *UsrDat, itg ThrIdx )
1150 {
1151  itg i, *IdxTab;
1152  fpn CurDis, SonMin[3], SonMax[3];
1153  LnkSct *lnk;
1154  MshThrSct *ThrMsh = &msh->thr[ ThrIdx ];
1155 
1156  if(oct->sub)
1157  {
1158  // Check each sons recursively as long as the minimum distance
1159  // between the vertex and the octant is lower than
1160  // the current distance from the closest entity
1161  for(i=0;i<8;i++)
1162  {
1163  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1164 
1165  if(DisVerOct(VerCrd, SonMin, SonMax) <= *MinDis)
1166  GetOctLnk( msh, typ, VerCrd, MinItm, MinDis, oct->son + i,
1167  SonMin, SonMax, UsrPrc, UsrDat, ThrIdx );
1168  }
1169  }
1170  else if((lnk = oct->lnk))
1171  {
1172  // When a leaf octant is reached, compute the distance
1173  // between its linked enities and the vertex
1174  do
1175  {
1176  if(lnk->typ != typ)
1177  continue;
1178 
1179  if(lnk->typ == LolTypVer)
1180  CurDis = DisPow(VerCrd, (fpn *)GetPtrItm(msh, LolTypVer, lnk->idx));
1181  else if(lnk->typ == LolTypEdg)
1182  {
1183  SetItm(msh, LolTypEdg, lnk->idx, 0, ThrIdx);
1184  CurDis = DisVerEdg(VerCrd, &ThrMsh->edg);
1185  }
1186  else if(lnk->typ == LolTypTri)
1187  {
1188  if(ThrMsh->TagTab[ lnk->idx ] == ThrMsh->tag)
1189  continue;
1190 
1191  ThrMsh->TagTab[ lnk->idx ] = ThrMsh->tag;
1192 
1193  if(UsrPrc && !UsrPrc(UsrDat, lnk->idx))
1194  continue;
1195 
1196  SetItm(msh, LolTypTri, lnk->idx, 0, ThrIdx);
1197  CurDis = DisVerTri(msh, VerCrd, &ThrMsh->tri);
1198  }
1199  else if(lnk->typ == LolTypQad)
1200  {
1201  if(ThrMsh->TagTab[ lnk->idx ] == ThrMsh->tag)
1202  continue;
1203 
1204  ThrMsh->TagTab[ lnk->idx ] = ThrMsh->tag;
1205 
1206  if(UsrPrc && !UsrPrc(UsrDat, lnk->idx))
1207  continue;
1208 
1209  SetItm(msh, LolTypQad, lnk->idx, 0, ThrIdx);
1210  CurDis = DisVerQad(msh, VerCrd, &ThrMsh->qad);
1211  }
1212  else if(lnk->typ == LolTypTet)
1213  {
1214  SetItm(msh, LolTypTet, lnk->idx, 0, ThrIdx);
1215  CurDis = DisVerTet(msh, VerCrd, &ThrMsh->tet);
1216  }
1217 
1218  if(CurDis < *MinDis)
1219  {
1220  *MinItm = lnk->idx;
1221  *MinDis = CurDis;
1222  }
1223  }while((lnk = lnk->nex));
1224  }
1225 }
1226 
1227 
1228 /*----------------------------------------------------------------------------*/
1229 /* Search for the nearest triangle instersected by a line */
1230 /*----------------------------------------------------------------------------*/
1231 
1232 static void IntRayOct( TreSct *tre, MshSct *msh, fpn *crd, fpn *tng,
1233  itg *MinItm, fpn *MinDis, OctSct *oct, fpn MinCrd[3],
1234  fpn MaxCrd[3], itg (UsrPrc)(void *, itg), void *UsrDat,
1235  itg ThrIdx )
1236 {
1237  itg i, *IdxTab;
1238  fpn CurDis, SonMin[3], SonMax[3];
1239  VerSct IntVer;
1240  LnkSct *lnk;
1241  MshThrSct *ThrMsh = &tre->msh->thr[ ThrIdx ];
1242 
1243  if(oct->sub)
1244  {
1245  // Check each sons recursively as long as the minimum distance
1246  // between the vertex and the octant is lower than
1247  // the current distance from the closest entity
1248  for(i=0;i<8;i++)
1249  {
1250  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1251 
1252  if(!LinIntBox(crd, tng, SonMin, SonMax, msh->eps))
1253  continue;
1254 
1255  IntRayOct( tre, msh, crd, tng, MinItm, MinDis, oct->son+i,
1256  SonMin, SonMax, UsrPrc, UsrDat, ThrIdx );
1257  }
1258  }
1259  else if((lnk = oct->lnk))
1260  {
1261  // When a leaf octant is reached, compute the intersection
1262  // between its linked triangles and the line
1263  do
1264  {
1265  if(lnk->typ != LolTypTri)
1266  continue;
1267 
1268  if(ThrMsh->TagTab[ lnk->idx ] == ThrMsh->tag)
1269  continue;
1270 
1271  ThrMsh->TagTab[ lnk->idx ] = ThrMsh->tag;
1272  SetItm(msh, LolTypTri, lnk->idx, 0, ThrIdx);
1273 
1274  if(UsrPrc && !UsrPrc(UsrDat, lnk->idx))
1275  continue;
1276 
1277  if(DotPrd(tng, ThrMsh->tri.nrm) != 0.)
1278  {
1279  LinIntPla( crd, tng, ThrMsh->tri.ver[0]->crd,
1280  ThrMsh->tri.nrm, IntVer.crd );
1281 
1282  if(VerInsTri(&ThrMsh->tri, &IntVer, msh->eps))
1283  {
1284  CurDis = DisPow(IntVer.crd, crd);
1285 
1286  if(CurDis < *MinDis)
1287  {
1288  *MinItm = lnk->idx;
1289  *MinDis = CurDis;
1290  }
1291  }
1292  }
1293  }while((lnk = lnk->nex));
1294  }
1295 }
1296 
1297 
1298 /*----------------------------------------------------------------------------*/
1299 /* Setup an arbitrary geometrical item structure from its index */
1300 /*----------------------------------------------------------------------------*/
1301 
1302 static void SetItm(MshSct *msh, itg typ, itg idx, itg flg, itg ThrIdx)
1303 {
1304  itg i, j, *IdxTab;
1305  const itg TetEdgFac[6][2] = { {2,3}, {1,3}, {1,2}, {0,3}, {0,2}, {0,1} };
1306  MshThrSct *ThrMsh = &msh->thr[ ThrIdx ];
1307 
1308  if(typ == LolTypVer)
1309  {
1310  ThrMsh->ver[0].idx = idx;
1311  CpyVec((fpn *)GetPtrItm(msh, typ, idx), ThrMsh->ver[0].crd);
1312  }
1313  else if(typ == LolTypEdg)
1314  {
1315  // Setup the temporary edge structure with this edge's ID
1316  ThrMsh->edg.idx = idx;
1317  IdxTab = (itg *)GetPtrItm(msh, typ, idx);
1318 
1319  for(i=0;i<2;i++)
1320  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, IdxTab[i]), ThrMsh->edg.ver[i]->crd);
1321 
1322  SetEdgTng(&ThrMsh->edg);
1323  }
1324  else if(typ == LolTypTri)
1325  {
1326  // Setup the temporary triangle structure with this triangle's ID
1327  ThrMsh->tri.idx = idx;
1328  IdxTab = (itg *)GetPtrItm(msh, typ, idx);
1329 
1330  for(i=0;i<3;i++)
1331  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, IdxTab[i]), ThrMsh->tri.ver[i]->crd);
1332 
1333  SetTriNrm(&ThrMsh->tri);
1334 
1335  // Set triangle edge tangents only on request
1336  if(flg & TngFlg)
1337  for(i=0;i<3;i++)
1338  SetEdgTng(&ThrMsh->tri.edg[i]);
1339 
1340  // Compute the aspect ratio on demand
1341  if(flg & AniFlg)
1342  ThrMsh->tri.ani = GetTriAni(&ThrMsh->tri);
1343  }
1344  else if(typ == LolTypQad)
1345  {
1346  // Setup the temporary triangle structure with this triangle's ID
1347  ThrMsh->qad.idx = idx;
1348  IdxTab = (itg *)GetPtrItm(msh, typ, idx);
1349 
1350  for(i=0;i<4;i++)
1351  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, IdxTab[i]), ThrMsh->qad.ver[i]->crd);
1352 
1353  // Set quad edge tangents only on request
1354  if(flg & TngFlg)
1355  for(i=0;i<4;i++)
1356  SetEdgTng(&ThrMsh->qad.edg[i]);
1357 
1358  for(i=0;i<2;i++)
1359  {
1360  SetTriNrm(&ThrMsh->qad.tri[i]);
1361 
1362  // Set triangle edge tangents only on request
1363  if(flg & TngFlg)
1364  for(j=0;j<3;j++)
1365  SetEdgTng(&ThrMsh->qad.tri[i].edg[j]);
1366 
1367  // Compute the aspect ratio on demand
1368  if(flg & AniFlg)
1369  ThrMsh->tri.ani = GetTriAni(&ThrMsh->qad.tri[i]);
1370  }
1371  }
1372  else if(typ == LolTypTet)
1373  {
1374  // Setup the temporary tetrahedron structure with this tet ID
1375  ThrMsh->tet.idx = idx;
1376  IdxTab = (itg *)GetPtrItm(msh, typ, idx);
1377 
1378  for(i=0;i<4;i++)
1379  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, IdxTab[i]), ThrMsh->tet.ver[i]->crd);
1380 
1381  for(i=0;i<4;i++)
1382  SetTriNrm(&ThrMsh->tet.tri[i]);
1383 
1384  // Set tet edge tangents only on request
1385  if(flg & TngFlg)
1386  for(i=0;i<6;i++)
1387  {
1388  SetEdgTng(&ThrMsh->tet.edg[i]);
1389  CpyVec(ThrMsh->tet.edg[i].tng, ThrMsh->tet.tri[ TetEdgFac[i][0] ].edg[0].tng);
1390  CpyVec(ThrMsh->tet.edg[i].tng, ThrMsh->tet.tri[ TetEdgFac[i][1] ].edg[1].tng);
1391  }
1392  }
1393 }
1394 
1395 
1396 /*----------------------------------------------------------------------------*/
1397 /* Save the local mesh entities into their backup position */
1398 /*----------------------------------------------------------------------------*/
1399 
1400 static void BakMshItm(MshSct *msh)
1401 {
1402  memcpy(&msh->thr[0].BakEdg, &msh->thr[0].edg, sizeof(EdgSct));
1403  memcpy(&msh->thr[0].BakTri, &msh->thr[0].tri, sizeof(TriSct));
1404  memcpy(&msh->thr[0].BakQad, &msh->thr[0].qad, sizeof(QadSct));
1405  memcpy(&msh->thr[0].BakTet, &msh->thr[0].tet, sizeof(TetSct));
1406  memcpy( msh->thr[0].BakVer, msh->thr[0].ver, 4 * sizeof(VerSct));
1407 }
1408 
1409 
1410 /*----------------------------------------------------------------------------*/
1411 /* Restore the local entities from the backup to the current position */
1412 /*----------------------------------------------------------------------------*/
1413 
1414 static void RstMshItm(MshSct *msh)
1415 {
1416  memcpy(&msh->thr[0].edg, &msh->thr[0].BakEdg, sizeof(EdgSct));
1417  memcpy(&msh->thr[0].tri, &msh->thr[0].BakTri, sizeof(TriSct));
1418  memcpy(&msh->thr[0].qad, &msh->thr[0].BakQad, sizeof(QadSct));
1419  memcpy(&msh->thr[0].tet, &msh->thr[0].BakTet, sizeof(TetSct));
1420  memcpy( msh->thr[0].ver, msh->thr[0].BakVer, 4 * sizeof(VerSct));
1421 }
1422 
1423 
1424 /*----------------------------------------------------------------------------*/
1425 /* Get a geometric item adress in the user's table */
1426 /*----------------------------------------------------------------------------*/
1427 
1428 static char *GetPtrItm(MshSct *msh, itg typ, itg idx)
1429 {
1430  return(msh->UsrPtr[ typ ] + (idx - msh->BasIdx) * msh->UsrSiz[ typ ]);
1431 }
1432 
1433 
1434 /*----------------------------------------------------------------------------*/
1435 /* Create a mesh with one hexa enclosing the whole surface */
1436 /*----------------------------------------------------------------------------*/
1437 
1438 static void SetMshBox(TreSct *box, MshSct *msh)
1439 {
1440  itg i, j;
1441  fpn MinCrd[3], MaxCrd[3], MidCrd[3], *CrdTab, siz;
1442 
1443  // Compute the bounding box (rectangular)
1444  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, msh->BasIdx), MinCrd);
1445  CpyVec((fpn *)GetPtrItm(msh, LolTypVer, msh->BasIdx), MaxCrd);
1446 
1447  for(i=1;i<msh->NmbItm[ LolTypVer ];i++)
1448  {
1449  CrdTab = (fpn *)GetPtrItm(msh, LolTypVer, i + msh->BasIdx);
1450 
1451  for(j=0;j<3;j++)
1452  {
1453  MinCrd[j] = MIN(MinCrd[j], CrdTab[j]);
1454  MaxCrd[j] = MAX(MaxCrd[j], CrdTab[j]);
1455  }
1456  }
1457 
1458  // Compute the bounding box
1459  siz = MAX(MaxCrd[0] - MinCrd[0], MaxCrd[1] - MinCrd[1]);
1460  siz = MAX(siz, MaxCrd[2] - MinCrd[2]);
1461  box->eps = siz * FLT_EPSILON;
1462  box->MinSiz = box->MaxSiz = siz;
1463 
1464  // Move the center 1/1000th away
1465  LinCmbVec3(.5, MinCrd, .5, MaxCrd, MidCrd);
1466  AddScaVec1(siz * FLT_EPSILON, MidCrd);
1467  AddScaVec2(-1.02 * siz / 2., MidCrd, box->bnd[0]);
1468  AddScaVec2( 1.02 * siz / 2., MidCrd, box->bnd[1]);
1469 }
1470 
1471 
1472 /*----------------------------------------------------------------------------*/
1473 /* Add a vertex to leaf octants */
1474 /*----------------------------------------------------------------------------*/
1475 
1476 static void AddVer( MshSct *msh, TreSct *tre, OctSct *oct,
1477  fpn MinCrd[3], fpn MaxCrd[3] )
1478 {
1479  itg i;
1480  fpn SonMin[3], SonMax[3];
1481 
1482  if(oct->sub)
1483  {
1484  for(i=0;i<8;i++)
1485  {
1486  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1487 
1488  if(VerInsOct(msh->thr[0].ver[0].crd, SonMin, SonMax))
1489  AddVer(msh, tre, oct->son+i, SonMin, SonMax);
1490  }
1491  }
1492  else
1493  {
1494  LnkItm(tre, oct, LolTypVer, msh->thr[0].ver[0].idx, 0);
1495 
1496  if((oct->lvl < tre->GrdLvl)
1497  || ((oct->NmbVer >= oct->MaxItm) && (oct->lvl < MaxOctLvl)) )
1498  {
1499  SubOct(msh, tre, oct, MinCrd, MaxCrd);
1500  }
1501  }
1502 }
1503 
1504 
1505 /*----------------------------------------------------------------------------*/
1506 /* Add an edge to leaf octants */
1507 /*----------------------------------------------------------------------------*/
1508 
1509 static void AddEdg( MshSct *msh, TreSct *tre, OctSct *oct,
1510  fpn MinCrd[3], fpn MaxCrd[3] )
1511 {
1512  itg i;
1513  fpn SonMin[3], SonMax[3];
1514 
1515  if(oct->sub)
1516  {
1517  for(i=0;i<8;i++)
1518  {
1519  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1520  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1521 
1522  if(EdgIntHex(&msh->thr[0].edg, &tre->thr[0].hex, tre->eps))
1523  AddEdg(msh, tre, oct->son+i, SonMin, SonMax);
1524  }
1525  }
1526  else
1527  {
1528  LnkItm(tre, oct, LolTypEdg, msh->thr[0].edg.idx, 0);
1529 
1530  if( (oct->lvl < tre->GrdLvl)
1531  || ((oct->NmbEdg >= oct->MaxItm) && (oct->lvl < MaxOctLvl)) )
1532  {
1533  SubOct(msh, tre, oct, MinCrd, MaxCrd);
1534  }
1535  }
1536 }
1537 
1538 
1539 /*----------------------------------------------------------------------------*/
1540 /* Add a triangle to leaf octants */
1541 /*----------------------------------------------------------------------------*/
1542 
1543 static void AddTri( MshSct *msh, TreSct *tre, OctSct *oct,
1544  fpn MinCrd[3], fpn MaxCrd[3] )
1545 {
1546  itg i;
1547  fpn SonMin[3], SonMax[3];
1548 
1549  if(oct->sub)
1550  {
1551  for(i=0;i<8;i++)
1552  {
1553  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1554  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1555 
1556  if(TriIntHex(&msh->thr[0].tri, &tre->thr[0].hex, tre->eps))
1557  AddTri(msh, tre, oct->son+i, SonMin, SonMax);
1558  }
1559  }
1560  else
1561  {
1562  LnkItm(tre, oct, LolTypTri, msh->thr[0].tri.idx, msh->thr[0].tri.ani);
1563 
1564  if( (oct->lvl < tre->GrdLvl)
1565  || ((oct->NmbFac >= oct->MaxItm) && (oct->lvl < MaxOctLvl)) )
1566  {
1567  SubOct(msh, tre, oct, MinCrd, MaxCrd);
1568  }
1569  }
1570 }
1571 
1572 
1573 /*----------------------------------------------------------------------------*/
1574 /* Add a quad to leaf octants */
1575 /*----------------------------------------------------------------------------*/
1576 
1577 static void AddQad( MshSct *msh, TreSct *tre, OctSct *oct,
1578  fpn MinCrd[3], fpn MaxCrd[3] )
1579 {
1580  itg i;
1581  fpn SonMin[3], SonMax[3];
1582 
1583  if(oct->sub)
1584  {
1585  for(i=0;i<8;i++)
1586  {
1587  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1588  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1589 
1590  if(QadIntHex(&msh->thr[0].qad, &tre->thr[0].hex, tre->eps))
1591  AddQad(msh, tre, oct->son+i, SonMin, SonMax);
1592  }
1593  }
1594  else
1595  {
1596  LnkItm(tre, oct, LolTypQad, msh->thr[0].qad.idx, 0);
1597 
1598  if( (oct->lvl < tre->GrdLvl)
1599  || ((oct->NmbFac >= oct->MaxItm) && (oct->lvl < MaxOctLvl)) )
1600  {
1601  SubOct(msh, tre, oct, MinCrd, MaxCrd);
1602  }
1603  }
1604 }
1605 
1606 
1607 /*----------------------------------------------------------------------------*/
1608 /* Add a tetrahedron to leaf octants */
1609 /*----------------------------------------------------------------------------*/
1610 
1611 static void AddTet( MshSct *msh, TreSct *tre, OctSct *oct,
1612  fpn MinCrd[3], fpn MaxCrd[3] )
1613 {
1614  itg i;
1615  fpn SonMin[3], SonMax[3];
1616 
1617  if(oct->sub)
1618  {
1619  for(i=0;i<8;i++)
1620  {
1621  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1622  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1623 
1624  if(TetIntHex(&msh->thr[0].tet, &tre->thr[0].hex, tre->eps))
1625  AddTet(msh, tre, oct->son+i, SonMin, SonMax);
1626  }
1627  }
1628  else
1629  {
1630  LnkItm(tre, oct, LolTypTet, msh->thr[0].tet.idx, msh->thr[0].tet.ani);
1631 
1632  if( (oct->lvl < tre->GrdLvl)
1633  || ((oct->NmbFac >= oct->MaxItm) && (oct->lvl < MaxOctLvl)) )
1634  {
1635  SubOct(msh, tre, oct, MinCrd, MaxCrd);
1636  }
1637  }
1638 }
1639 
1640 
1641 /*----------------------------------------------------------------------------*/
1642 /* Subdivide an octant and its content to its sons */
1643 /*----------------------------------------------------------------------------*/
1644 
1645 static void SubOct( MshSct *msh, TreSct *tre, OctSct *oct,
1646  fpn MinCrd[3], fpn MaxCrd[3] )
1647 {
1648  itg i, j, *IdxTab;
1649  fpn SonMin[3], SonMax[3];
1650  LnkSct *lnk , *OctLnk = oct->lnk;
1651  OctSct *son;
1652 
1653  // If there is no more free octants, allocate a new bloc
1654  if(!tre->NmbFreOct)
1655  {
1656  tre->CurOctBlk = NewMem(tre, MemBlkSiz * 8 * sizeof(OctSct));
1657  tre->NmbFreOct = MemBlkSiz;
1658  }
1659 
1660  // The octant points on its first son, the other are consectutive in memory
1661  oct->son = &tre->CurOctBlk[ (MemBlkSiz - tre->NmbFreOct--) * 8 ];
1662  oct->sub = 1;
1663  tre->NmbOct+=8;
1664 
1665  // Initialize each sons
1666  for(i=0;i<8;i++)
1667  {
1668  son = oct->son+i;
1669  son->lnk = NULL;
1670  son->son = NULL;
1671  son->NmbVer = son->NmbEdg = son->NmbFac = son->NmbVol = 0;
1672  son->ani = 1;
1673  son->MaxItm = MaxItmOct;
1674  son->sub = 0;
1675  son->lvl = oct->lvl + 1;
1676  }
1677 
1678  // Update octree min octant size and max level
1679  tre->MinSiz = MIN(tre->MinSiz, (MaxCrd[0] - MinCrd[0])/2.);
1680  tre->MaxLvl = MAX(tre->MaxLvl, oct->lvl+1);
1681 
1682  // Backup the curent mesh local entities that are being tested
1683  // as the propagation process will call SetItm an all linked items
1684  BakMshItm(msh);
1685 
1686  // Now unlink every items from the father and add them to its sons
1687  while((lnk = OctLnk))
1688  {
1689  if(lnk->typ == LolTypVer)
1690  {
1691  // Check inclusion of vertices among the 8 sons
1692  SetItm(msh, LolTypVer, lnk->idx, 0, 0);
1693 
1694  for(i=0;i<8;i++)
1695  {
1696  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1697 
1698  if(VerInsOct(msh->thr[0].ver[0].crd, SonMin, SonMax))
1699  LnkItm(tre, oct->son+i, LolTypVer, lnk->idx, 0);
1700  }
1701  }
1702  else if(lnk->typ == LolTypEdg)
1703  {
1704  // Check the intersection between edge and the 8 sons
1705  SetItm(msh, LolTypEdg, lnk->idx, 0, 0);
1706 
1707  for(i=0;i<8;i++)
1708  {
1709  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1710  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1711 
1712  if(EdgIntHex(&msh->thr[0].edg, &tre->thr[0].hex, tre->eps))
1713  LnkItm(tre, oct->son+i, LolTypEdg, lnk->idx, 0);
1714  }
1715  }
1716  else if(lnk->typ == LolTypTri)
1717  {
1718  // Check the intersection between the triangle and the 8 sons
1719  SetItm(msh, LolTypTri, lnk->idx, TngFlg | AniFlg, 0);
1720 
1721  for(i=0;i<8;i++)
1722  {
1723  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1724  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1725 
1726  if(TriIntHex(&msh->thr[0].tri, &tre->thr[0].hex, tre->eps))
1727  LnkItm(tre, oct->son+i, LolTypTri, lnk->idx, oct->ani);
1728  }
1729  }
1730  else if(lnk->typ == LolTypQad)
1731  {
1732  // Check the intersection between the quad and the 8 sons
1733  SetItm(msh, LolTypQad, lnk->idx, TngFlg | AniFlg, 0);
1734 
1735  for(i=0;i<8;i++)
1736  {
1737  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1738  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1739 
1740  if(QadIntHex(&msh->thr[0].qad, &tre->thr[0].hex, tre->eps))
1741  LnkItm(tre, oct->son+i, LolTypQad, lnk->idx, oct->ani);
1742  }
1743  }
1744  else if(lnk->typ == LolTypTet)
1745  {
1746  // Check the intersection between the tet and the 8 sons
1747  SetItm(msh, LolTypTet, lnk->idx, TngFlg, 0);
1748 
1749  for(i=0;i<8;i++)
1750  {
1751  SetSonCrd(i, SonMin, SonMax, MinCrd, MaxCrd);
1752  SetTmpHex(&tre->thr[0].hex, SonMin, SonMax);
1753 
1754  if(TetIntHex(&msh->thr[0].tet, &tre->thr[0].hex, tre->eps))
1755  LnkItm(tre, oct->son+i, LolTypTet, lnk->idx, oct->ani);
1756  }
1757  }
1758 
1759  OctLnk = lnk->nex;
1760  lnk->nex = tre->NexFreLnk;
1761  tre->NexFreLnk = lnk;
1762  }
1763 
1764  // Put back the current mesh enities that are being inserted
1765  RstMshItm(msh);
1766 }
1767 
1768 
1769 /*----------------------------------------------------------------------------*/
1770 /* Add an entity to an octant linked list */
1771 /*----------------------------------------------------------------------------*/
1772 
1773 static void LnkItm(TreSct *tre, OctSct *oct, itg typ, itg idx, char ani)
1774 {
1775  itg i;
1776  LnkSct *lnk;
1777 
1778 
1779  lnk = oct->lnk;
1780  while(lnk)
1781  {
1782  if(lnk->typ == typ && lnk->idx == idx)
1783  return;
1784  lnk = lnk->nex;
1785  }
1786 
1787 
1788  // In case nore more link container are availbable, allocate a new bloc
1789  if(!tre->NexFreLnk)
1790  {
1791  tre->NexFreLnk = NewMem(tre, MemBlkSiz * sizeof(LnkSct));
1792 
1793  for(i=0;i<MemBlkSiz;i++)
1794  tre->NexFreLnk[i].nex = &tre->NexFreLnk[ i+1 ];
1795 
1796  tre->NexFreLnk[ MemBlkSiz - 1 ].nex = NULL;
1797  }
1798 
1799  // Get the next free link container and set it with the given item
1800  lnk = tre->NexFreLnk;
1801  tre->NexFreLnk = lnk->nex;
1802  lnk->typ = typ;
1803  lnk->idx = idx;
1804  lnk->nex = oct->lnk;
1805  oct->lnk = lnk;
1806 
1807  // Update the octant local counters
1808  if(typ == LolTypVer)
1809  oct->NmbVer++;
1810  else if(typ == LolTypEdg)
1811  oct->NmbEdg++;
1812  else if( (typ == LolTypTri) || (typ == LolTypQad) )
1813  {
1814  // Update the octant mean anisotropic factor
1815  // and deduce a max number of item form it
1816  if(typ == LolTypTri)
1817  {
1818  oct->ani = (oct->ani * oct->NmbFac + ani) / (oct->NmbFac + 1);
1819  oct->MaxItm = MIN(MaxItmOct * ani, 255);
1820  }
1821 
1822  oct->NmbFac++;
1823  }
1824  else
1825  oct->NmbVol++;
1826 }
1827 
1828 
1829 /*----------------------------------------------------------------------------*/
1830 /* Build an octant strucutre from the two corner points */
1831 /*----------------------------------------------------------------------------*/
1832 
1833 static void SetSonCrd( itg SonIdx, fpn SonMin[3], fpn SonMax[3],
1834  fpn MinCrd[3], fpn MaxCrd[3] )
1835 {
1836  fpn MidCrd[3];
1837 
1838  LinCmbVec3(.5, MinCrd, .5, MaxCrd, MidCrd);
1839 
1840  switch(SonIdx)
1841  {
1842  case 0 : {
1843  SonMin[0] = MinCrd[0];
1844  SonMin[1] = MinCrd[1];
1845  SonMin[2] = MinCrd[2];
1846  SonMax[0] = MidCrd[0];
1847  SonMax[1] = MidCrd[1];
1848  SonMax[2] = MidCrd[2];
1849  }return;
1850 
1851  case 1 : {
1852  SonMin[0] = MidCrd[0];
1853  SonMin[1] = MinCrd[1];
1854  SonMin[2] = MinCrd[2];
1855  SonMax[0] = MaxCrd[0];
1856  SonMax[1] = MidCrd[1];
1857  SonMax[2] = MidCrd[2];
1858  }return;
1859 
1860  case 2 : {
1861  SonMin[0] = MinCrd[0];
1862  SonMin[1] = MidCrd[1];
1863  SonMin[2] = MinCrd[2];
1864  SonMax[0] = MidCrd[0];
1865  SonMax[1] = MaxCrd[1];
1866  SonMax[2] = MidCrd[2];
1867  }return;
1868 
1869  case 3 : {
1870  SonMin[0] = MidCrd[0];
1871  SonMin[1] = MidCrd[1];
1872  SonMin[2] = MinCrd[2];
1873  SonMax[0] = MaxCrd[0];
1874  SonMax[1] = MaxCrd[1];
1875  SonMax[2] = MidCrd[2];
1876  }return;
1877 
1878  case 4 : {
1879  SonMin[0] = MinCrd[0];
1880  SonMin[1] = MinCrd[1];
1881  SonMin[2] = MidCrd[2];
1882  SonMax[0] = MidCrd[0];
1883  SonMax[1] = MidCrd[1];
1884  SonMax[2] = MaxCrd[2];
1885  }return;
1886 
1887  case 5 : {
1888  SonMin[0] = MidCrd[0];
1889  SonMin[1] = MinCrd[1];
1890  SonMin[2] = MidCrd[2];
1891  SonMax[0] = MaxCrd[0];
1892  SonMax[1] = MidCrd[1];
1893  SonMax[2] = MaxCrd[2];
1894  }return;
1895 
1896  case 6 : {
1897  SonMin[0] = MinCrd[0];
1898  SonMin[1] = MidCrd[1];
1899  SonMin[2] = MidCrd[2];
1900  SonMax[0] = MidCrd[0];
1901  SonMax[1] = MaxCrd[1];
1902  SonMax[2] = MaxCrd[2];
1903  }return;
1904 
1905  case 7 : {
1906  SonMin[0] = MidCrd[0];
1907  SonMin[1] = MidCrd[1];
1908  SonMin[2] = MidCrd[2];
1909  SonMax[0] = MaxCrd[0];
1910  SonMax[1] = MaxCrd[1];
1911  SonMax[2] = MaxCrd[2];
1912  }return;
1913  }
1914 }
1915 
1916 
1917 /*----------------------------------------------------------------------------*/
1918 /* Test if a vertex is inside an octant */
1919 /*----------------------------------------------------------------------------*/
1920 
1921 static itg VerInsOct(fpn VerCrd[3], fpn MinCrd[3], fpn MaxCrd[3])
1922 {
1923  itg i;
1924 
1925  for(i=0;i<3;i++)
1926  if( (VerCrd[i] > MaxCrd[i]) || (VerCrd[i] < MinCrd[i]) )
1927  return(0);
1928 
1929  return(1);
1930 }
1931 
1932 
1933 /*----------------------------------------------------------------------------*/
1934 /* Setup a temporary test hex from the bounding box limits */
1935 /*----------------------------------------------------------------------------*/
1936 
1937 static void SetTmpHex(HexSct *hex, fpn MinCrd[3], fpn MaxCrd[3])
1938 {
1939  hex->ver[0]->crd[0] = MinCrd[0];
1940  hex->ver[0]->crd[1] = MinCrd[1];
1941  hex->ver[0]->crd[2] = MaxCrd[2];
1942  hex->ver[1]->crd[0] = MaxCrd[0];
1943  hex->ver[1]->crd[1] = MinCrd[1];
1944  hex->ver[1]->crd[2] = MaxCrd[2];
1945  hex->ver[2]->crd[0] = MaxCrd[0];
1946  hex->ver[2]->crd[1] = MinCrd[1];
1947  hex->ver[2]->crd[2] = MinCrd[2];
1948  hex->ver[3]->crd[0] = MinCrd[0];
1949  hex->ver[3]->crd[1] = MinCrd[1];
1950  hex->ver[3]->crd[2] = MinCrd[2];
1951  hex->ver[4]->crd[0] = MinCrd[0];
1952  hex->ver[4]->crd[1] = MaxCrd[1];
1953  hex->ver[4]->crd[2] = MaxCrd[2];
1954  hex->ver[5]->crd[0] = MaxCrd[0];
1955  hex->ver[5]->crd[1] = MaxCrd[1];
1956  hex->ver[5]->crd[2] = MaxCrd[2];
1957  hex->ver[6]->crd[0] = MaxCrd[0];
1958  hex->ver[6]->crd[1] = MaxCrd[1];
1959  hex->ver[6]->crd[2] = MinCrd[2];
1960  hex->ver[7]->crd[0] = MinCrd[0];
1961  hex->ver[7]->crd[1] = MaxCrd[1];
1962  hex->ver[7]->crd[2] = MinCrd[2];
1963 }
1964 
1965 
1966 /*----------------------------------------------------------------------------*/
1967 /* Compute/test the intersection between a line and a cube */
1968 /*----------------------------------------------------------------------------*/
1969 
1970 static itg LinIntBox(fpn *LinCrd, fpn *LinTng,
1971  fpn *BoxMin, fpn *BoxMax, fpn eps)
1972 {
1973  itg i;
1974  fpn IntCrd[3], NrmTab[3][3] = { {1,0,0}, {0,1,0}, {0,0,1} };
1975 
1976  // Compute the intersections between the edge and the cube's faces
1977  for(i=0;i<3;i++)
1978  {
1979  if(LinTng[i] == 0.)
1980  continue;
1981 
1982  LinIntPla(LinCrd, LinTng, BoxMin, NrmTab[i], IntCrd);
1983 
1984  if(VerInsBox(IntCrd, BoxMin, BoxMax, eps))
1985  return(1);
1986 
1987  LinIntPla(LinCrd, LinTng, BoxMax, NrmTab[i], IntCrd);
1988 
1989  if(VerInsBox(IntCrd, BoxMin, BoxMax, eps))
1990  return(1);
1991  }
1992 
1993  return(0);
1994 }
1995 
1996 
1997 /*----------------------------------------------------------------------------*/
1998 /* Compute the intersection between a line and a plane */
1999 /*----------------------------------------------------------------------------*/
2000 
2001 static void LinIntPla( fpn *LinCrd, fpn *LinTng,
2002  fpn *PlaCrd, fpn *PlaNrm, fpn *IntCrd)
2003 {
2004  fpn u[3];
2005 
2006  SubVec3(LinCrd, PlaCrd, u);
2007  LinCmbVec3( 1., LinCrd, -DotPrd(PlaNrm, u) / DotPrd(PlaNrm, LinTng),
2008  LinTng, IntCrd);
2009 }
2010 
2011 
2012 /*----------------------------------------------------------------------------*/
2013 /* Test if a vertex stands inside a bounding box + epsilon */
2014 /*----------------------------------------------------------------------------*/
2015 
2016 static fpn VerInsBox(fpn *VerCrd, fpn *BoxMin, fpn *BoxMax, fpn eps)
2017 {
2018  if((VerCrd[0] < BoxMin[0] - eps) || (VerCrd[1] < BoxMin[1] - eps)
2019  || (VerCrd[2] < BoxMin[2] - eps) || (VerCrd[0] > BoxMax[0] + eps)
2020  || (VerCrd[1] > BoxMax[1] + eps) || (VerCrd[2] > BoxMax[2] + eps) )
2021  {
2022  return(0);
2023  }
2024 
2025  return(1);
2026 }
2027 
2028 
2029 /*----------------------------------------------------------------------------*/
2030 /* Compute/test the intersection between an edge and a hex */
2031 /*----------------------------------------------------------------------------*/
2032 
2033 static itg EdgIntHex(EdgSct *edg, HexSct *hex, fpn eps)
2034 {
2035  itg i;
2036  VerSct IntVer;
2037 
2038  // Test if an edge's vertex is included in the octant
2039  if(VerInsHex(edg->ver[0], hex) || VerInsHex(edg->ver[1], hex))
2040  return(1);
2041 
2042  // Compute the intersections between the edge and the hex's faces
2043  for(i=0;i<6;i++)
2044  if(EdgIntQad(hex, i, edg, &IntVer, eps))
2045  return(1);
2046 
2047  return(0);
2048 }
2049 
2050 
2051 /*----------------------------------------------------------------------------*/
2052 /* Test if an octant is intersected by a triangle */
2053 /*----------------------------------------------------------------------------*/
2054 
2055 static itg TriIntHex(TriSct *tri, HexSct *hex, fpn eps)
2056 {
2057  itg i, j, pos, neg;
2058  fpn CurDis;
2059  VerSct IntVer;
2060 
2061  // If there is no intersection between the bounding box
2062  // of the triangle and the octant it is no use to test,
2063  // the triangle doesn't intersect the octant
2064  if(( (tri->ver[0]->crd[0] < hex->ver[3]->crd[0])
2065  && (tri->ver[1]->crd[0] < hex->ver[3]->crd[0])
2066  && (tri->ver[2]->crd[0] < hex->ver[3]->crd[0]) )
2067  || ( (tri->ver[0]->crd[0] > hex->ver[5]->crd[0])
2068  && (tri->ver[1]->crd[0] > hex->ver[5]->crd[0])
2069  && (tri->ver[2]->crd[0] > hex->ver[5]->crd[0]) )
2070  || ( (tri->ver[0]->crd[1] < hex->ver[3]->crd[1])
2071  && (tri->ver[1]->crd[1] < hex->ver[3]->crd[1])
2072  && (tri->ver[2]->crd[1] < hex->ver[3]->crd[1]) )
2073  || ( (tri->ver[0]->crd[1] > hex->ver[5]->crd[1])
2074  && (tri->ver[1]->crd[1] > hex->ver[5]->crd[1])
2075  && (tri->ver[2]->crd[1] > hex->ver[5]->crd[1]) )
2076  || ( (tri->ver[0]->crd[2] < hex->ver[3]->crd[2])
2077  && (tri->ver[1]->crd[2] < hex->ver[3]->crd[2])
2078  && (tri->ver[2]->crd[2] < hex->ver[3]->crd[2]) )
2079  || ( (tri->ver[0]->crd[2] > hex->ver[5]->crd[2])
2080  && (tri->ver[1]->crd[2] > hex->ver[5]->crd[2])
2081  && (tri->ver[2]->crd[2] > hex->ver[5]->crd[2]) ) )
2082  {
2083  return(0);
2084  }
2085 
2086  // Test if a triangle's vertex is included in the octant
2087  for(i=0;i<3;i++)
2088  if(VerInsHex(tri->ver[i], hex))
2089  return(1);
2090 
2091  // Check whether the triangle plane intersects the octant
2092  pos = neg = 0;
2093 
2094  for(i=0;i<8;i++)
2095  {
2096  CurDis = DisVerPla(hex->ver[i]->crd, tri->ver[0]->crd, tri->nrm);
2097 
2098  if(CurDis < -eps)
2099  neg = 1;
2100  else if(CurDis > eps)
2101  pos = 1;
2102  else
2103  pos = neg = 1;
2104  }
2105 
2106  if(!pos || !neg)
2107  return(0);
2108 
2109  // Compute the intersections between the triangle edges and the hex faces
2110  for(i=0;i<6;i++)
2111  for(j=0;j<3;j++)
2112  if(EdgIntQad(hex, i, &tri->edg[j], &IntVer, eps))
2113  return(1);
2114 
2115  // Compute the intersections between the triangle and the hex edges
2116  for(i=0;i<12;i++)
2117  if(EdgIntTri(tri, &hex->edg[i], &IntVer, eps))
2118  return(1);
2119 
2120  return(0);
2121 }
2122 
2123 
2124 /*----------------------------------------------------------------------------*/
2125 /* Test if an octant is intersected by a quad */
2126 /*----------------------------------------------------------------------------*/
2127 
2128 static itg QadIntHex(QadSct *qad, HexSct *hex, fpn eps)
2129 {
2130  if(TriIntHex(&qad->tri[0], hex, eps) || TriIntHex(&qad->tri[1], hex, eps))
2131  return(1);
2132  else
2133  return(0);
2134 }
2135 
2136 
2137 /*----------------------------------------------------------------------------*/
2138 /* Test if an octant is intersected by a tetrahedron */
2139 /*----------------------------------------------------------------------------*/
2140 
2141 static itg TetIntHex(TetSct *tet, HexSct *hex, fpn eps)
2142 {
2143  itg i, j, pos, neg;
2144  fpn CurDis;
2145  VerSct IntVer;
2146 
2147  // If there is no intersection between the bounding box
2148  // of the tet and the octant it is no use to test,
2149  // the tet doesn't intersect the octant
2150  if(( (tet->ver[0]->crd[0] < hex->ver[3]->crd[0])
2151  && (tet->ver[1]->crd[0] < hex->ver[3]->crd[0])
2152  && (tet->ver[2]->crd[0] < hex->ver[3]->crd[0])
2153  && (tet->ver[3]->crd[0] < hex->ver[3]->crd[0]) )
2154  || ( (tet->ver[0]->crd[0] > hex->ver[5]->crd[0])
2155  && (tet->ver[1]->crd[0] > hex->ver[5]->crd[0])
2156  && (tet->ver[2]->crd[0] > hex->ver[5]->crd[0])
2157  && (tet->ver[3]->crd[0] > hex->ver[5]->crd[0]) )
2158  || ( (tet->ver[0]->crd[1] < hex->ver[3]->crd[1])
2159  && (tet->ver[1]->crd[1] < hex->ver[3]->crd[1])
2160  && (tet->ver[2]->crd[1] < hex->ver[3]->crd[1])
2161  && (tet->ver[3]->crd[1] < hex->ver[3]->crd[1]) )
2162  || ( (tet->ver[0]->crd[1] > hex->ver[5]->crd[1])
2163  && (tet->ver[1]->crd[1] > hex->ver[5]->crd[1])
2164  && (tet->ver[2]->crd[1] > hex->ver[5]->crd[1])
2165  && (tet->ver[3]->crd[1] > hex->ver[5]->crd[1]) )
2166  || ( (tet->ver[0]->crd[2] < hex->ver[3]->crd[2])
2167  && (tet->ver[1]->crd[2] < hex->ver[3]->crd[2])
2168  && (tet->ver[2]->crd[2] < hex->ver[3]->crd[2])
2169  && (tet->ver[3]->crd[2] < hex->ver[3]->crd[2]) )
2170  || ( (tet->ver[0]->crd[2] > hex->ver[5]->crd[2])
2171  && (tet->ver[1]->crd[2] > hex->ver[5]->crd[2])
2172  && (tet->ver[2]->crd[2] > hex->ver[5]->crd[2])
2173  && (tet->ver[3]->crd[2] > hex->ver[5]->crd[2]) ) )
2174  {
2175  return(0);
2176  }
2177 
2178  // Test if a tet's vertex is included in the octant
2179  for(i=0;i<4;i++)
2180  if(VerInsHex(tet->ver[i], hex))
2181  return(1);
2182 
2183  // Test if a oct's vertex is included in the octant
2184  for(i=0;i<8;i++)
2185  if(VerInsTet(hex->ver[i], tet, eps))
2186  return(1);
2187 
2188  // Compute the intersections between the tet edges and the hex faces
2189  for(i=0;i<6;i++)
2190  for(j=0;j<6;j++)
2191  if(EdgIntQad(hex, i, &tet->edg[j], &IntVer, eps))
2192  return(1);
2193 
2194  // Compute the intersections between the tet's face and the hex's edges
2195  for(i=0;i<12;i++)
2196  for(j=0;j<4;j++)
2197  if(EdgIntTri(&tet->tri[j], &hex->edg[i], &IntVer, eps))
2198  return(1);
2199 
2200  return(0);
2201 }
2202 
2203 
2204 /*----------------------------------------------------------------------------*/
2205 /* Test if a vertex is included in an inflated tet */
2206 /*----------------------------------------------------------------------------*/
2207 
2208 itg VerInsTet(VerSct *ver, TetSct *tet, fpn eps)
2209 {
2210  const itg TetEdg[6][2] = { {0,1}, {1,2}, {2,0}, {3,0}, {3,1}, {3,2} };
2211  const itg TetFac[4][3] = { {3,2,1}, {0,2,3}, {3,1,0}, {0,1,2} };
2212  itg i, j, ins = 1;
2213  TetSct SubTet;
2214 
2215  // Create 4 subtets from the vertex and a face of the tet.
2216  // Compute the volume of each subtet, if all of them are positive,
2217  // it means that the vertex is indide the tet
2218  SubTet.ver[3] = ver;
2219 
2220  for(i=0;i<4;i++)
2221  {
2222  for(j=0;j<3;j++)
2223  SubTet.ver[j] = tet->ver[ TetFac[i][j] ];
2224 
2225  if(GetVolTet(&SubTet) <= 0)
2226  {
2227  ins = 0;
2228  break;
2229  }
2230  }
2231 
2232  if(ins)
2233  return(1);
2234 
2235  // If the vertex is not inside the tet, it may be close to one of its faces
2236  for(i=0;i<4;i++)
2237  if(VerInsTri(&tet->tri[i], ver, eps))
2238  return(1);
2239 
2240  return(0);
2241 }
2242 
2243 
2244 /*----------------------------------------------------------------------------*/
2245 /* Test if a vertex is inside a hex */
2246 /*----------------------------------------------------------------------------*/
2247 
2248 static itg VerInsHex(VerSct *ver, HexSct *hex)
2249 {
2250  itg i;
2251 
2252  for(i=0;i<3;i++)
2253  if( (ver->crd[i] > hex->ver[5]->crd[i])
2254  || (ver->crd[i] < hex->ver[3]->crd[i]) )
2255  {
2256  return(0);
2257  }
2258 
2259  return(1);
2260 }
2261 
2262 
2263 /*----------------------------------------------------------------------------*/
2264 /* Test if an edge intersects a quad */
2265 /*----------------------------------------------------------------------------*/
2266 
2267 static itg EdgIntQad(HexSct *hex, itg FacIdx, EdgSct *edg,
2268  VerSct *IntVer, fpn eps)
2269 {
2270  itg i, NmbVer = 0;
2271  fpn sgn[2];
2272  VerSct *ver=NULL;
2273  EdgSct edg2;
2274  QadSct *qad = &hex->qad[ FacIdx ];
2275 
2276  // Compute the distance between the edge's vertices and the quad's plane
2277  for(i=0;i<2;i++)
2278  {
2279  sgn[i] = DisVerPla(edg->ver[i]->crd, qad->ver[0]->crd, qad->nrm);
2280 
2281  if(fabs(sgn[i]) < eps)
2282  {
2283  ver = edg->ver[i];
2284  NmbVer++;
2285  }
2286  }
2287 
2288  // This leads to 3 possible cases
2289  switch(NmbVer)
2290  {
2291  // Both vertices are far from the plane :
2292  // if they stand in the same side of the plane, there is no intersection,
2293  // otherwise we compute the intersection between the edge and the plane
2294  // and test wether it falls inside the quad or not
2295  case 0 :
2296  {
2297  // Test if both vertices stand on opposite sides of the plane
2298  if(sgn[0] * sgn[1] < 0)
2299  {
2300  LinCmbVec3( fabs(sgn[0]) / (fabs(sgn[0])
2301  + fabs(sgn[1])), edg->ver[1]->crd,
2302  fabs(sgn[1]) / (fabs(sgn[0])
2303  + fabs(sgn[1])), edg->ver[0]->crd, IntVer->crd);
2304 
2305  return(VerInsHex(IntVer, hex));
2306  }
2307  }break;
2308 
2309  case 1 :
2310  {
2311  if(VerInsHex(ver, hex))
2312  {
2313  CpyVec(ver->crd, IntVer->crd);
2314  return(1);
2315  }
2316  }break;
2317 
2318  case 2 :
2319  {
2320  for(i=0;i<4;i++)
2321  {
2322  edg2.ver[0] = qad->ver[i];
2323  edg2.ver[1] = qad->ver[ (i+1)%4 ];
2324  SetEdgTng(&edg2);
2325 
2326  if(EdgIntEdg(edg, &edg2, IntVer, eps))
2327  return(1);
2328  }
2329  }break;
2330  }
2331 
2332  return(0);
2333 }
2334 
2335 
2336 /*----------------------------------------------------------------------------*/
2337 /* Test if an edge intersects a triangle */
2338 /*----------------------------------------------------------------------------*/
2339 
2340 static itg EdgIntTri(TriSct *tri, EdgSct *edg, VerSct *IntVer, fpn eps)
2341 {
2342  itg i, NmbVer = 0;
2343  fpn sgn[2];
2344  VerSct *ver=NULL;
2345  EdgSct edg2;
2346 
2347  // Compute the distance between the edge's vertices and the triangle's plane
2348  for(i=0;i<2;i++)
2349  {
2350  sgn[i] = DisVerPla(edg->ver[i]->crd, tri->ver[0]->crd, tri->nrm);
2351 
2352  if(fabs(sgn[i]) < eps)
2353  {
2354  ver = edg->ver[i];
2355  NmbVer++;
2356  }
2357  }
2358 
2359  // This leads to 3 possible cases
2360  switch(NmbVer)
2361  {
2362  // Both vertices are far from the plane :
2363  // if they stand in the same side of the plane, there is no intersection,
2364  // otherwise we compute the intersection between the edge and the plane
2365  // and test wether it falls inside the triangle or not
2366  case 0 :
2367  {
2368  // Test if both vertices stand on opposite sides of the plane
2369  if(sgn[0] * sgn[1] < 0)
2370  {
2371  LinCmbVec3( fabs(sgn[0]) / (fabs(sgn[0])
2372  + fabs(sgn[1])), edg->ver[1]->crd,
2373  fabs(sgn[1]) / (fabs(sgn[0])
2374  + fabs(sgn[1])), edg->ver[0]->crd, IntVer->crd);
2375 
2376  return(VerInsTri(tri, IntVer, eps));
2377  }
2378  }break;
2379 
2380  case 1 :
2381  {
2382  if(VerInsTri(tri, ver, eps))
2383  {
2384  CpyVec(ver->crd, IntVer->crd);
2385  return(1);
2386  }
2387  }break;
2388 
2389  case 2 :
2390  {
2391  for(i=0;i<3;i++)
2392  {
2393  edg2.ver[0] = tri->ver[i];
2394  edg2.ver[1] = tri->ver[ (i+1)%3 ];
2395  SetEdgTng(&edg2);
2396 
2397  if(EdgIntEdg(edg, &edg2, IntVer, eps))
2398  return(1);
2399  }
2400  }break;
2401  }
2402 
2403  return(0);
2404 }
2405 
2406 
2407 /*----------------------------------------------------------------------------*/
2408 /* Test if a vertex is included in an inflated triangle */
2409 /*----------------------------------------------------------------------------*/
2410 
2411 static itg VerInsTri(TriSct *tri, VerSct *ver, fpn eps)
2412 {
2413  itg i, ins = 1;
2414  fpn vec[3][3], nrm[3];
2415  VerSct img;
2416  EdgSct edg;
2417 
2418  // Project the vertex on the triangle's plane and check the distance
2419  if(PrjVerPla(ver->crd, tri->ver[0]->crd, tri->nrm, img.crd) > eps)
2420  return(0);
2421 
2422  // Compare the normals of the three sub triangles against
2423  // the original triangle's normal. If one of them is opposite,
2424  // the vertex does not belong to the triangle
2425  for(i=0;i<3;i++)
2426  SubVec3(tri->ver[i]->crd, img.crd, vec[i]);
2427 
2428  for(i=0;i<3;i++)
2429  {
2430  CrsPrd(vec[ (i+1)%3 ], vec[i], nrm);
2431 
2432  if(DotPrd(nrm, tri->nrm) <= 0)
2433  {
2434  ins = 0;
2435  break;
2436  }
2437  }
2438 
2439  if(ins)
2440  return(1);
2441 
2442  // If the vertex is not inside the triangle,
2443  // it may be close to one of its edges
2444  for(i=0;i<3;i++)
2445  {
2446  edg.ver[0] = tri->ver[i];
2447  edg.ver[1] = tri->ver[ (i+1)%3 ];
2448  SetEdgTng(&edg);
2449 
2450  if(VerInsEdg(&edg, &img, eps))
2451  return(1);
2452  }
2453 
2454  return(0);
2455 }
2456 
2457 
2458 /*----------------------------------------------------------------------------*/
2459 /* Compute the intersection between coplanar edges */
2460 /*----------------------------------------------------------------------------*/
2461 
2462 static itg EdgIntEdg(EdgSct *edg1, EdgSct *edg2, VerSct *IntVer, fpn eps)
2463 {
2464  itg i, NmbVer = 0;
2465  fpn siz[2];
2466  VerSct img, *ver=NULL;
2467 
2468  // Compute the distance between the edge1's vertices
2469  // and the edge2 support line
2470  for(i=0;i<2;i++)
2471  {
2472  PrjVerLin(edg2->ver[i]->crd, edg1->ver[0]->crd, edg1->tng, img.crd);
2473  siz[i] = dis(edg2->ver[i]->crd, img.crd);
2474 
2475  if(siz[i] < eps)
2476  NmbVer++;
2477  }
2478 
2479  // This leads to 3 possible cases
2480  if(NmbVer < 2)
2481  {
2482  LinCmbVec3( siz[0]/(siz[0]+siz[1]), edg2->ver[1]->crd,
2483  siz[1]/(siz[0]+siz[1]), edg2->ver[0]->crd, IntVer->crd);
2484 
2485  return(VerInsEdg(edg1, IntVer, eps));
2486  }
2487  else
2488  {
2489  ver = NULL;
2490 
2491  if(VerInsEdg(edg2, edg1->ver[0], eps))
2492  ver = edg1->ver[0];
2493  else if(VerInsEdg(edg2, edg1->ver[1], eps))
2494  ver = edg1->ver[1];
2495  else if(VerInsEdg(edg1, edg2->ver[0], eps))
2496  ver = edg2->ver[0];
2497  else if(VerInsEdg(edg1, edg2->ver[1], eps))
2498  ver = edg2->ver[1];
2499 
2500  if(ver)
2501  {
2502  CpyVec(ver->crd, IntVer->crd);
2503  return(1);
2504  }
2505  }
2506 
2507  return(0);
2508 }
2509 
2510 
2511 /*----------------------------------------------------------------------------*/
2512 /* Test if a vertex belongs to an edge */
2513 /*----------------------------------------------------------------------------*/
2514 
2515 static itg VerInsEdg(EdgSct *edg, VerSct *ver, fpn eps)
2516 {
2517  itg i;
2518  fpn u[3], v[3];
2519  VerSct img;
2520 
2521  // Project the vertex on the edge's support line
2522  PrjVerLin(ver->crd, edg->ver[0]->crd, edg->tng, img.crd);
2523 
2524  if(DisPow(ver->crd, img.crd) > POW(eps))
2525  return(0);
2526 
2527  // Check if the image belongs to the edge
2528  SubVec3(img.crd, edg->ver[0]->crd, u);
2529  SubVec3(img.crd, edg->ver[1]->crd, v);
2530 
2531  if(DotPrd(u, v) < 0)
2532  return(1);
2533 
2534  // If the vertex does not belong to the edge,
2535  // it may be close to one of its vertices
2536  for(i=0;i<2;i++)
2537  if(DisPow(img.crd, edg->ver[i]->crd) < POW(eps))
2538  return(1);
2539 
2540  return(0);
2541 }
2542 
2543 
2544 /*----------------------------------------------------------------------------*/
2545 /* Compute the distance between a vertex and a triangle */
2546 /*----------------------------------------------------------------------------*/
2547 
2548 static fpn DisVerTri(MshSct *msh, fpn VerCrd[3], TriSct *tri)
2549 {
2550  fpn ImgCrd[3], TmpCrd[3], u[3], v[3], w[3], nrm[3], SubVol[3], TotVol;
2551 
2552  // Project the vertex on the triangle's plane
2553  PrjVerPla(VerCrd, tri->ver[0]->crd, tri->nrm, ImgCrd);
2554 
2555  // Compute the vectors stemming from the projection to the triangle's nodes
2556  SubVec3(tri->ver[0]->crd, ImgCrd, u);
2557  SubVec3(tri->ver[1]->crd, ImgCrd, v);
2558  SubVec3(tri->ver[2]->crd, ImgCrd, w);
2559 
2560  // Compute the three tets' volumes
2561  CrsPrd(v, w, nrm);
2562  SubVol[0] = -DotPrd(nrm, tri->nrm);
2563  SubVol[0] = MAX(SubVol[0], 0.);
2564 
2565  CrsPrd(w, u, nrm);
2566  SubVol[1] = -DotPrd(nrm, tri->nrm);
2567  SubVol[1] = MAX(SubVol[1], 0.);
2568 
2569  CrsPrd(u, v, nrm);
2570  SubVol[2] = -DotPrd(nrm, tri->nrm);
2571  SubVol[2] = MAX(SubVol[2], 0.);
2572 
2573  // Compute the closest position with the barycentric coordinates
2574  TotVol = SubVol[0] + SubVol[1] + SubVol[2];
2575  MulVec2(SubVol[0] / TotVol, tri->ver[0]->crd, ImgCrd);
2576  MulVec2(SubVol[1] / TotVol, tri->ver[1]->crd, TmpCrd);
2577  AddVec2(TmpCrd, ImgCrd);
2578  MulVec2(SubVol[2] / TotVol, tri->ver[2]->crd, TmpCrd);
2579  AddVec2(TmpCrd, ImgCrd);
2580 
2581  // Return the square of the distance
2582  return(DisPow(VerCrd, ImgCrd));
2583 }
2584 
2585 
2586 /*----------------------------------------------------------------------------*/
2587 /* Compute the distance between a vertex and a quad */
2588 /*----------------------------------------------------------------------------*/
2589 
2590 static fpn DisVerQad(MshSct *msh, fpn VerCrd[3], QadSct *qad)
2591 {
2592  return( MIN(DisVerTri(msh, VerCrd, &qad->tri[0]),
2593  DisVerTri(msh, VerCrd, &qad->tri[1])) );
2594 }
2595 
2596 
2597 /*----------------------------------------------------------------------------*/
2598 /* Compute the distance between a vertex and a tetrahedron */
2599 /*----------------------------------------------------------------------------*/
2600 
2601 static fpn DisVerTet(MshSct *msh, fpn *VerCrd, TetSct *tet)
2602 {
2603  itg i;
2604  fpn CurDis, MinDis = DBL_MAX;
2605  VerSct TmpVer;
2606 
2607  CpyVec(VerCrd, TmpVer.crd);
2608 
2609  if(VerInsTet(&TmpVer, tet, msh->eps))
2610  return(0.);
2611 
2612  for(i=0;i<4;i++)
2613  {
2614  CurDis = DisPow(VerCrd, tet->ver[i]->crd);
2615 
2616  if(CurDis < MinDis)
2617  MinDis = CurDis;
2618  }
2619 
2620  return(MinDis);
2621 }
2622 
2623 
2624 /*----------------------------------------------------------------------------*/
2625 /* Compute the triangle's surface */
2626 /*----------------------------------------------------------------------------*/
2627 
2628 static fpn GetTriSrf(TriSct *tri)
2629 {
2630  fpn nrm[3];
2631 
2632  // Compute the cross-product vector and get its size
2633  GetTriVec(tri, nrm);
2634  return(GetNrmVec(nrm) / 2.);
2635 }
2636 
2637 
2638 /*----------------------------------------------------------------------------*/
2639 /* fpn precision computed volume on fpn coordinates */
2640 /*----------------------------------------------------------------------------*/
2641 
2642 static fpn GetVolTet(TetSct *tet)
2643 {
2644  fpn c[9];
2645 
2646  // Create the linear system, each column is filled with a vector
2647  c[0] = tet->ver[1]->crd[0] - tet->ver[0]->crd[0];
2648  c[1] = tet->ver[2]->crd[0] - tet->ver[0]->crd[0];
2649  c[2] = tet->ver[3]->crd[0] - tet->ver[0]->crd[0];
2650 
2651  c[3] = tet->ver[1]->crd[1] - tet->ver[0]->crd[1];
2652  c[4] = tet->ver[2]->crd[1] - tet->ver[0]->crd[1];
2653  c[5] = tet->ver[3]->crd[1] - tet->ver[0]->crd[1];
2654 
2655  c[6] = tet->ver[1]->crd[2] - tet->ver[0]->crd[2];
2656  c[7] = tet->ver[2]->crd[2] - tet->ver[0]->crd[2];
2657  c[8] = tet->ver[3]->crd[2] - tet->ver[0]->crd[2];
2658 
2659  // Return the "determinant" of the matrix
2660  return( c[0] * (c[4]*c[8] - c[5]*c[7])
2661  + c[1] * (c[5]*c[6] - c[3]*c[8])
2662  + c[2] * (c[3]*c[7] - c[4]*c[6]) );
2663 }
2664 
2665 
2666 /*----------------------------------------------------------------------------*/
2667 /* Compute the distance between a vertex and an edge */
2668 /*----------------------------------------------------------------------------*/
2669 
2670 static fpn DisVerEdg(fpn VerCrd[3], EdgSct *edg)
2671 {
2672  fpn dis0, dis1, TotSiz = 0.;
2673  VerSct img;
2674 
2675  PrjVerLin(VerCrd, edg->ver[0]->crd, edg->tng, img.crd);
2676  TotSiz = dis(edg->ver[0]->crd, img.crd) + dis(edg->ver[1]->crd, img.crd);
2677 
2678  if( (TotSiz - edg->siz) < .00001 * (TotSiz + edg->siz))
2679  return(DisPow(VerCrd, img.crd));
2680 
2681  dis0 = DisPow(VerCrd, edg->ver[0]->crd);
2682  dis1 = DisPow(VerCrd, edg->ver[1]->crd);
2683 
2684  return(MIN(dis0, dis1));
2685 }
2686 
2687 
2688 /*----------------------------------------------------------------------------*/
2689 /* Compute the triangle's normal vector */
2690 /*----------------------------------------------------------------------------*/
2691 
2692 static void GetTriVec(TriSct *tri, fpn w[3])
2693 {
2694  fpn u[3], v[3];
2695 
2696  // Compute the vector product
2697  SubVec3(tri->ver[1]->crd, tri->ver[0]->crd, u);
2698  SubVec3(tri->ver[2]->crd, tri->ver[0]->crd, v);
2699  CrsPrd(v, u, w);
2700 }
2701 
2702 
2703 /*----------------------------------------------------------------------------*/
2704 /* Compute and set the triangle normal vector */
2705 /*----------------------------------------------------------------------------*/
2706 
2707 static void SetTriNrm(TriSct *tri)
2708 {
2709  // Compute the cross-product vector and normalize it
2710  GetTriVec(tri, tri->nrm);
2711  NrmVec(tri->nrm);
2712 }
2713 
2714 
2715 /*----------------------------------------------------------------------------*/
2716 /* Compute and store the edge's unit tangent */
2717 /*----------------------------------------------------------------------------*/
2718 
2719 static void SetEdgTng(EdgSct *edg)
2720 {
2721  SubVec3(edg->ver[1]->crd, edg->ver[0]->crd, edg->tng);
2722  edg->siz = GetNrmVec(edg->tng);
2723 
2724  if(edg->siz)
2725  MulVec1(1./edg->siz, edg->tng);
2726 }
2727 
2728 
2729 /*----------------------------------------------------------------------------*/
2730 /* Compute the normal projection of a point on a line */
2731 /*----------------------------------------------------------------------------*/
2732 
2733 static void PrjVerLin( fpn VerCrd[3], fpn LinCrd[3],
2734  fpn LinTng[3], fpn ImgCrd[3] )
2735 {
2736  fpn dp, u[3];
2737 
2738  // Compute the scalar product of the projected vector and the edge's vector
2739  SubVec3(VerCrd, LinCrd, u);
2740  dp = DotPrd(u, LinTng);
2741  LinCmbVec3(1., LinCrd, dp, LinTng, ImgCrd);
2742 }
2743 
2744 
2745 /*----------------------------------------------------------------------------*/
2746 /* Compute the normal projection of a vertex on a plane */
2747 /*----------------------------------------------------------------------------*/
2748 
2749 static fpn PrjVerPla(fpn VerCrd[3], fpn PlaCrd[3],
2750  fpn PlaNrm[3], fpn ImgCrd[3] )
2751 {
2752  fpn DisPla, u[3];
2753 
2754  // Compute the scalar product between the unit normal N and a vector V
2755  // defined by the vertex to project and the base plane vertex
2756  SubVec3(PlaCrd, VerCrd, u);
2757  DisPla = DotPrd(u, PlaNrm);
2758  MulVec2(DisPla, PlaNrm, ImgCrd);
2759  AddVec2(VerCrd, ImgCrd);
2760 
2761  // Return the absolute ditance on the fly
2762  return(fabs(DisPla));
2763 }
2764 
2765 
2766 /*----------------------------------------------------------------------------*/
2767 /* Compute the distance between a vertex and a plane */
2768 /*----------------------------------------------------------------------------*/
2769 
2770 static fpn DisVerPla(fpn VerCrd[3], fpn PlaCrd[3], fpn PlaNrm[3])
2771 {
2772  fpn vec[3];
2773 
2774  SubVec3(VerCrd, PlaCrd, vec);
2775  return(DotPrd(vec, PlaNrm));
2776 }
2777 
2778 
2779 /*----------------------------------------------------------------------------*/
2780 /* Test the intersection of two bounding boxex + epsilon */
2781 /*----------------------------------------------------------------------------*/
2782 
2783 static itg BoxIntBox(fpn box1[2][3], fpn box2[2][3], fpn eps)
2784 {
2785  if(( ((box1[0][0] > box2[0][0] - eps) && (box1[0][0] < box2[1][0] + eps))
2786  || ((box1[1][0] > box2[0][0] - eps) && (box1[1][0] < box2[1][0] + eps))
2787  || ((box1[0][0] < box2[0][0] ) && (box1[1][0] > box2[1][0] )))
2788  && ( ((box1[0][1] > box2[0][1] - eps) && (box1[0][1] < box2[1][1] + eps))
2789  || ((box1[1][1] > box2[0][1] - eps) && (box1[1][1] < box2[1][1] + eps))
2790  || ((box1[0][1] < box2[0][1] ) && (box1[1][1] > box2[1][1] )))
2791  && ( ((box1[0][2] > box2[0][2] - eps) && (box1[0][2] < box2[1][2] + eps))
2792  || ((box1[1][2] > box2[0][2] - eps) && (box1[1][2] < box2[1][2] + eps))
2793  || ((box1[0][2] < box2[0][2] ) && (box1[1][2] > box2[1][2] ))) )
2794  {
2795  return(1);
2796  }
2797 
2798  return(0);
2799 }
2800 
2801 
2802 /*----------------------------------------------------------------------------*/
2803 /* Compute a triangle aspect ration */
2804 /*----------------------------------------------------------------------------*/
2805 
2806 static fpn GetTriAni(TriSct *tri)
2807 {
2808  fpn srf, len, MaxLen;
2809 
2810  srf = GetTriSrf(tri);
2811  MaxLen = len = DisPow(tri->ver[0]->crd, tri->ver[1]->crd);
2812  len = DisPow(tri->ver[1]->crd, tri->ver[2]->crd);
2813  MaxLen = MAX(len, MaxLen);
2814  len = DisPow(tri->ver[2]->crd, tri->ver[0]->crd);
2815  MaxLen = MAX(len, MaxLen);
2816 
2817  return(sqrt(MaxLen / srf));
2818 }
2819 
2820 
2821 /*----------------------------------------------------------------------------*/
2822 /* Various basic operations on vectors */
2823 /*----------------------------------------------------------------------------*/
2824 
2825 // Euclidian distance
2826 static fpn dis(fpn a[3], fpn b[3])
2827 {
2828  itg i;
2829  fpn siz = 0;
2830 
2831  for(i=0;i<3;i++)
2832  siz += POW(b[i] - a[i]);
2833 
2834  return(sqrt(siz));
2835 }
2836 
2837 // Euclidian distance to the power of two
2838 static fpn DisPow(fpn a[3], fpn b[3])
2839 {
2840  itg i;
2841  fpn siz = 0;
2842 
2843  for(i=0;i<3;i++)
2844  siz += POW(b[i] - a[i]);
2845 
2846  return(siz);
2847 }
2848 
2849 // V = V - U
2850 static void SubVec2(fpn u[3], fpn v[3])
2851 {
2852  itg i;
2853 
2854  for(i=0;i<3;i++)
2855  v[i] -= u[i];
2856 }
2857 
2858 // W = U - V
2859 static void SubVec3(fpn u[3], fpn v[3], fpn w[3])
2860 {
2861  itg i;
2862 
2863  for(i=0;i<3;i++)
2864  w[i] = u[i] - v[i];
2865 }
2866 
2867 // Euclidian norm
2868 static void NrmVec(fpn u[3])
2869 {
2870  itg i;
2871  fpn dp = 0;
2872 
2873  for(i=0;i<3;i++)
2874  dp += u[i] * u[i];
2875 
2876  if(dp < DBL_MIN)
2877  return;
2878 
2879  dp = 1. / sqrt(dp);
2880 
2881  for(i=0;i<3;i++)
2882  u[i] *= dp;
2883 }
2884 
2885 // Dot Product
2886 static fpn DotPrd(fpn u[3], fpn v[3])
2887 {
2888  itg i;
2889  fpn dp = 0;
2890 
2891  for(i=0;i<3;i++)
2892  dp += u[i] * v[i];
2893 
2894  return(dp);
2895 }
2896 
2897 // Cross product
2898 static void CrsPrd(fpn u[3], fpn v[3], fpn w[3])
2899 {
2900  w[0] = u[1] * v[2] - u[2] * v[1];
2901  w[1] = u[2] * v[0] - u[0] * v[2];
2902  w[2] = u[0] * v[1] - u[1] * v[0];
2903 }
2904 
2905 // Linear combinaison: W = a*U + b*V
2906 static void LinCmbVec3(fpn w1, fpn v1[3], fpn w2, fpn v2[3], fpn v3[3])
2907 {
2908  itg i;
2909 
2910  for(i=0;i<3;i++)
2911  v3[i] = w1 * v1[i] + w2 * v2[i];
2912 }
2913 
2914 // U = 0
2915 static void ClrVec(fpn u[3])
2916 {
2917  itg i;
2918 
2919  for(i=0;i<3;i++)
2920  u[i] = 0.;
2921 }
2922 
2923 // V = U
2924 static void CpyVec(fpn u[3], fpn v[3])
2925 {
2926  itg i;
2927 
2928  for(i=0;i<3;i++)
2929  v[i] = u[i];
2930 }
2931 
2932 // V = V + U
2933 static void AddVec2(fpn u[3], fpn v[3])
2934 {
2935  itg i;
2936 
2937  for(i=0;i<3;i++)
2938  v[i] += u[i];
2939 }
2940 
2941 // U = U + [s;s;s]
2942 static void AddScaVec1(fpn s, fpn u[3])
2943 {
2944  itg i;
2945 
2946  for(i=0;i<3;i++)
2947  u[i] += s;
2948 }
2949 
2950 // V = U + [s;s;s]
2951 static void AddScaVec2(fpn s, fpn u[3], fpn v[3])
2952 {
2953  itg i;
2954 
2955  for(i=0;i<3;i++)
2956  v[i] = u[i] + s;
2957 }
2958 
2959 // U = w*U
2960 static void MulVec1(const fpn w, fpn u[3])
2961 {
2962  u[0] *= w;
2963  u[1] *= w;
2964  u[2] *= w;
2965 }
2966 
2967 // V = w*U
2968 static void MulVec2(fpn w, fpn u[3], fpn v[3])
2969 {
2970  itg i;
2971 
2972  for(i=0;i<3;i++)
2973  v[i] = w * u[i];
2974 }
2975 
2976 // Return Euclidian norm
2977 static fpn GetNrmVec(fpn u[3])
2978 {
2979  return(sqrt(POW(u[0]) + POW(u[1]) + POW(u[2])));
2980 }
2981 
2982 
2983 /*----------------------------------------------------------------------------*/
2984 /* Allocate a chunk of memory and link it to memory allocator's list */
2985 /*----------------------------------------------------------------------------*/
2986 
2987 static void *NewMem(TreSct *tre, size_t siz)
2988 {
2989  MemSct *mem;
2990 
2991  mem = malloc(sizeof(MemSct));
2992  assert(mem);
2993  mem->adr = malloc(siz);
2994  assert(mem->adr);
2995  mem->siz = siz;
2996  mem->nex = tre->NexMem;
2997  tre->NexMem = mem;
2998  tre->MemUse += siz;
2999 
3000  return(mem->adr);
3001 }
3002 
3003 
3004 /*----------------------------------------------------------------------------*/
3005 /* Cylct through the linked list of allocated block and free them all */
3006 /*----------------------------------------------------------------------------*/
3007 
3008 static void FreAllMem(TreSct *tre)
3009 {
3010  MemSct *mem = tre->NexMem, *nex;
3011 
3012  while(mem)
3013  {
3014  tre->MemUse -= mem->siz;
3015  nex = mem->nex;
3016  free(mem->adr);
3017  free(mem);
3018  mem = nex;
3019  }
3020 }
3021 
3022 
3023 /*----------------------------------------------------------------------------*/
3024 /* Fortran 77 API */
3025 /*----------------------------------------------------------------------------*/
3026 
3027 int64_t call(lolnewoctree)(itg *NmbVer, fpn *VerTab1, fpn *VerTab2,
3028  itg *NmbEdg, itg *EdgTab1, itg *EdgTab2,
3029  itg *NmbTri, itg *TriTab1, itg *TriTab2,
3030  itg *NmbQad, itg *QadTab1, itg *QadTab2,
3031  itg *NmbTet, itg *TetTab1, itg *TetTab2,
3032  itg *NmbPyr, itg *PyrTab1, itg *PyrTab2,
3033  itg *NmbPri, itg *PriTab1, itg *PriTab2,
3034  itg *NmbHex, itg *HexTab1, itg *HexTab2,
3035  itg *BasIdx, itg *NmbThr)
3036 {
3037  return(LolNewOctree( *NmbVer, VerTab1, VerTab2, *NmbEdg, EdgTab1, EdgTab2,
3038  *NmbTri, TriTab1, TriTab2, *NmbQad, QadTab1, QadTab2,
3039  *NmbTet, TetTab1, TetTab2, *NmbPyr, PyrTab1, PyrTab2,
3040  *NmbPri, PriTab1, PriTab2, *NmbHex, HexTab1, HexTab2,
3041  *BasIdx, *NmbThr ));
3042 }
3043 
3044 int64_t call(lolfreeoctree)(int64_t *OctIdx)
3045 {
3046  return(LolFreeOctree(*OctIdx));
3047 }
3048 
3049 itg call(lolgetboundingbox)( int64_t *OctIdx, itg *typ, itg *MaxItm,
3050  itg *ItmTab, fpn *MinCrd, fpn *MaxCrd, itg *ThrIdx )
3051 {
3052  return(LolGetBoundingBox(*OctIdx, *typ, *MaxItm, ItmTab, MinCrd, MaxCrd, *ThrIdx));
3053 }
3054 
3055 itg call(lolgetnearest)(int64_t *OctIdx, itg *typ, fpn *MinCrd, fpn *MinDis,
3056  fpn *MaxDis, void *UsrPrc, void *UsrDat, itg *ThrIdx)
3057 {
3058  return(LolGetNearest(*OctIdx, *typ, MinCrd, MinDis,
3059  *MaxDis, UsrPrc, UsrDat, *ThrIdx));
3060 }
LolNewOctree
int64_t LolNewOctree(itg NmbVer, fpn *PtrCrd1, fpn *PtrCrd2, itg NmbEdg, itg *PtrEdg1, itg *PtrEdg2, itg NmbTri, itg *PtrTri1, itg *PtrTri2, itg NmbQad, itg *PtrQad1, itg *PtrQad2, itg NmbTet, itg *PtrTet1, itg *PtrTet2, itg NmbPyr, itg *PtrPyr1, itg *PtrPyr2, itg NmbPri, itg *PtrPri1, itg *PtrPri2, itg NmbHex, itg *PtrHex1, itg *PtrHex2, itg BasIdx, itg NmbThr)
Definition: libol1.c:343
MulVec1
static void MulVec1(fpn, fpn *)
BucSct::oct
OctSct * oct
Definition: libol1.c:202
MemSctPtr::adr
void * adr
Definition: libol1.c:163
QadSct::idx
itg idx
Definition: libol1.c:118
lnk
Definition: findLinks.cpp:17
MshSct
Definition: libol1.c:182
EdgIntTri
static itg EdgIntTri(TriSct *, EdgSct *, VerSct *, fpn)
Definition: libol1.c:2340
TreSct::oct
OctSct oct
Definition: libol1.c:221
PrjVerLin
static void PrjVerLin(fpn *, fpn *, fpn *, fpn *)
call
#define call(x)
Definition: libol1.h:104
LnkSctPtr::nex
struct LnkSctPtr * nex
Definition: libol1.c:157
TreSct::bnd
fpn bnd[2][3]
Definition: libol1.c:220
OctSctPtr
Definition: libol1.c:191
TriSct
Definition: libol1.c:104
EdgIntQad
static itg EdgIntQad(HexSct *, itg, EdgSct *, VerSct *, fpn)
Definition: libol1.c:2267
TreSct::NmbFreOct
itg NmbFreOct
Definition: libol1.c:218
GetBucBox
static void GetBucBox(TreSct *, BucSct *, fpn *, fpn *)
MshThrSct::hex
HexSct hex
Definition: libol1.c:176
AddScaVec1
static void AddScaVec1(fpn, fpn *)
lolnewoctree
int64_t call() lolnewoctree(itg *NmbVer, fpn *VerTab1, fpn *VerTab2, itg *NmbEdg, itg *EdgTab1, itg *EdgTab2, itg *NmbTri, itg *TriTab1, itg *TriTab2, itg *NmbQad, itg *QadTab1, itg *QadTab2, itg *NmbTet, itg *TetTab1, itg *TetTab2, itg *NmbPyr, itg *PyrTab1, itg *PyrTab2, itg *NmbPri, itg *PriTab1, itg *PriTab2, itg *NmbHex, itg *HexTab1, itg *HexTab2, itg *BasIdx, itg *NmbThr)
Definition: libol1.c:3027
TreSct::NmbBuc
itg NmbBuc
Definition: libol1.c:218
libol1.h
MshThrSct
Definition: libol1.c:168
LinIntBox
static itg LinIntBox(fpn *, fpn *, fpn *, fpn *, fpn)
Definition: libol1.c:1970
OctThrSct::tag
itg tag
Definition: libol1.c:211
AddVec2
static void AddVec2(fpn *, fpn *)
MinGrdLvl
#define MinGrdLvl
Definition: libol1.c:74
MshThrSct::qad
QadSct qad
Definition: libol1.c:172
TriSct::ani
float ani
Definition: libol1.c:109
TngFlg
#define TngFlg
Definition: libol1.c:77
DisPow
static fpn DisPow(fpn *, fpn *)
SetSonCrd
static void SetSonCrd(itg, fpn *, fpn *, fpn *, fpn *)
GetNrmVec
static fpn GetNrmVec(fpn *)
SetMshBox
static void SetMshBox(TreSct *, MshSct *)
Definition: libol1.c:1438
POW
#define POW(a)
Definition: libol1.c:82
MAX
#define MAX(a, b)
Definition: libol1.c:81
MshThrSct::BakQad
QadSct BakQad
Definition: libol1.c:172
TreSct::CurOctBlk
OctSct * CurOctBlk
Definition: libol1.c:221
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
DisVerTri
static fpn DisVerTri(MshSct *, fpn *, TriSct *)
TreSct::NexMem
MemSct * NexMem
Definition: libol1.c:224
SubOct
static void SubOct(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
GetBucNgb
static BucSct * GetBucNgb(TreSct *, BucSct *, itg)
Definition: libol1.c:1091
QadSct::nrm
fpn nrm[3]
Definition: libol1.c:117
box
Definition: gl2gif.cpp:311
GetTriSrf
static fpn GetTriSrf(TriSct *)
Definition: libol1.c:2628
DisVerOct
static fpn DisVerOct(fpn *, fpn *, fpn *)
AddTet
static void AddTet(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
BucSct::idx
itg idx
Definition: libol1.c:203
EdgIntHex
static itg EdgIntHex(EdgSct *, HexSct *, fpn)
Definition: libol1.c:2033
NewMem
static void * NewMem(TreSct *, size_t)
Definition: libol1.c:2987
GetTriVec
static void GetTriVec(TriSct *, fpn *)
TreSct::MinSiz
fpn MinSiz
Definition: libol1.c:220
TriSct::ver
VerSct * ver[3]
Definition: libol1.c:106
MulVec2
static void MulVec2(fpn, fpn *, fpn *)
DisVerPla
static fpn DisVerPla(fpn *, fpn *, fpn *)
HexSct
Definition: libol1.c:147
IntRayOct
static void IntRayOct(TreSct *, MshSct *, fpn *, fpn *, itg *, fpn *, OctSct *, fpn *, fpn *, itg(void *, itg), void *, itg)
LolGetBoundingBox
itg LolGetBoundingBox(int64_t OctIdx, itg typ, itg MaxItm, itg *ItmTab, fpn MinCrd[3], fpn MaxCrd[3], itg ThrIdx)
Definition: libol1.c:639
LolFreeOctree
size_t LolFreeOctree(int64_t OctIdx)
Definition: libol1.c:623
CUB
#define CUB(a)
Definition: libol1.c:83
VerInsBox
static fpn VerInsBox(fpn *, fpn *, fpn *, fpn)
Definition: libol1.c:2016
AddVer
static void AddVer(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
QadSct::ver
VerSct * ver[4]
Definition: libol1.c:114
MshThrSct::BakTri
TriSct BakTri
Definition: libol1.c:171
OctSctPtr::son
struct OctSctPtr * son
Definition: libol1.c:195
BucSct::pos
char pos[3]
Definition: libol1.c:204
w1
const double w1
Definition: GaussQuadratureHex.cpp:18
TetSct
Definition: libol1.c:122
OctThrSct::ThrTag
itg * ThrTag
Definition: libol1.c:211
VerInsTri
static itg VerInsTri(TriSct *, VerSct *, fpn)
Definition: libol1.c:2411
TreSct::grd
BucSct * grd
Definition: libol1.c:222
LnkItm
static void LnkItm(TreSct *, OctSct *, itg, itg, char)
Definition: libol1.c:1773
ClrVec
static void ClrVec(fpn *)
TreSct::eps
fpn eps
Definition: libol1.c:220
EdgSct::idx
itg idx
Definition: libol1.c:100
lolgetboundingbox
itg call() lolgetboundingbox(int64_t *OctIdx, itg *typ, itg *MaxItm, itg *ItmTab, fpn *MinCrd, fpn *MaxCrd, itg *ThrIdx)
Definition: libol1.c:3049
SubVec3
static void SubVec3(fpn *, fpn *, fpn *)
CrsPrd
static void CrsPrd(fpn *, fpn *, fpn *)
LnkSctPtr
Definition: libol1.c:155
OctSctPtr::sub
unsigned char sub
Definition: libol1.c:197
MshThrSct::edg
EdgSct edg
Definition: libol1.c:170
DisVerTet
static fpn DisVerTet(MshSct *, fpn *, TetSct *)
Definition: libol1.c:2601
OctSctPtr::NmbEdg
unsigned char NmbEdg
Definition: libol1.c:197
TetSct::ver
VerSct * ver[4]
Definition: libol1.c:123
MshSct::BasIdx
itg BasIdx
Definition: libol1.c:186
AniFlg
#define AniFlg
Definition: libol1.c:78
SetEdgTng
static void SetEdgTng(EdgSct *)
Definition: libol1.c:2719
VerInsHex
static itg VerInsHex(VerSct *, HexSct *)
Definition: libol1.c:2248
tvpe
static const itg tvpe[12][2]
Definition: libol1.c:333
SetTriNrm
static void SetTriNrm(TriSct *)
Definition: libol1.c:2707
TetSct::ani
float ani
Definition: libol1.c:127
MaxItmOct
#define MaxItmOct
Definition: libol1.c:72
LolTypVer
@ LolTypVer
Definition: libol1.h:55
MshThrSct::tri
TriSct tri
Definition: libol1.c:171
VerInsOct
static itg VerInsOct(fpn *, fpn *, fpn *)
GetOctLnk
static void GetOctLnk(MshSct *, itg, fpn *, itg *, fpn *, OctSct *, fpn *, fpn *, itg(void *, itg), void *, itg)
AniTri
static void AniTri(MshSct *, itg)
HexSct::edg
EdgSct edg[12]
Definition: libol1.c:149
MshThrSct::tet
TetSct tet
Definition: libol1.c:173
lolgetnearest
itg call() lolgetnearest(int64_t *OctIdx, itg *typ, fpn *MinCrd, fpn *MinDis, fpn *MaxDis, void *UsrPrc, void *UsrDat, itg *ThrIdx)
Definition: libol1.c:3055
TreSct::GrdLvl
itg GrdLvl
Definition: libol1.c:218
TreSct::NmbThr
itg NmbThr
Definition: libol1.c:218
VerSct::crd
fpn crd[3]
Definition: libol1.c:93
TreSct::thr
OctThrSct thr[MaxThr]
Definition: libol1.c:217
TriSct::nrm
fpn nrm[3]
Definition: libol1.c:107
fpn
#define fpn
Definition: libol1.h:76
LolGetNearest
itg LolGetNearest(int64_t OctIdx, itg typ, fpn *VerCrd, fpn *MinDis, fpn MaxDis, itg(UsrPrc)(void *, itg), void *UsrDat, itg ThrIdx)
Definition: libol1.c:786
GetBox
static void GetBox(TreSct *, OctSct *, itg, itg *, itg, itg *, char *, fpn[2][3], fpn, fpn *, fpn *, itg)
MaxThr
#define MaxThr
Definition: libol1.c:79
TreSct::MaxLvl
itg MaxLvl
Definition: libol1.c:218
OctSctPtr::NmbFac
unsigned char NmbFac
Definition: libol1.c:197
LnkSctPtr::typ
itg typ
Definition: libol1.c:156
FreAllMem
static void FreAllMem(TreSct *)
Definition: libol1.c:3008
OctSctPtr::NmbVer
unsigned char NmbVer
Definition: libol1.c:197
MshThrSct::pri
PriSct pri
Definition: libol1.c:175
MemSctPtr
Definition: libol1.c:161
DotPrd
static fpn DotPrd(fpn *, fpn *)
TreSct::BucSiz
fpn BucSiz
Definition: libol1.c:220
EdgSct::ver
VerSct * ver[2]
Definition: libol1.c:98
QadIntHex
static itg QadIntHex(QadSct *, HexSct *, fpn)
Definition: libol1.c:2128
GetPtrItm
static char * GetPtrItm(MshSct *, itg, itg)
Definition: libol1.c:1428
BucSct
Definition: libol1.c:201
VerSct::idx
itg idx
Definition: libol1.c:92
LolTypEdg
@ LolTypEdg
Definition: libol1.h:55
OctThrSct::ThrStk
BucSct ** ThrStk
Definition: libol1.c:212
LolProjectVertex
itg LolProjectVertex(int64_t OctIdx, fpn *VerCrd, itg typ, itg MinItm, fpn *MinCrd, itg ThrIdx)
Definition: libol1.c:930
ItmPerBuc
#define ItmPerBuc
Definition: libol1.c:75
GetVolTet
static fpn GetVolTet(TetSct *)
Definition: libol1.c:2642
VerInsTet
static itg VerInsTet(VerSct *, TetSct *, fpn)
Definition: libol1.c:2208
QadSct::tri
TriSct tri[2]
Definition: libol1.c:116
TreSct::MemUse
size_t MemUse
Definition: libol1.c:219
TriSct::idx
itg idx
Definition: libol1.c:108
itg
#define itg
Definition: libol1.h:69
TetEdg
static const itg TetEdg[6][2]
Definition: libol1.c:330
TetIntHex
static itg TetIntHex(TetSct *, HexSct *, fpn)
Definition: libol1.c:2141
HexSct::qad
QadSct qad[6]
Definition: libol1.c:150
SetItm
static void SetItm(MshSct *, itg, itg, itg, itg)
Definition: libol1.c:1302
TreSct::NmbOct
itg NmbOct
Definition: libol1.c:218
DisVerQad
static fpn DisVerQad(MshSct *, fpn *, QadSct *)
RstMshItm
static void RstMshItm(MshSct *)
Definition: libol1.c:1414
GetCrd
static OctSct * GetCrd(OctSct *, itg, fpn *, fpn *, fpn *)
AddEdg
static void AddEdg(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
MshThrSct::tag
itg tag
Definition: libol1.c:178
MshThrSct::FlgTab
char * FlgTab
Definition: libol1.c:177
LolIntersectSurface
itg LolIntersectSurface(int64_t OctIdx, fpn *VerCrd, fpn *VerTng, fpn *MinDis, fpn MaxDis, itg(UsrPrc)(void *, itg), void *UsrDat, itg ThrIdx)
Definition: libol1.c:857
TetSct::tri
TriSct tri[4]
Definition: libol1.c:125
PyrSct
Definition: libol1.c:131
LnkSctPtr::idx
itg idx
Definition: libol1.c:156
MshThrSct::TagTab
itg * TagTab
Definition: libol1.c:178
OctThrSct::hex
HexSct hex
Definition: libol1.c:210
OctSctPtr::lnk
LnkSct * lnk
Definition: libol1.c:194
MshThrSct::ver
VerSct ver[8]
Definition: libol1.c:169
TreSct
Definition: libol1.c:216
LnkSct
struct LnkSctPtr LnkSct
OctSctPtr::MaxItm
unsigned char MaxItm
Definition: libol1.c:197
TreSct::msh
MshSct * msh
Definition: libol1.c:225
QadSct::edg
EdgSct edg[4]
Definition: libol1.c:115
MaxOctLvl
#define MaxOctLvl
Definition: libol1.c:73
MemSctPtr::nex
struct MemSctPtr * nex
Definition: libol1.c:164
OctThrSct
Definition: libol1.c:208
HexSct::ver
VerSct * ver[8]
Definition: libol1.c:148
MshSct::UsrSiz
size_t UsrSiz[LolNmbTyp]
Definition: libol1.c:184
dis
static fpn dis(fpn *, fpn *)
SetTmpHex
static void SetTmpHex(HexSct *, fpn *, fpn *)
CpyVec
static void CpyVec(fpn *, fpn *)
OctSctPtr::lvl
unsigned char lvl
Definition: libol1.c:197
TetSct::edg
EdgSct edg[6]
Definition: libol1.c:124
EdgSct
Definition: libol1.c:97
MshThrSct::BakTet
TetSct BakTet
Definition: libol1.c:173
BakMshItm
static void BakMshItm(MshSct *)
Definition: libol1.c:1400
LolTypTet
@ LolTypTet
Definition: libol1.h:56
HexSct::idx
itg idx
Definition: libol1.c:151
SubVec2
static void SubVec2(fpn *, fpn *)
MshSct::UsrPtr
char * UsrPtr[LolNmbTyp]
Definition: libol1.c:187
MemSct
struct MemSctPtr MemSct
EdgIntEdg
static itg EdgIntEdg(EdgSct *, EdgSct *, VerSct *, fpn)
Definition: libol1.c:2462
AddTri
static void AddTri(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
MshThrSct::BakVer
VerSct BakVer[8]
Definition: libol1.c:169
MIN
#define MIN(a, b)
Definition: libol1.c:80
LolTypQad
@ LolTypQad
Definition: libol1.h:55
AddQad
static void AddQad(MshSct *, TreSct *, OctSct *, fpn *, fpn *)
PriSct::idx
itg idx
Definition: libol1.c:143
lolfreeoctree
int64_t call() lolfreeoctree(int64_t *OctIdx)
Definition: libol1.c:3044
MshThrSct::BakEdg
EdgSct BakEdg
Definition: libol1.c:170
tvpf
static const itg tvpf[6][4]
Definition: libol1.c:335
TreSct::NexFreLnk
LnkSct * NexFreLnk
Definition: libol1.c:223
TetFacEdg
static const itg TetFacEdg[4][3]
Definition: libol1.c:332
MshThrSct::pyr
PyrSct pyr
Definition: libol1.c:174
OctSct
struct OctSctPtr OctSct
TetSct::idx
itg idx
Definition: libol1.c:126
TriSct::edg
EdgSct edg[3]
Definition: libol1.c:105
BoxIntBox
static itg BoxIntBox(fpn[2][3], fpn[2][3], fpn)
Definition: libol1.c:2783
TetFac
static const itg TetFac[4][3]
Definition: libol1.c:331
AddScaVec2
static void AddScaVec2(fpn, fpn *, fpn *)
MshSct::thr
MshThrSct thr[MaxThr]
Definition: libol1.c:183
VerSct
Definition: libol1.c:91
PyrSct::idx
itg idx
Definition: libol1.c:135
TriIntHex
static itg TriIntHex(TriSct *, HexSct *, fpn)
Definition: libol1.c:2055
PriSct
Definition: libol1.c:139
MshSct::eps
fpn eps
Definition: libol1.c:185
GetTriAni
static fpn GetTriAni(TriSct *)
Definition: libol1.c:2806
LinIntPla
static void LinIntPla(fpn *, fpn *, fpn *, fpn *, fpn *)
Definition: libol1.c:2001
VerInsEdg
static itg VerInsEdg(EdgSct *, VerSct *, fpn)
Definition: libol1.c:2515
OctSctPtr::NmbVol
unsigned char NmbVol
Definition: libol1.c:197
OctSctPtr::ani
unsigned char ani
Definition: libol1.c:197
LolTypTri
@ LolTypTri
Definition: libol1.h:55
LinCmbVec3
static void LinCmbVec3(fpn, fpn *, fpn, fpn *, fpn *)
MemBlkSiz
#define MemBlkSiz
Definition: libol1.c:76
LolNmbTyp
@ LolNmbTyp
Definition: libol1.h:56
DisVerEdg
static fpn DisVerEdg(fpn *, EdgSct *)
MshSct::NmbItm
size_t NmbItm[LolNmbTyp]
Definition: libol1.c:184
MemSctPtr::siz
size_t siz
Definition: libol1.c:162
NrmVec
static void NrmVec(fpn *)
EdgSct::siz
fpn siz
Definition: libol1.c:99
QadSct
Definition: libol1.c:113
OctThrSct::ver
VerSct ver[8]
Definition: libol1.c:209
EdgSct::tng
fpn tng[3]
Definition: libol1.c:99
PrjVerPla
static fpn PrjVerPla(fpn *, fpn *, fpn *, fpn *)