ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
adtLocalSearch.F90
Go to the documentation of this file.
1 module adtlocalsearch
2  !
3  ! Module which contains the subroutines to perform the local
4  ! searches, i.e. the tree traversals.
5  !
6  use constants
7  use adtutils, only: adts, adtleaftype, adtbboxtargettype, stack, &
8  nstack, adtterminate, reallocbboxtargettypeplus, reallocplus, &
9  qsortbboxtargets
10  use adtdata, only: adttype
11  implicit none
12 
13  !=================================================================
14 
15 contains
16 
17  !===============================================================
18 
19  subroutine containmenttreesearch(ADT, coor, &
20  intInfo, uvw, &
21  arrDonor, nCoor, &
22  nInterpol)
23  !
24  ! This routine performs the actual containment search in the
25  ! local tree. It is a local routine in the sense that no
26  ! communication is involved.
27  ! Subroutine intent(in) arguments.
28  ! --------------------------------
29  ! ADT: ADT type whose ADT must be searched
30  ! nCoor: Number of coordinates for which the element must
31  ! be determined.
32  ! coor: The coordinates of these points.
33  ! nInterpol: Number of variables to be interpolated.
34  ! arrDonor: Array with the donor data; needed to obtain the
35  ! interpolated data.
36  ! Subroutine intent(out) arguments.
37  ! ---------------------------------
38  ! intInfo: 2D integer array, in which the following output
39  ! be be stored:
40  ! intInfo(1,:): processor ID of the processor where
41  ! the element is stored. This of course
42  ! is myID. If no element is found this
43  ! value is set to -1.
44  ! intInfo(2,:): The element type of the element.
45  ! intInfo(3,:): The element ID of the element in the
46  ! connectivity.
47  ! uvw: 2D floating point array to store the parametric
48  ! coordinates of the point in the transformed element
49  ! and the interpolated data:
50  ! uvw(1, :): Parametric u-weight.
51  ! uvw(2, :): Parametric v-weight.
52  ! uvw(3, :): Parametric w-weight.
53  ! uvw(4:,:): Interpolated solution, if desired. It is
54  ! possible to call this routine with
55  ! nInterpol == 0.
56  !
57  implicit none
58  !
59  ! Subroutine arguments.
60  !
61  type(adttype), intent(inout) :: ADT
62  integer(kind=intType), intent(in) :: nCoor
63  integer(kind=intType), intent(in) :: nInterpol
64 
65  real(kind=realtype), dimension(:, :), intent(in) :: coor
66  real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
67 
68  integer(kind=intType), dimension(:, :), intent(out) :: intInfo
69  real(kind=realtype), dimension(:, :), intent(out) :: uvw
70  !
71  ! Local variables.
72  !
73  integer(kind=intType), dimension(:), pointer :: BB
74  integer(kind=intType), dimension(:), pointer :: frontLeaves
75  integer(kind=intType), dimension(:), pointer :: frontLeavesNew
76  integer(kind=intType) :: nAllocBB, nAllocFront
77  integer(kind=intType) :: ierr, nn
78  logical :: failed
79 
80  ! Initial allocation of the arrays for the tree traversal.
81 
82  nallocbb = 10
83  nallocfront = 25
84 
85  allocate (bb(nallocbb), frontleaves(nallocfront), &
86  frontleavesnew(nallocfront), stat=ierr)
87  if (ierr /= 0) &
88  call adtterminate(adt, "containmentTreeSearch", &
89  "Memory allocation failure for BB, &
90  &frontLeaves and frontLeavesNew.")
91 
92  ! Loop over the number of coordinates to be treated.
93 
94  coorloop: do nn = 1, ncoor
95 
96  call containmenttreesearchsinglepoint(adt, coor(:, nn), &
97  intinfo(:, nn), uvw(:, nn), arrdonor, ninterpol, bb, &
98  frontleaves, frontleavesnew, failed)
99 
100  end do coorloop
101 
102  ! Release the memory allocated in this routine.
103 
104  deallocate (bb, frontleaves, frontleavesnew, stat=ierr)
105  if (ierr /= 0) &
106  call adtterminate(adt, "containmentTreeSearch", &
107  "Deallocation failure for BB, &
108  &frontLeaves and frontLeavesNew.")
109 
110  end subroutine containmenttreesearch
111 
112  subroutine containmenttreesearchsinglepoint(ADT, coor, &
113  intInfo, uvw, arrDonor, nInterpol, BB, &
114  frontLeaves, frontLeavesNew, failed)
115  !
116  ! This routine is replaces the original containment
117  ! tree search, however, it has been optimized for searching a
118  ! single query point at a time. Specifically, this means that
119  ! there are no allocatable arrays inside this routine; they
120  ! must be supplied exterally. Also since only a single point
121  ! is searched we can fix some of the dimensions. This routine
122  ! requires pointers for BB, frontLeaves and frontLeavesNew to
123  ! be passed to the routine. Since this routine is called a
124  ! very large number of times, these cannot be allocated
125  ! dynamically inside for speed purposes.
126  ! Subroutine intent(in) arguments.
127  ! --------------------------------
128  ! ADT: ADT type whose ADT must be searched
129  ! coor: The coordinate of the point to be searched.
130  ! nInterpol: Number of variables to be interpolated.
131  ! arrDonor: Array with the donor data; needed to obtain the
132  ! interpolated data.
133  ! Subroutine intent(out) arguments.
134  ! ---------------------------------
135  ! intInfo: 1D integer array of length three , in which the
136  ! following output
137  ! be be stored:
138  ! * intInfo(1): processor ID of the processor where
139  ! the element is stored. This of course
140  ! is myID. If no element is found this
141  ! value is set to -1.
142  ! * intInfo(2): The element type of the element.
143  ! * intInfo(3): The element ID of the element in the
144  ! connectivity.
145  ! uvw: 1D floating point array to store the parametric
146  ! coordinates of the point in the transformed element
147  ! and the interpolated data:
148  ! uvw(1): Parametric u-weight.
149  ! uvw(2): Parametric v-weight.
150  ! uvw(3): Parametric w-weight.
151  ! uvw(4:): Interpolated solution(s), if desired. It is
152  ! possible to call this routine with
153  ! nInterpol == 0.
154  !
155  use constants
156  implicit none
157  !
158  ! Subroutine arguments.
159  !
160  type(adttype), intent(inout) :: ADT
161  integer(kind=intType), intent(in) :: nInterpol
162 
163  real(kind=realtype), dimension(3), intent(in) :: coor
164  real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
165 
166  integer(kind=intType), dimension(3), intent(out) :: intInfo
167  real(kind=realtype), dimension(:), intent(out) :: uvw
168  logical, intent(out) :: failed
169  integer(kind=intType), dimension(:), pointer :: BB
170  integer(kind=intType), dimension(:), pointer :: frontLeaves
171  integer(kind=intType), dimension(:), pointer :: frontLeavesNew
172 
173  !
174  ! Local parameters used in the Newton algorithm.
175  !
176  integer(kind=intType), parameter :: iterMax = 15
177  real(kind=realtype), parameter :: adteps = 1.e-25_realtype
178  real(kind=realtype), parameter :: thresconv = 1.e-10_realtype
179  !
180  ! Local variables.
181  !
182  integer :: ierr
183 
184  integer(kind=intType) :: ii, kk, ll, mm, nn
185  integer(kind=intType) :: nBB, nFrontLeaves, nFrontLeavesNew
186  integer(kind=intType) :: nAllocBB, nAllocFront
187  integer(kind=intType) :: i, nNodeElement
188 
189  integer(kind=intType), dimension(8) :: n
190 
191  real(kind=realtype) :: u, v, w, uv, uw, vw, wvu, du, dv, dw
192  real(kind=realtype) :: oneminusu, oneminusv, oneminusw
193  real(kind=realtype) :: oneminusuminusv
194  real(kind=realtype) :: a11, a12, a13, a21, a22, a23
195  real(kind=realtype) :: a31, a32, a33, val
196 
197  real(kind=realtype), dimension(3) :: x, f
198  real(kind=realtype), dimension(8) :: weight
199  real(kind=realtype), dimension(3, 2:8) :: xn
200 
201  real(kind=realtype), dimension(:, :), pointer :: xbbox
202 
203  logical :: elementFound
204 
205  type(adtleaftype), dimension(:), pointer :: ADTree
206 
207  ! Set some pointers to make the code more readable.
208 
209  xbbox => adt%xBBox
210  adtree => adt%ADTree
211 
212  ! Determine the sizes from the arrays we have been passed
213 
214  nallocbb = size(bb)
215  nallocfront = size(frontleaves)
216 
217  ! Initialize the processor ID to -1 to indicate that no
218  ! corresponding volume element is found.
219 
220  intinfo(1) = -1
221  failed = .true.
222  !
223  ! Part 1. Traverse the tree and determine the target
224  ! bounding boxes, which may contain the element.
225  !
226  ! Start at the root, i.e. set the front leaf to the root leaf.
227  ! Also initialize the number of possible bounding boxes to 0.
228 
229  nbb = 0
230 
231  nfrontleaves = 1
232  frontleaves(1) = 1
233 
234  treetraversalloop: do
235 
236  ! Initialize the number of leaves for the new front, i.e.
237  ! the front of the next round, to 0.
238 
239  nfrontleavesnew = 0
240 
241  ! Loop over the leaves of the current front.
242 
243  currentfrontloop: do ii = 1, nfrontleaves
244 
245  ! Store the ID of the leaf a bit easier and loop over
246  ! its two children.
247 
248  ll = frontleaves(ii)
249 
250  childrenloop: do mm = 1, 2
251 
252  ! Determine whether this child contains a bounding box
253  ! or a leaf of the next level.
254 
255  kk = adtree(ll)%children(mm)
256 
257  terminaltest: if (kk < 0) then
258 
259  ! Child contains a bounding box. Check if the
260  ! coordinate is inside the bounding box.
261 
262  kk = -kk
263  if (coor(1) >= xbbox(1, kk) .and. &
264  coor(1) <= xbbox(4, kk) .and. &
265  coor(2) >= xbbox(2, kk) .and. &
266  coor(2) <= xbbox(5, kk) .and. &
267  coor(3) >= xbbox(3, kk) .and. &
268  coor(3) <= xbbox(6, kk)) then
269 
270  ! Coordinate is inside the bounding box. Store the
271  ! bounding box in the list of possible candidates.
272 
273  if (nbb == nallocbb) &
274  call reallocplus(bb, nallocbb, 100, adt)
275 
276  nbb = nbb + 1
277  bb(nbb) = kk
278  end if
279 
280  else terminaltest
281 
282  ! Child contains a leaf. Check if the coordinate is
283  ! inside the bounding box of the leaf.
284 
285  if (coor(1) >= adtree(kk)%xMin(1) .and. &
286  coor(1) <= adtree(kk)%xMax(4) .and. &
287  coor(2) >= adtree(kk)%xMin(2) .and. &
288  coor(2) <= adtree(kk)%xMax(5) .and. &
289  coor(3) >= adtree(kk)%xMin(3) .and. &
290  coor(3) <= adtree(kk)%xMax(6)) then
291 
292  ! Coordinate is inside the leaf. Store the leaf in
293  ! the list for the new front.
294 
295  if (nfrontleavesnew == nallocfront) then
296  i = nallocfront
297  call reallocplus(frontleavesnew, i, 250, adt)
298  call reallocplus(frontleaves, nallocfront, 250, adt)
299  end if
300 
301  nfrontleavesnew = nfrontleavesnew + 1
302  frontleavesnew(nfrontleavesnew) = kk
303 
304  end if
305 
306  end if terminaltest
307 
308  end do childrenloop
309 
310  end do currentfrontloop
311 
312  ! End of the loop over the current front. If the new front
313  ! is empty the entire tree has been traversed and an exit is
314  ! made from the corresponding loop.
315 
316  if (nfrontleavesnew == 0) exit treetraversalloop
317 
318  ! Copy the data of the new front leaves into the current
319  ! front for the next round.
320 
321  nfrontleaves = nfrontleavesnew
322  do ll = 1, nfrontleaves
323  frontleaves(ll) = frontleavesnew(ll)
324  end do
325 
326  end do treetraversalloop
327  !
328  ! Part 2: Loop over the selected bounding boxes and check if
329  ! the corresponding elements contain the point.
330  !
331  elementfound = .false.
332 
333  bboxloop: do mm = 1, nbb
334 
335  ! Determine the element type stored in this bounding box.
336 
337  kk = bb(mm)
338  select case (adt%elementType(kk))
339 
340  case (adttetrahedron)
341 
342  ! Element is a tetrahedron.
343  ! Compute the coordinates relative to node 1.
344 
345  ll = adt%elementID(kk)
346  n(1) = adt%tetraConn(1, ll)
347 
348  do i = 2, 4
349  n(i) = adt%tetraConn(i, ll)
350 
351  xn(1, i) = adt%coor(1, n(i)) - adt%coor(1, n(1))
352  xn(2, i) = adt%coor(2, n(i)) - adt%coor(2, n(1))
353  xn(3, i) = adt%coor(3, n(i)) - adt%coor(3, n(1))
354  end do
355 
356  x(1) = coor(1) - adt%coor(1, n(1))
357  x(2) = coor(2) - adt%coor(2, n(1))
358  x(3) = coor(3) - adt%coor(3, n(1))
359 
360  ! Determine the matrix for the linear transformation
361  ! from the standard element to the current element.
362 
363  a11 = xn(1, 2); a12 = xn(1, 3); a13 = xn(1, 4)
364  a21 = xn(2, 2); a22 = xn(2, 3); a23 = xn(2, 4)
365  a31 = xn(3, 2); a32 = xn(3, 3); a33 = xn(3, 4)
366 
367  ! Compute the determinant. Make sure that it is not zero
368  ! and invert the value.
369 
370  val = a11 * (a22 * a33 - a32 * a23) + a21 * (a13 * a32 - a12 * a33) &
371  + a31 * (a12 * a23 - a13 * a22)
372  val = sign(one, val) / max(abs(val), adteps)
373 
374  ! Compute the u, v, w weights for the given coordinate.
375 
376  u = val * ((a22 * a33 - a23 * a32) * x(1) &
377  + (a13 * a32 - a12 * a33) * x(2) &
378  + (a12 * a23 - a13 * a22) * x(3))
379  v = val * ((a23 * a31 - a21 * a33) * x(1) &
380  + (a11 * a33 - a13 * a31) * x(2) &
381  + (a13 * a21 - a11 * a23) * x(3))
382  w = val * ((a21 * a32 - a22 * a31) * x(1) &
383  + (a12 * a31 - a11 * a32) * x(2) &
384  + (a11 * a22 - a12 * a21) * x(3))
385 
386  ! Check if the coordinate is inside the tetrahedron.
387  ! If so, set elementFound to .true. and determine the
388  ! interpolation weights.
389 
390  if (u >= zero .and. v >= zero .and. &
391  w >= zero .and. (u + v + w) <= one) then
392  elementfound = .true.
393  failed = .false. !No iteration for tetrahedrons
394 
395  ! Set the number of interpolation nodes to 4 and
396  ! determine the interpolation weights.
397 
398  nnodeelement = 4
399 
400  weight(1) = one - u - v - w
401  weight(2) = u
402  weight(3) = v
403  weight(4) = w
404  end if
405 
406  !=========================================================
407 
408  case (adtpyramid)
409 
410  ! Element is a pyramid.
411  ! Compute the coordinates relative to node 1.
412 
413  ll = adt%elementID(kk)
414  n(1) = adt%pyraConn(1, ll)
415 
416  do i = 2, 5
417  n(i) = adt%pyraConn(i, ll)
418 
419  xn(1, i) = adt%coor(1, n(i)) - adt%coor(1, n(1))
420  xn(2, i) = adt%coor(2, n(i)) - adt%coor(2, n(1))
421  xn(3, i) = adt%coor(3, n(i)) - adt%coor(3, n(1))
422  end do
423 
424  x(1) = coor(1) - adt%coor(1, n(1))
425  x(2) = coor(2) - adt%coor(2, n(1))
426  x(3) = coor(3) - adt%coor(3, n(1))
427 
428  ! Modify the coordinates of node 3, such that it
429  ! corresponds to the weights of the u*v term in the
430  ! transformation.
431 
432  xn(1, 3) = xn(1, 3) - xn(1, 2) - xn(1, 4)
433  xn(2, 3) = xn(2, 3) - xn(2, 2) - xn(2, 4)
434  xn(3, 3) = xn(3, 3) - xn(3, 2) - xn(3, 4)
435 
436  ! Set the starting values of u, v and w such that it is
437  ! somewhere in the middle of the element. In this way the
438  ! Jacobian matrix is always regular, even if the element
439  ! is degenerate.
440 
441  u = half; v = half; w = half
442 
443  ! The Newton algorithm to determine the parametric
444  ! weights u, v and w for the given coordinate.
445 
446  newtonpyra: do ll = 1, itermax
447 
448  ! Compute the RHS.
449 
450  uv = u * v
451  oneminusw = one - w
452 
453  f(1) = oneminusw * (xn(1, 2) * u + xn(1, 4) * v + xn(1, 3) * uv) &
454  + xn(1, 5) * w - x(1)
455  f(2) = oneminusw * (xn(2, 2) * u + xn(2, 4) * v + xn(2, 3) * uv) &
456  + xn(2, 5) * w - x(2)
457  f(3) = oneminusw * (xn(3, 2) * u + xn(3, 4) * v + xn(3, 3) * uv) &
458  + xn(3, 5) * w - x(3)
459 
460  ! Compute the Jacobian.
461 
462  a11 = oneminusw * (xn(1, 2) + xn(1, 3) * v)
463  a12 = oneminusw * (xn(1, 4) + xn(1, 3) * u)
464  a13 = xn(1, 5) - xn(1, 2) * u - xn(1, 4) * v - xn(1, 3) * uv
465 
466  a21 = oneminusw * (xn(2, 2) + xn(2, 3) * v)
467  a22 = oneminusw * (xn(2, 4) + xn(2, 3) * u)
468  a23 = xn(2, 5) - xn(2, 2) * u - xn(2, 4) * v - xn(2, 3) * uv
469 
470  a31 = oneminusw * (xn(3, 2) + xn(3, 3) * v)
471  a32 = oneminusw * (xn(3, 4) + xn(3, 3) * u)
472  a33 = xn(3, 5) - xn(3, 2) * u - xn(3, 4) * v - xn(3, 3) * uv
473 
474  ! Compute the determinant. Make sure that it is not zero
475  ! and invert the value. The cut off is needed to be able
476  ! to handle exceptional cases for degenerate elements.
477 
478  val = a11 * (a22 * a33 - a32 * a23) + a21 * (a13 * a32 - a12 * a33) &
479  + a31 * (a12 * a23 - a13 * a22)
480  val = sign(one, val) / max(abs(val), adteps)
481 
482  ! Compute the new values of u, v and w.
483 
484  du = val * ((a22 * a33 - a23 * a32) * f(1) &
485  + (a13 * a32 - a12 * a33) * f(2) &
486  + (a12 * a23 - a13 * a22) * f(3))
487  dv = val * ((a23 * a31 - a21 * a33) * f(1) &
488  + (a11 * a33 - a13 * a31) * f(2) &
489  + (a13 * a21 - a11 * a23) * f(3))
490  dw = val * ((a21 * a32 - a22 * a31) * f(1) &
491  + (a12 * a31 - a11 * a32) * f(2) &
492  + (a11 * a22 - a12 * a21) * f(3))
493 
494  u = u - du; v = v - dv; w = w - dw
495 
496  ! Exit the loop if the update of the parametric
497  ! weights is below the threshold
498 
499  val = sqrt(du * du + dv * dv + dw * dw)
500  if (val <= thresconv) then
501  failed = .false.
502  exit newtonpyra
503  end if
504  end do newtonpyra
505 
506  ! Check if the coordinate is inside the pyramid.
507  ! If so, set elementFound to .true. and determine the
508  ! interpolation weights.
509 
510  if (u >= zero .and. v >= zero .and. &
511  w >= zero .and. (u + w) <= one .and. &
512  (v + w) <= one .and. .not. failed) then
513  elementfound = .true.
514 
515  ! Set the number of interpolation nodes to 5 and
516  ! determine the interpolation weights.
517 
518  nnodeelement = 5
519 
520  oneminusu = one - u
521  oneminusv = one - v
522  oneminusw = one - w
523 
524  weight(1) = oneminusu * oneminusv * oneminusw
525  weight(2) = u * oneminusv * oneminusw
526  weight(3) = u * v * oneminusw
527  weight(4) = oneminusu * v * oneminusw
528  weight(5) = w
529  end if
530 
531  !=========================================================
532 
533  case (adtprism)
534 
535  ! Element is a prism.
536  ! Compute the coordinates relative to node 1.
537 
538  ll = adt%elementID(kk)
539  n(1) = adt%prismsConn(1, ll)
540 
541  do i = 2, 6
542  n(i) = adt%prismsConn(i, ll)
543 
544  xn(1, i) = adt%coor(1, n(i)) - adt%coor(1, n(1))
545  xn(2, i) = adt%coor(2, n(i)) - adt%coor(2, n(1))
546  xn(3, i) = adt%coor(3, n(i)) - adt%coor(3, n(1))
547  end do
548 
549  x(1) = coor(1) - adt%coor(1, n(1))
550  x(2) = coor(2) - adt%coor(2, n(1))
551  x(3) = coor(3) - adt%coor(3, n(1))
552 
553  ! Modify the coordinates of node 5 and 6, such that they
554  ! correspond to the weights of the u*w and v*w term in the
555  ! transformation respectively.
556 
557  xn(1, 5) = xn(1, 5) - xn(1, 2) - xn(1, 4)
558  xn(2, 5) = xn(2, 5) - xn(2, 2) - xn(2, 4)
559  xn(3, 5) = xn(3, 5) - xn(3, 2) - xn(3, 4)
560 
561  xn(1, 6) = xn(1, 6) - xn(1, 3) - xn(1, 4)
562  xn(2, 6) = xn(2, 6) - xn(2, 3) - xn(2, 4)
563  xn(3, 6) = xn(3, 6) - xn(3, 3) - xn(3, 4)
564 
565  ! Set the starting values of u, v and w such that it is
566  ! somewhere in the middle of the element. In this way the
567  ! Jacobian matrix is always regular, even if the element
568  ! is degenerate.
569 
570  u = fourth; v = fourth; w = half
571 
572  ! The Newton algorithm to determine the parametric
573  ! weights u, v and w for the given coordinate.
574 
575  newtonprisms: do ll = 1, itermax
576 
577  ! Compute the RHS.
578 
579  uw = u * w; vw = v * w
580 
581  f(1) = xn(1, 2) * u + xn(1, 3) * v + xn(1, 4) * w &
582  + xn(1, 5) * uw + xn(1, 6) * vw - x(1)
583  f(2) = xn(2, 2) * u + xn(2, 3) * v + xn(2, 4) * w &
584  + xn(2, 5) * uw + xn(2, 6) * vw - x(2)
585  f(3) = xn(3, 2) * u + xn(3, 3) * v + xn(3, 4) * w &
586  + xn(3, 5) * uw + xn(3, 6) * vw - x(3)
587 
588  ! Compute the Jacobian.
589 
590  a11 = xn(1, 2) + xn(1, 5) * w
591  a12 = xn(1, 3) + xn(1, 6) * w
592  a13 = xn(1, 4) + xn(1, 5) * u + xn(1, 6) * v
593 
594  a21 = xn(2, 2) + xn(2, 5) * w
595  a22 = xn(2, 3) + xn(2, 6) * w
596  a23 = xn(2, 4) + xn(2, 5) * u + xn(2, 6) * v
597 
598  a31 = xn(3, 2) + xn(3, 5) * w
599  a32 = xn(3, 3) + xn(3, 6) * w
600  a33 = xn(3, 4) + xn(3, 5) * u + xn(3, 6) * v
601 
602  ! Compute the determinant. Make sure that it is not zero
603  ! and invert the value. The cut off is needed to be able
604  ! to handle exceptional cases for degenerate elements.
605 
606  val = a11 * (a22 * a33 - a32 * a23) + a21 * (a13 * a32 - a12 * a33) &
607  + a31 * (a12 * a23 - a13 * a22)
608  val = sign(one, val) / max(abs(val), adteps)
609 
610  ! Compute the new values of u, v and w.
611 
612  du = val * ((a22 * a33 - a23 * a32) * f(1) &
613  + (a13 * a32 - a12 * a33) * f(2) &
614  + (a12 * a23 - a13 * a22) * f(3))
615  dv = val * ((a23 * a31 - a21 * a33) * f(1) &
616  + (a11 * a33 - a13 * a31) * f(2) &
617  + (a13 * a21 - a11 * a23) * f(3))
618  dw = val * ((a21 * a32 - a22 * a31) * f(1) &
619  + (a12 * a31 - a11 * a32) * f(2) &
620  + (a11 * a22 - a12 * a21) * f(3))
621 
622  u = u - du; v = v - dv; w = w - dw
623 
624  ! Exit the loop if the update of the parametric
625  ! weights is below the threshold
626 
627  val = sqrt(du * du + dv * dv + dw * dw)
628  if (val <= thresconv) then
629  failed = .false.
630  exit newtonprisms
631  end if
632 
633  end do newtonprisms
634 
635  ! Check if the coordinate is inside the prism.
636  ! If so, set elementFound to .true. and determine the
637  ! interpolation weights.
638 
639  if (u >= zero .and. v >= zero .and. &
640  w >= zero .and. w <= one .and. &
641  (u + v) <= one .and. .not. failed) then
642  elementfound = .true.
643 
644  ! Set the number of interpolation nodes to 6 and
645  ! determine the interpolation weights.
646 
647  nnodeelement = 6
648 
649  oneminusuminusv = one - u - v
650  oneminusw = one - w
651 
652  weight(1) = oneminusuminusv * oneminusw
653  weight(2) = u * oneminusw
654  weight(3) = v * oneminusw
655  weight(4) = oneminusuminusv * w
656  weight(5) = u * w
657  weight(6) = v * w
658  end if
659 
660  !=========================================================
661 
662  case (adthexahedron)
663 
664  ! Element is a hexahedron.
665  ! Compute the coordinates relative to node 1.
666 
667  ll = adt%elementID(kk)
668  n(1) = adt%hexaConn(1, ll)
669 
670  do i = 2, 8
671  n(i) = adt%hexaConn(i, ll)
672 
673  xn(1, i) = adt%coor(1, n(i)) - adt%coor(1, n(1))
674  xn(2, i) = adt%coor(2, n(i)) - adt%coor(2, n(1))
675  xn(3, i) = adt%coor(3, n(i)) - adt%coor(3, n(1))
676  end do
677 
678  x(1) = coor(1) - adt%coor(1, n(1))
679  x(2) = coor(2) - adt%coor(2, n(1))
680  x(3) = coor(3) - adt%coor(3, n(1))
681 
682  ! Modify the coordinates of node 3, 6, 8 and 7 such that
683  ! they correspond to the weights of the u*v, u*w, v*w and
684  ! u*v*w term in the transformation respectively.
685 
686  xn(1, 7) = xn(1, 7) + xn(1, 2) + xn(1, 4) + xn(1, 5) &
687  - xn(1, 3) - xn(1, 6) - xn(1, 8)
688  xn(2, 7) = xn(2, 7) + xn(2, 2) + xn(2, 4) + xn(2, 5) &
689  - xn(2, 3) - xn(2, 6) - xn(2, 8)
690  xn(3, 7) = xn(3, 7) + xn(3, 2) + xn(3, 4) + xn(3, 5) &
691  - xn(3, 3) - xn(3, 6) - xn(3, 8)
692 
693  xn(1, 3) = xn(1, 3) - xn(1, 2) - xn(1, 4)
694  xn(2, 3) = xn(2, 3) - xn(2, 2) - xn(2, 4)
695  xn(3, 3) = xn(3, 3) - xn(3, 2) - xn(3, 4)
696 
697  xn(1, 6) = xn(1, 6) - xn(1, 2) - xn(1, 5)
698  xn(2, 6) = xn(2, 6) - xn(2, 2) - xn(2, 5)
699  xn(3, 6) = xn(3, 6) - xn(3, 2) - xn(3, 5)
700 
701  xn(1, 8) = xn(1, 8) - xn(1, 4) - xn(1, 5)
702  xn(2, 8) = xn(2, 8) - xn(2, 4) - xn(2, 5)
703  xn(3, 8) = xn(3, 8) - xn(3, 4) - xn(3, 5)
704 
705  ! Set the starting values of u, v and w such that it is
706  ! somewhere in the middle of the element. In this way the
707  ! Jacobian matrix is always regular, even if the element
708  ! is degenerate.
709 
710  u = half; v = half; w = half
711 
712  ! The Newton algorithm to determine the parametric
713  ! weights u, v and w for the given coordinate.
714 
715  newtonhexa: do ll = 1, itermax
716 
717  ! Compute the RHS.
718 
719  uv = u * v; uw = u * w; vw = v * w; wvu = u * v * w
720 
721  f(1) = xn(1, 2) * u + xn(1, 4) * v + xn(1, 5) * w &
722  + xn(1, 3) * uv + xn(1, 6) * uw + xn(1, 8) * vw &
723  + xn(1, 7) * wvu - x(1)
724  f(2) = xn(2, 2) * u + xn(2, 4) * v + xn(2, 5) * w &
725  + xn(2, 3) * uv + xn(2, 6) * uw + xn(2, 8) * vw &
726  + xn(2, 7) * wvu - x(2)
727  f(3) = xn(3, 2) * u + xn(3, 4) * v + xn(3, 5) * w &
728  + xn(3, 3) * uv + xn(3, 6) * uw + xn(3, 8) * vw &
729  + xn(3, 7) * wvu - x(3)
730 
731  ! Compute the Jacobian.
732 
733  a11 = xn(1, 2) + xn(1, 3) * v + xn(1, 6) * w + xn(1, 7) * vw
734  a12 = xn(1, 4) + xn(1, 3) * u + xn(1, 8) * w + xn(1, 7) * uw
735  a13 = xn(1, 5) + xn(1, 6) * u + xn(1, 8) * v + xn(1, 7) * uv
736 
737  a21 = xn(2, 2) + xn(2, 3) * v + xn(2, 6) * w + xn(2, 7) * vw
738  a22 = xn(2, 4) + xn(2, 3) * u + xn(2, 8) * w + xn(2, 7) * uw
739  a23 = xn(2, 5) + xn(2, 6) * u + xn(2, 8) * v + xn(2, 7) * uv
740 
741  a31 = xn(3, 2) + xn(3, 3) * v + xn(3, 6) * w + xn(3, 7) * vw
742  a32 = xn(3, 4) + xn(3, 3) * u + xn(3, 8) * w + xn(3, 7) * uw
743  a33 = xn(3, 5) + xn(3, 6) * u + xn(3, 8) * v + xn(3, 7) * uv
744 
745  ! Compute the determinant. Make sure that it is not zero
746  ! and invert the value. The cut off is needed to be able
747  ! to handle exceptional cases for degenerate elements.
748 
749  val = a11 * (a22 * a33 - a32 * a23) + a21 * (a13 * a32 - a12 * a33) &
750  + a31 * (a12 * a23 - a13 * a22)
751  val = sign(one, val) / max(abs(val), adteps)
752 
753  ! Compute the new values of u, v and w.
754 
755  du = val * ((a22 * a33 - a23 * a32) * f(1) &
756  + (a13 * a32 - a12 * a33) * f(2) &
757  + (a12 * a23 - a13 * a22) * f(3))
758  dv = val * ((a23 * a31 - a21 * a33) * f(1) &
759  + (a11 * a33 - a13 * a31) * f(2) &
760  + (a13 * a21 - a11 * a23) * f(3))
761  dw = val * ((a21 * a32 - a22 * a31) * f(1) &
762  + (a12 * a31 - a11 * a32) * f(2) &
763  + (a11 * a22 - a12 * a21) * f(3))
764 
765  u = u - du; v = v - dv; w = w - dw
766 
767  ! Exit the loop if the update of the parametric
768  ! weights is below the threshold
769 
770  val = sqrt(du * du + dv * dv + dw * dw)
771  if (val <= thresconv) then
772  failed = .false.
773  exit newtonhexa
774  end if
775 
776  end do newtonhexa
777 
778  ! Check if the coordinate is inside the hexahedron.
779  ! If so, set elementFound to .true. and determine the
780  ! interpolation weights.
781 
782  if (u >= zero .and. u <= one .and. &
783  v >= zero .and. v <= one .and. &
784  w >= zero .and. w <= one .and. .not. failed) then
785  elementfound = .true.
786 
787  ! Set the number of interpolation nodes to 8 and
788  ! determine the interpolation weights.
789 
790  nnodeelement = 8
791 
792  oneminusu = one - u
793  oneminusv = one - v
794  oneminusw = one - w
795 
796  weight(1) = oneminusu * oneminusv * oneminusw
797  weight(2) = u * oneminusv * oneminusw
798  weight(3) = u * v * oneminusw
799  weight(4) = oneminusu * v * oneminusw
800  weight(5) = oneminusu * oneminusv * w
801  weight(6) = u * oneminusv * w
802  weight(7) = u * v * w
803  weight(8) = oneminusu * v * w
804  end if
805 
806  end select
807 
808  ! If the coordinate is inside the element store all the
809  ! necessary information and exit the loop over the target
810  ! bounding boxes.
811 
812  if (elementfound) then
813 
814  ! The processor, element type and local element ID.
815 
816  intinfo(1) = adt%myID
817  intinfo(2) = adt%elementType(kk)
818  intinfo(3) = adt%elementID(kk)
819 
820  ! The parametric weights.
821 
822  uvw(1) = u
823  uvw(2) = v
824  uvw(3) = w
825 
826  ! The interpolated solution.
827 
828  do ll = 1, ninterpol
829  ii = 3 + ll
830  uvw(ii) = weight(1) * arrdonor(ll, n(1))
831  do i = 2, nnodeelement
832  uvw(ii) = uvw(ii) + weight(i) * arrdonor(ll, n(i))
833  end do
834  end do
835 
836  ! And exit the loop over the bounding boxes.
837 
838  exit bboxloop
839  end if
840 
841  end do bboxloop
842 
843  end subroutine containmenttreesearchsinglepoint
844 
845  subroutine mindistancetreesearch(ADT, coor, &
846  intInfo, uvw, &
847  arrDonor, nCoor, &
848  nInterpol)
849  !
850  ! This routine performs the actual minimum distance search in
851  ! the local tree. It is a local routine in the sense that no
852  ! communication is involved.
853  ! Subroutine intent(in) arguments.
854  ! --------------------------------
855  ! ADT: ADT type whose ADT must be searched
856  ! nCoor: Number of coordinates for which the element must
857  ! be determined.
858  ! coor: The coordinates and the currently stored minimum
859  ! distance squared of these points:
860  ! coor(1,;): Coordinate 1.
861  ! coor(2,;): Coordinate 2.
862  ! coor(3,;): Coordinate 3.
863  ! coor(4,;): The currently stored minimum distance
864  ! squared.
865  ! nInterpol: Number of variables to be interpolated.
866  ! arrDonor: Array with the donor data; needed to obtain the
867  ! interpolated data.
868  ! Subroutine intent(out) arguments.
869  ! ---------------------------------
870  ! intInfo: 2D integer array, in which the following output
871  ! will be stored:
872  ! intInfo(1,:): processor ID of the processor where
873  ! the element is stored. This of course
874  ! is myID. If no element is found this
875  ! value is set to -1.
876  ! intInfo(2,:): The element type of the element.
877  ! intInfo(3,:): The element ID of the element in the
878  ! the connectivity.
879  ! uvw: 2D floating point array to store the parametric
880  ! coordinates of the point in the transformed element
881  ! as well as the new distance squared and the
882  ! interpolated solution:
883  ! uvw(1, :): Parametric u-weight.
884  ! uvw(2, :): Parametric v-weight.
885  ! uvw(3, :): Parametric w-weight.
886  ! uvw(4, :): The new distance squared.
887  ! uvw(5:,:): Interpolated solution, if desired. It is
888  ! possible to call this routine with
889  ! nInterpol == 0.
890  !
891  implicit none
892  !
893  ! Subroutine arguments.
894  !
895  type(adttype), intent(inout) :: ADT
896  integer(kind=intType), intent(in) :: nCoor
897  integer(kind=intType), intent(in) :: nInterpol
898 
899  real(kind=realtype), dimension(:, :), intent(in) :: coor
900  real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
901 
902  integer(kind=intType), dimension(:, :), intent(out) :: intInfo
903  real(kind=realtype), dimension(:, :), intent(out) :: uvw
904  !
905  ! Local parameters used in the Newton algorithm.
906  !
907  integer(kind=intType), parameter :: iterMax = 15
908  real(kind=realtype), parameter :: adteps = 1.e-25_realtype
909  real(kind=realtype), parameter :: thresconv = 1.e-10_realtype
910  !
911  ! Local variables.
912  !
913  integer :: ierr, nn
914  integer(kind=intType) :: nAllocBB, nAllocFront, nStack
915  integer(kind=intType), dimension(:), pointer :: frontLeaves
916  integer(kind=intType), dimension(:), pointer :: frontLeavesNew
917  type(adtbboxtargettype), dimension(:), pointer :: BB
918 
919  ! Initial allocation of the arrays for the tree traversal as well
920  ! as the stack array used in the qsort routine. The latter is
921  ! done, because the qsort routine is called for every coordinate
922  ! and therefore it is more efficient to allocate the stack once
923  ! rather than over and over again. The disadvantage of course is
924  ! that an essentially local variable, stack, is now stored in
925  ! adtData.
926 
927  nallocbb = 10
928  nallocfront = 25
929  nstack = 100
930 
931  allocate (stack(nstack), bb(nallocbb), frontleaves(nallocfront), &
932  frontleavesnew(nallocfront), stat=ierr)
933  if (ierr /= 0) &
934  call adtterminate(adt, "minDistanceTreeSearch", &
935  "Memory allocation failure for stack, BB, &
936  &etc.")
937 
938  ! Loop over the number of coordinates to be treated.
939 
940  coorloop: do nn = 1, ncoor
941 
942  call mindistancetreesearchsinglepoint(adt, coor(:, nn), &
943  intinfo(:, nn), uvw(:, nn), arrdonor, ninterpol, bb, &
944  frontleaves, frontleavesnew)
945 
946  end do coorloop
947 
948  ! Release the memory allocated in this routine.
949 
950  deallocate (stack, bb, frontleaves, frontleavesnew, stat=ierr)
951  if (ierr /= 0) &
952  call adtterminate(adt, "minDistanceTreeSearch", &
953  "Deallocation failure for stack, BB, etc.")
954 
955  end subroutine mindistancetreesearch
956 
957  subroutine mindistancetreesearchsinglepoint(ADT, coor, intInfo, &
958  uvw, arrDonor, nInterpol, BB, frontLeaves, frontLeavesNew)
959  !
960  ! This routine performs the actual minimum distance search for
961  ! a single point on the local tree. It is local in the sens
962  ! that no communication is involved. This routine does the
963  ! actual search. The minDistanceTreeSearch is just a wrapper
964  ! around this routine. The reason for the split is that the
965  ! overset mesh connectivity requires efficient calling with
966  ! a single coordinate. Therefore, this rouine does not
967  ! allocate/deallocate any variables.
968  ! Subroutine intent(in) arguments.
969  ! --------------------------------
970  ! ADT: ADT type whose ADT must be searched
971  ! coor: The coordinates and the currently stored minimum
972  ! distance squared of these points:
973  ! coor(1): Coordinate 1.
974  ! coor(2): Coordinate 2.
975  ! coor(3): Coordinate 3.
976  ! coor(4): The currently stored minimum distance
977  ! squared.
978  ! nInterpol: Number of variables to be interpolated.
979  ! arrDonor: Array with the donor data; needed to obtain the
980  ! interpolated data.
981  ! Subroutine intent(out) arguments.
982  ! ---------------------------------
983  ! intInfo: 1D integer array, in which the following output
984  ! will be stored:
985  ! intInfo(1): processor ID of the processor where
986  ! the element is stored. This of course
987  ! is myID. If no element is found this
988  ! value is set to -1.
989  ! intInfo(2): The element type of the element.
990  ! intInfo(3): The element ID of the element in the
991  ! the connectivity.
992  ! uvw: 2D floating point array to store the parametric
993  ! coordinates of the point in the transformed element
994  ! as well as the new distance squared and the
995  ! interpolated solution:
996  ! uvw(1): Parametric u-weight.
997  ! uvw(2): Parametric v-weight.
998  ! uvw(3): Parametric w-weight.
999  ! uvw(4): The new distance squared.
1000  ! uvw(5): Interpolated solution, if desired. It is
1001  ! possible to call this routine with
1002  ! nInterpol == 0.
1003  !
1004  implicit none
1005  !
1006  ! Subroutine arguments.
1007  !
1008  type(adttype), intent(inout) :: ADT
1009  integer(kind=intType), intent(in) :: nInterpol
1010 
1011  real(kind=realtype), dimension(4), intent(in) :: coor
1012  real(kind=realtype), dimension(:, :), intent(in) :: arrdonor
1013 
1014  integer(kind=intType), dimension(3), intent(out) :: intInfo
1015  real(kind=realtype), dimension(5), intent(out) :: uvw
1016  integer(kind=intType), dimension(:), pointer :: frontLeaves
1017  integer(kind=intType), dimension(:), pointer :: frontLeavesNew
1018  type(adtbboxtargettype), dimension(:), pointer :: BB
1019  !
1020  ! Local parameters used in the Newton algorithm.
1021  !
1022  integer(kind=intType), parameter :: iterMax = 15
1023  real(kind=realtype), parameter :: adteps = 1.e-25_realtype
1024  real(kind=realtype), parameter :: thresconv = 1.e-10_realtype
1025  !
1026  ! Local variables.
1027  !
1028  integer :: ierr
1029 
1030  integer(kind=intType) :: ii, kk, ll, mm, nn, activeLeaf
1031  integer(kind=intType) :: nBB, nFrontLeaves, nFrontLeavesNew
1032  integer(kind=intType) :: nAllocBB, nAllocFront, nNodeElement
1033  integer(kind=intType) :: i, kkk
1034 
1035  integer(kind=intType), dimension(8) :: n, m
1036 
1037  real(kind=realtype) :: dx, dy, dz, d1, d2, invlen, val
1038  real(kind=realtype) :: u, v, w, uv, uold, vold, vn, du, dv
1039  real(kind=realtype) :: uu, vv, ww
1040 
1041  real(kind=realtype), dimension(2) :: dd
1042  real(kind=realtype), dimension(3) :: x1, x21, x41, x3142, xf
1043  real(kind=realtype), dimension(3) :: vf, vt, a, b, norm, an, bn
1044  real(kind=realtype), dimension(3) :: chi
1045  real(kind=realtype), dimension(8) :: weight
1046 
1047  real(kind=realtype), dimension(:, :), pointer :: xbbox
1048 
1049  logical :: elementFound
1050  type(adtleaftype), dimension(:), pointer :: ADTree
1051 
1052  ! Set some pointers to make the code more readable.
1053 
1054  xbbox => adt%xBBox
1055  adtree => adt%ADTree
1056 
1057  ! Initial allocation of the arrays for the tree traversal as well
1058  ! as the stack array used in the qsort routine. The latter is
1059  ! done, because the qsort routine is called for every coordinate
1060  ! and therefore it is more efficient to allocate the stack once
1061  ! rather than over and over again. The disadvantage of course is
1062  ! that an essentially local variable, stack, is now stored in
1063  ! adtData.
1064 
1065  nallocbb = size(bb)
1066  nallocfront = size(frontleaves)
1067  nstack = size(stack)
1068 
1069  ! Initialize the processor ID to -1 to indicate that no
1070  ! corresponding volume element is found and the new minimum
1071  ! distance squared to the old value.
1072 
1073  intinfo(1) = -1
1074  uvw(4) = coor(4)
1075  !
1076  ! Part 1. Determine the possible minimum distance squared to
1077  ! the root leaf. If larger than the current distance
1078  ! there is no need to search this tree.
1079  !
1080  if (coor(1) < adtree(1)%xMin(1)) then
1081  dx = coor(1) - adtree(1)%xMin(1)
1082  else if (coor(1) > adtree(1)%xMax(4)) then
1083  dx = coor(1) - adtree(1)%xMax(4)
1084  else
1085  dx = zero
1086  end if
1087 
1088  if (coor(2) < adtree(1)%xMin(2)) then
1089  dy = coor(2) - adtree(1)%xMin(2)
1090  else if (coor(2) > adtree(1)%xMax(5)) then
1091  dy = coor(2) - adtree(1)%xMax(5)
1092  else
1093  dy = zero
1094  end if
1095 
1096  if (coor(3) < adtree(1)%xMin(3)) then
1097  dz = coor(3) - adtree(1)%xMin(3)
1098  else if (coor(3) > adtree(1)%xMax(6)) then
1099  dz = coor(3) - adtree(1)%xMax(6)
1100  else
1101  dz = zero
1102  end if
1103 
1104  ! Continue with the next coordinate if the possible distance
1105  ! squared to the root leaf is larger than the currently stored
1106  ! value.
1107 
1108  if ((dx * dx + dy * dy + dz * dz) >= uvw(4)) return
1109  !
1110  ! Part 2. Find a likely bounding box, which minimizes the
1111  ! guaranteed distance.
1112  !
1113  activeleaf = 1
1114 
1115  ! Traverse the tree until a terminal leaf is found.
1116 
1117  treetraversal1: do
1118 
1119  ! Exit the loop when a terminal leaf has been found.
1120  ! This is indicated by a negative value of activeLeaf.
1121 
1122  if (activeleaf < 0) exit treetraversal1
1123 
1124  ! Determine the guaranteed distance squared for both children.
1125 
1126  do mm = 1, 2
1127 
1128  ! Determine whether the child contains a bounding box or
1129  ! another leaf of the tree.
1130 
1131  ll = adtree(activeleaf)%children(mm)
1132  if (ll > 0) then
1133 
1134  ! Child contains a leaf. Determine the guaranteed distance
1135  ! vector to the leaf.
1136 
1137  d1 = abs(coor(1) - adtree(ll)%xMin(1))
1138  d2 = abs(coor(1) - adtree(ll)%xMax(4))
1139  dx = max(d1, d2)
1140 
1141  d1 = abs(coor(2) - adtree(ll)%xMin(2))
1142  d2 = abs(coor(2) - adtree(ll)%xMax(5))
1143  dy = max(d1, d2)
1144 
1145  d1 = abs(coor(3) - adtree(ll)%xMin(3))
1146  d2 = abs(coor(3) - adtree(ll)%xMax(6))
1147  dz = max(d1, d2)
1148 
1149  else
1150 
1151  ! Child contains a bounding box. Determine the guaranteed
1152  ! distance vector to it.
1153 
1154  ll = -ll
1155 
1156  d1 = abs(coor(1) - xbbox(1, ll))
1157  d2 = abs(coor(1) - xbbox(4, ll))
1158  dx = max(d1, d2)
1159 
1160  d1 = abs(coor(2) - xbbox(2, ll))
1161  d2 = abs(coor(2) - xbbox(5, ll))
1162  dy = max(d1, d2)
1163 
1164  d1 = abs(coor(3) - xbbox(3, ll))
1165  d2 = abs(coor(3) - xbbox(6, ll))
1166  dz = max(d1, d2)
1167 
1168  end if
1169 
1170  ! Compute the guaranteed distance squared for this child.
1171 
1172  dd(mm) = dx * dx + dy * dy + dz * dz
1173 
1174  end do
1175 
1176  ! Determine which will be the next leaf in the tree traversal.
1177  ! This will be the leaf which has the minimum guaranteed
1178  ! distance. In case of ties take the left leaf, because this
1179  ! leaf may have more children.
1180 
1181  if (dd(1) <= dd(2)) then
1182  activeleaf = adtree(activeleaf)%children(1)
1183  else
1184  activeleaf = adtree(activeleaf)%children(2)
1185  end if
1186 
1187  end do treetraversal1
1188 
1189  ! Store the minimum of the just computed guaranteed distance
1190  ! squared and the currently stored value in uvw.
1191 
1192  uvw(4) = min(uvw(4), dd(1), dd(2))
1193  !
1194  ! Part 3. Find the bounding boxes whose possible minimum
1195  ! distance is less than the currently stored value.
1196  !
1197  ! In part 1 it was already tested that the possible distance
1198  ! squared of the root leaf was less than the current value.
1199  ! Therefore initialize the current front to the root leaf and
1200  ! set the number of bounding boxes to 0.
1201 
1202  nbb = 0
1203 
1204  nfrontleaves = 1
1205  frontleaves(1) = 1
1206 
1207  ! Second tree traversal. Now to find all possible bounding
1208  ! box candidates.
1209 
1210  treetraversal2: do
1211 
1212  ! Initialize the number of leaves for the new front, i.e.
1213  ! the front of the next round, to 0.
1214 
1215  nfrontleavesnew = 0
1216 
1217  ! Loop over the leaves of the current front.
1218 
1219  currentfrontloop: do ii = 1, nfrontleaves
1220 
1221  ! Store the ID of the leaf a bit easier and loop over
1222  ! its two children.
1223 
1224  ll = frontleaves(ii)
1225 
1226  childrenloop: do mm = 1, 2
1227 
1228  ! Determine whether this child contains a bounding box
1229  ! or a leaf of the next level.
1230 
1231  kk = adtree(ll)%children(mm)
1232  terminaltest: if (kk < 0) then
1233 
1234  ! Child contains a bounding box. Determine the possible
1235  ! minimum distance squared to this bounding box.
1236 
1237  kk = -kk
1238 
1239  if (coor(1) < xbbox(1, kk)) then
1240  dx = coor(1) - xbbox(1, kk)
1241  else if (coor(1) > xbbox(4, kk)) then
1242  dx = coor(1) - xbbox(4, kk)
1243  else
1244  dx = zero
1245  end if
1246 
1247  if (coor(2) < xbbox(2, kk)) then
1248  dy = coor(2) - xbbox(2, kk)
1249  else if (coor(2) > xbbox(5, kk)) then
1250  dy = coor(2) - xbbox(5, kk)
1251  else
1252  dy = zero
1253  end if
1254 
1255  if (coor(3) < xbbox(3, kk)) then
1256  dz = coor(3) - xbbox(3, kk)
1257  else if (coor(3) > xbbox(6, kk)) then
1258  dz = coor(3) - xbbox(6, kk)
1259  else
1260  dz = zero
1261  end if
1262 
1263  d2 = dx * dx + dy * dy + dz * dz
1264 
1265  ! If this distance squared is less than the current
1266  ! value, store this bounding box as a target.
1267 
1268  teststorebbox: if (d2 < uvw(4)) then
1269 
1270  ! Check if the memory must be reallocated.
1271 
1272  if (nbb == nallocbb) &
1273  call reallocbboxtargettypeplus(bb, nallocbb, &
1274  100, adt)
1275 
1276  ! Update the counter and store the data.
1277 
1278  nbb = nbb + 1
1279  bb(nbb)%ID = kk
1280  bb(nbb)%posDist2 = d2
1281 
1282  ! Although in step 2, i.e. the first tree traversal,
1283  ! the guaranteed distance squared to a bounding box
1284  ! has already been computed, this has been done only
1285  ! for a likely candidate and not for all the possible
1286  ! candidates. As this test is relatively cheap, do it
1287  ! now for this bounding box.
1288 
1289  d1 = abs(coor(1) - xbbox(1, kk))
1290  d2 = abs(coor(1) - xbbox(4, kk))
1291  dx = max(d1, d2)
1292 
1293  d1 = abs(coor(2) - xbbox(2, kk))
1294  d2 = abs(coor(2) - xbbox(5, kk))
1295  dy = max(d1, d2)
1296 
1297  d1 = abs(coor(3) - xbbox(3, kk))
1298  d2 = abs(coor(3) - xbbox(6, kk))
1299  dz = max(d1, d2)
1300 
1301  d2 = dx * dx + dy * dy + dz * dz
1302  uvw(4) = min(uvw(4), d2)
1303 
1304  end if teststorebbox
1305 
1306  else terminaltest
1307 
1308  ! Child contains a leaf. Compute the possible minimum
1309  ! distance squared to the current coordinate.
1310 
1311  if (coor(1) < adtree(kk)%xMin(1)) then
1312  dx = coor(1) - adtree(kk)%xMin(1)
1313  else if (coor(1) > adtree(kk)%xMax(4)) then
1314  dx = coor(1) - adtree(kk)%xMax(4)
1315  else
1316  dx = zero
1317  end if
1318 
1319  if (coor(2) < adtree(kk)%xMin(2)) then
1320  dy = coor(2) - adtree(kk)%xMin(2)
1321  else if (coor(2) > adtree(kk)%xMax(5)) then
1322  dy = coor(2) - adtree(kk)%xMax(5)
1323  else
1324  dy = zero
1325  end if
1326 
1327  if (coor(3) < adtree(kk)%xMin(3)) then
1328  dz = coor(3) - adtree(kk)%xMin(3)
1329  else if (coor(3) > adtree(kk)%xMax(6)) then
1330  dz = coor(3) - adtree(kk)%xMax(6)
1331  else
1332  dz = zero
1333  end if
1334 
1335  d2 = dx * dx + dy * dy + dz * dz
1336 
1337  ! If this distance squared is less than the current
1338  ! value, store this leaf in the new front.
1339 
1340  teststoreleave: if (d2 < uvw(4)) then
1341 
1342  ! Check if enough memory has been allocated and
1343  ! store the leaf.
1344 
1345  if (nfrontleavesnew == nallocfront) then
1346  i = nallocfront
1347  call reallocplus(frontleavesnew, i, 250, adt)
1348  call reallocplus(frontleaves, nallocfront, 250, adt)
1349  end if
1350 
1351  nfrontleavesnew = nfrontleavesnew + 1
1352  frontleavesnew(nfrontleavesnew) = kk
1353 
1354  ! Compute the guaranteed distance squared to this leaf.
1355  ! It may be less than the currently stored value.
1356 
1357  d1 = abs(coor(1) - adtree(kk)%xMin(1))
1358  d2 = abs(coor(1) - adtree(kk)%xMax(4))
1359  dx = max(d1, d2)
1360 
1361  d1 = abs(coor(2) - adtree(kk)%xMin(2))
1362  d2 = abs(coor(2) - adtree(kk)%xMax(5))
1363  dy = max(d1, d2)
1364 
1365  d1 = abs(coor(3) - adtree(kk)%xMin(3))
1366  d2 = abs(coor(3) - adtree(kk)%xMax(6))
1367  dz = max(d1, d2)
1368 
1369  d2 = dx * dx + dy * dy + dz * dz
1370  uvw(4) = min(uvw(4), d2)
1371 
1372  end if teststoreleave
1373 
1374  end if terminaltest
1375  end do childrenloop
1376  end do currentfrontloop
1377 
1378  ! End of the loop over the current front. If the new front
1379  ! is empty the entire tree has been traversed and an exit is
1380  ! made from the corresponding loop.
1381 
1382  if (nfrontleavesnew == 0) exit treetraversal2
1383 
1384  ! Copy the data of the new front leaves into the current
1385  ! front for the next round.
1386 
1387  nfrontleaves = nfrontleavesnew
1388  do ll = 1, nfrontleaves
1389  frontleaves(ll) = frontleavesnew(ll)
1390  end do
1391 
1392  end do treetraversal2
1393 
1394  ! Sort the target bounding boxes in increasing order such that
1395  ! the one with the smallest possible distance is first.
1396 
1397  call qsortbboxtargets(bb, nbb, adt)
1398  !
1399  ! Part 4: Loop over the selected bounding boxes and check if
1400  ! the corresponding element minimizes the distance.
1401  !
1402  elementfound = .false.
1403 
1404  bboxloop: do mm = 1, nbb
1405 
1406  ! Exit the loop if the possible minimum distance of this
1407  ! bounding box is not smaller than the current value.
1408  ! Remember that BB has been sorted in increasing order.
1409 
1410  if (uvw(4) <= bb(mm)%posDist2) exit bboxloop
1411 
1412  ! Determine the element type stored in this bounding box.
1413 
1414  kk = bb(mm)%ID
1415  select case (adt%elementType(kk))
1416 
1417  case (adttriangle)
1418  call adtterminate(adt, "minDistanceTreeSearch", &
1419  "Minimum distance search for &
1420  &triangles not implemented yet")
1421 
1422  !=========================================================
1423 
1424  case (adtquadrilateral)
1425 
1426  ! Temporary implementation. I'm waiting for Juan to come
1427  ! up with his more sophisticated algorithm.
1428 
1429  ! This is a surface element, so set the parametric weight
1430  ! w to zero.
1431 
1432  w = zero
1433 
1434  ! Determine the 4 vectors which completely describe
1435  ! the quadrilateral face
1436 
1437  ll = adt%elementID(kk)
1438  n(1) = adt%quadsConn(1, ll)
1439  n(2) = adt%quadsConn(2, ll)
1440  n(3) = adt%quadsConn(3, ll)
1441  n(4) = adt%quadsConn(4, ll)
1442 
1443  x1(1) = adt%coor(1, n(1))
1444  x1(2) = adt%coor(2, n(1))
1445  x1(3) = adt%coor(3, n(1))
1446 
1447  x21(1) = adt%coor(1, n(2)) - x1(1)
1448  x21(2) = adt%coor(2, n(2)) - x1(2)
1449  x21(3) = adt%coor(3, n(2)) - x1(3)
1450 
1451  x41(1) = adt%coor(1, n(4)) - x1(1)
1452  x41(2) = adt%coor(2, n(4)) - x1(2)
1453  x41(3) = adt%coor(3, n(4)) - x1(3)
1454 
1455  x3142(1) = adt%coor(1, n(3)) - x1(1) - x21(1) - x41(1)
1456  x3142(2) = adt%coor(2, n(3)) - x1(2) - x21(2) - x41(2)
1457  x3142(3) = adt%coor(3, n(3)) - x1(3) - x21(3) - x41(3)
1458 
1459  ! Initialize u and v to 0.5 and determine the
1460  ! corresponding coordinates on the face, which is the
1461  ! centroid.
1462 
1463  u = half
1464  v = half
1465  uv = u * v
1466 
1467  xf(1) = x1(1) + u * x21(1) + v * x41(1) + uv * x3142(1)
1468  xf(2) = x1(2) + u * x21(2) + v * x41(2) + uv * x3142(2)
1469  xf(3) = x1(3) + u * x21(3) + v * x41(3) + uv * x3142(3)
1470 
1471  ! Newton loop to determine the point on the surface,
1472  ! which minimizes the distance to the given coordinate.
1473 
1474  newtonquads: do ll = 1, itermax
1475 
1476  ! Store the current values of u and v for a stop
1477  ! criterion later on.
1478 
1479  uold = u
1480  vold = v
1481 
1482  ! Determine the vector vf from xf to given coordinate.
1483 
1484  vf(1) = coor(1) - xf(1)
1485  vf(2) = coor(2) - xf(2)
1486  vf(3) = coor(3) - xf(3)
1487 
1488  ! Determine the tangent vectors in u- and v-direction.
1489  ! Store these in a and b respectively.
1490 
1491  a(1) = x21(1) + v * x3142(1)
1492  a(2) = x21(2) + v * x3142(2)
1493  a(3) = x21(3) + v * x3142(3)
1494 
1495  b(1) = x41(1) + u * x3142(1)
1496  b(2) = x41(2) + u * x3142(2)
1497  b(3) = x41(3) + u * x3142(3)
1498 
1499  ! Determine the normal vector of the face by taking the
1500  ! cross product of a and b. Afterwards this vector will
1501  ! be scaled to a unit vector.
1502 
1503  norm(1) = a(2) * b(3) - a(3) * b(2)
1504  norm(2) = a(3) * b(1) - a(1) * b(3)
1505  norm(3) = a(1) * b(2) - a(2) * b(1)
1506 
1507  invlen = one / max(adteps, sqrt(norm(1) * norm(1) &
1508  + norm(2) * norm(2) &
1509  + norm(3) * norm(3)))
1510 
1511  norm(1) = norm(1) * invlen
1512  norm(2) = norm(2) * invlen
1513  norm(3) = norm(3) * invlen
1514 
1515  ! Determine the projection of the vector vf onto
1516  ! the face.
1517 
1518  vn = vf(1) * norm(1) + vf(2) * norm(2) + vf(3) * norm(3)
1519  vt(1) = vf(1) - vn * norm(1)
1520  vt(2) = vf(2) - vn * norm(2)
1521  vt(3) = vf(3) - vn * norm(3)
1522 
1523  ! The vector vt points from the current point on the
1524  ! face to the new point. However this new point lies on
1525  ! the plane determined by the vectors a and b, but not
1526  ! necessarily on the face itself. The new point on the
1527  ! face is obtained by projecting the point in the a-b
1528  ! plane onto the face. this can be done by determining
1529  ! the coefficients du and dv, such that vt = du*a + dv*b.
1530  ! To solve du and dv the vectors normal to a and b
1531  ! inside the plane ab are needed.
1532 
1533  an(1) = a(2) * norm(3) - a(3) * norm(2)
1534  an(2) = a(3) * norm(1) - a(1) * norm(3)
1535  an(3) = a(1) * norm(2) - a(2) * norm(1)
1536 
1537  bn(1) = b(2) * norm(3) - b(3) * norm(2)
1538  bn(2) = b(3) * norm(1) - b(1) * norm(3)
1539  bn(3) = b(1) * norm(2) - b(2) * norm(1)
1540 
1541  ! Solve du and dv. the clipping of vn should not be
1542  ! active, as this would mean that the vectors a and b
1543  ! are parallel. This corresponds to a quad degenerated
1544  ! to a line, which should not occur in the surface mesh.
1545 
1546  vn = a(1) * bn(1) + a(2) * bn(2) + a(3) * bn(3)
1547  vn = sign(max(adteps, abs(vn)), vn)
1548  du = (vt(1) * bn(1) + vt(2) * bn(2) + vt(3) * bn(3)) / vn
1549 
1550  vn = b(1) * an(1) + b(2) * an(2) + b(3) * an(3)
1551  vn = sign(max(adteps, abs(vn)), vn)
1552  dv = (vt(1) * an(1) + vt(2) * an(2) + vt(3) * an(3)) / vn
1553 
1554  ! Determine the new parameter values uu and vv. These
1555  ! are limited to 0 <= (uu,vv) <= 1.
1556 
1557  u = u + du; u = min(one, max(zero, u))
1558  v = v + dv; v = min(one, max(zero, v))
1559 
1560  ! Determine the final values of the corrections.
1561 
1562  du = abs(u - uold)
1563  dv = abs(v - vold)
1564 
1565  ! Determine the new coordinates of the point xf.
1566 
1567  uv = u * v
1568  xf(1) = x1(1) + u * x21(1) + v * x41(1) + uv * x3142(1)
1569  xf(2) = x1(2) + u * x21(2) + v * x41(2) + uv * x3142(2)
1570  xf(3) = x1(3) + u * x21(3) + v * x41(3) + uv * x3142(3)
1571 
1572  ! Exit the loop if the update of the parametric
1573  ! weights is below the threshold
1574 
1575  val = sqrt(du * du + dv * dv)
1576  if (val <= thresconv) exit newtonquads
1577 
1578  end do newtonquads
1579 
1580  ! Compute the distance squared between the given
1581  ! coordinate and the point xf.
1582 
1583  dx = coor(1) - xf(1)
1584  dy = coor(2) - xf(2)
1585  dz = coor(3) - xf(3)
1586 
1587  val = dx * dx + dy * dy + dz * dz
1588 
1589  ! If the distance squared is less than the current value
1590  ! store the wall distance and interpolation info and
1591  ! indicate that an element was found.
1592 
1593  if (val < uvw(4)) then
1594  uvw(4) = val
1595  nnodeelement = 4
1596  elementfound = .true.
1597 
1598  kkk = kk; uu = u; vv = v; ww = w
1599  m(1) = n(1); m(2) = n(2); m(3) = n(3); m(4) = n(4)
1600 
1601  weight(1) = (one - u) * (one - v)
1602  weight(2) = u * (one - v)
1603  weight(3) = u * v
1604  weight(4) = (one - u) * v
1605  end if
1606 
1607  !=========================================================
1608 
1609  case (adttetrahedron)
1610  call adtterminate(adt, "minDistanceTreeSearch", &
1611  "Minimum distance search for &
1612  &tetrahedra not implemented yet")
1613 
1614  !===========================================================
1615 
1616  case (adtpyramid)
1617  call adtterminate(adt, "minDistanceTreeSearch", &
1618  "Minimum distance search for &
1619  &pyramids not implemented yet")
1620 
1621  !===========================================================
1622 
1623  case (adtprism)
1624  call adtterminate(adt, "minDistanceTreeSearch", &
1625  "Minimum distance search for &
1626  &prisms not implemented yet")
1627 
1628  !===========================================================
1629 
1630  case (adthexahedron)
1631 
1632  ! Determine the element ID and the corresponding
1633  ! 8 node ID's.
1634 
1635  ll = adt%elementID(kk)
1636  n(1) = adt%hexaConn(1, ll)
1637  n(2) = adt%hexaConn(2, ll)
1638  n(3) = adt%hexaConn(3, ll)
1639  n(4) = adt%hexaConn(4, ll)
1640  n(5) = adt%hexaConn(5, ll)
1641  n(6) = adt%hexaConn(6, ll)
1642  n(7) = adt%hexaConn(7, ll)
1643  n(8) = adt%hexaConn(8, ll)
1644 
1645  ! Call the subroutine minD2Hexa to do the work.
1646 
1647  call mind2hexa(coor(1:3), &
1648  adt%coor(1:3, n(1)), &
1649  adt%coor(1:3, n(2)), &
1650  adt%coor(1:3, n(3)), &
1651  adt%coor(1:3, n(4)), &
1652  adt%coor(1:3, n(5)), &
1653  adt%coor(1:3, n(6)), &
1654  adt%coor(1:3, n(7)), &
1655  adt%coor(1:3, n(8)), &
1656  val, chi, ierr)
1657 
1658  ! If the distance squared is less than the current value
1659  ! store the wall distance and interpolation info and
1660  ! indicate that an element was found.
1661 
1662  if (val < uvw(4)) then
1663  uvw(4) = val
1664  nnodeelement = 8
1665  elementfound = .true.
1666 
1667  kkk = kk;
1668  uu = chi(1); vv = chi(2); ww = chi(3)
1669 
1670  m(1) = n(1); m(2) = n(2); m(3) = n(3); m(4) = n(4)
1671  m(5) = n(5); m(6) = n(6); m(7) = n(7); m(8) = n(8)
1672 
1673  weight(1) = (one - uu) * (one - vv) * (one - ww)
1674  weight(2) = uu * (one - vv) * (one - ww)
1675  weight(3) = uu * vv * (one - ww)
1676  weight(4) = (one - uu) * vv * (one - ww)
1677  weight(5) = (one - uu) * (one - vv) * ww
1678  weight(6) = uu * (one - vv) * ww
1679  weight(7) = uu * vv * ww
1680  weight(8) = (one - uu) * vv * ww
1681  end if
1682 
1683  end select
1684 
1685  end do bboxloop
1686 
1687  ! Check if an element was found. As all the minimum distance
1688  ! searches are initialized by the calling routine (to support
1689  ! periodic searches) this is not always the case.
1690 
1691  if (elementfound) then
1692 
1693  ! Store the interpolation information for this point.
1694  ! First the integer info, i.e. the processor ID, element type
1695  ! and local element ID.
1696 
1697  intinfo(1) = adt%myID
1698  intinfo(2) = adt%elementType(kkk)
1699  intinfo(3) = adt%elementID(kkk)
1700 
1701  ! The parametric weights. Note that the wall distance
1702  ! squared, stored in the 4th position of uvw, already
1703  ! contains the correct value.
1704 
1705  uvw(1) = uu
1706  uvw(2) = vv
1707  uvw(3) = ww
1708 
1709  ! The interpolated solution, if needed.
1710 
1711  do ll = 1, ninterpol
1712  ii = 4 + ll
1713  uvw(ii) = weight(1) * arrdonor(ll, m(1))
1714  do i = 2, nnodeelement
1715  uvw(ii) = uvw(ii) + weight(i) * arrdonor(ll, m(i))
1716  end do
1717  end do
1718 
1719  end if
1720 
1721  end subroutine mindistancetreesearchsinglepoint
1722 
1723  subroutine intersectiontreesearchsinglepoint(ADT, coor, &
1724  intInfo, BB, frontLeaves, frontLeavesNew)
1725  !
1726  ! This routine is used in the ray casting approach to determine
1727  ! if a given ray intersects any of the surface elements. The
1728  ! purpose is to determine if a point of interest is inside
1729  ! or outside the (closed) surface defined by the ADT.
1730  ! Subroutine intent(in) arguments.
1731  ! --------------------------------
1732  ! ADT: ADT type whose ADT must be searched
1733  ! coor(3): The coordinate of the point to be searched.
1734  ! Subroutine intent(out) arguments.
1735  ! ---------------------------------
1736  ! intInfo: Intersection info. The number of intersections we
1737  ! * found.
1738  !
1739  implicit none
1740  !
1741  ! Subroutine arguments.
1742  !
1743  type(adttype), intent(inout) :: ADT
1744 
1745  real(kind=realtype), dimension(3), intent(in) :: coor
1746  integer(kind=intType), intent(out) :: intInfo
1747 
1748  integer(kind=intType), dimension(:), pointer :: BB
1749  integer(kind=intType), dimension(:), pointer :: frontLeaves
1750  integer(kind=intType), dimension(:), pointer :: frontLeavesNew
1751  !
1752  ! Local variables.
1753  !
1754  integer :: ierr
1755 
1756  integer(kind=intType) :: ii, kk, ll, mm, nn
1757  integer(kind=intType) :: nBB, nFrontLeaves, nFrontLeavesNew
1758  integer(kind=intType) :: nAllocBB, nAllocFront
1759  integer(kind=intType) :: i, nNodeElement
1760  real(kind=realtype), dimension(:, :), pointer :: xbbox
1761  logical :: elementFound
1762  type(adtleaftype), dimension(:), pointer :: ADTree
1763 
1764  ! Set some pointers to make the code more readable.
1765 
1766  xbbox => adt%xBBox
1767  adtree => adt%ADTree
1768 
1769  ! Determine the sizes from the arrays we have been passed
1770 
1771  nallocbb = size(bb)
1772  nallocfront = size(frontleaves)
1773 
1774  ! Initialize the number of possible intersections to 0
1775 
1776  intinfo = 0
1777  !
1778  ! Part 1. Traverse the tree and determine the target
1779  ! bounding boxes, which may contain the intersection
1780  !
1781  ! Start at the root, i.e. set the front leaf to the root leaf.
1782  ! Also initialize the number of possible bounding boxes to 0.
1783 
1784  nbb = 0
1785 
1786  nfrontleaves = 1
1787  frontleaves(1) = 1
1788 
1789  treetraversalloop: do
1790 
1791  ! Initialize the number of leaves for the new front, i.e.
1792  ! the front of the next round, to 0.
1793 
1794  nfrontleavesnew = 0
1795 
1796  ! Loop over the leaves of the current front.
1797 
1798  currentfrontloop: do ii = 1, nfrontleaves
1799 
1800  ! Store the ID of the leaf a bit easier and loop over
1801  ! its two children.
1802 
1803  ll = frontleaves(ii)
1804 
1805  childrenloop: do mm = 1, 2
1806 
1807  ! Determine whether this child contains a bounding box
1808  ! or a leaf of the next level.
1809 
1810  kk = adtree(ll)%children(mm)
1811  terminaltest: if (kk < 0) then
1812 
1813  ! Child contains a bounding box. Check if the
1814  ! ray is inside the bounding box.
1815 
1816  kk = -kk
1817  if (coor(1) <= xbbox(4, kk) .and. &
1818  coor(2) >= xbbox(2, kk) .and. &
1819  coor(2) <= xbbox(5, kk) .and. &
1820  coor(3) >= xbbox(3, kk) .and. &
1821  coor(3) <= xbbox(6, kk)) then
1822 
1823  ! Ray intersectst he bounding box. That's all we
1824  ! wanted to know:
1825  intinfo = 1
1826  exit treetraversalloop
1827  end if
1828  else terminaltest
1829 
1830  ! Child contains a leaf. Check if the coordinate is
1831  ! inside the bounding box of the leaf.
1832 
1833  if (coor(1) <= adtree(kk)%xMax(4) .and. &
1834  coor(2) >= adtree(kk)%xMin(2) .and. &
1835  coor(2) <= adtree(kk)%xMax(5) .and. &
1836  coor(3) >= adtree(kk)%xMin(3) .and. &
1837  coor(3) <= adtree(kk)%xMax(6)) then
1838 
1839  ! Coordinate is inside the leaf. Store the leaf in
1840  ! the list for the new front.
1841 
1842  if (nfrontleavesnew == nallocfront) then
1843  i = nallocfront
1844  call reallocplus(frontleavesnew, i, 250, adt)
1845  call reallocplus(frontleaves, nallocfront, 250, adt)
1846  end if
1847 
1848  nfrontleavesnew = nfrontleavesnew + 1
1849  frontleavesnew(nfrontleavesnew) = kk
1850  end if
1851 
1852  end if terminaltest
1853 
1854  end do childrenloop
1855 
1856  end do currentfrontloop
1857 
1858  ! End of the loop over the current front. If the new front
1859  ! is empty the entire tree has been traversed and an exit is
1860  ! made from the corresponding loop.
1861 
1862  if (nfrontleavesnew == 0) then
1863  exit treetraversalloop
1864  end if
1865  ! Copy the data of the new front leaves into the current
1866  ! front for the next round.
1867 
1868  nfrontleaves = nfrontleavesnew
1869  do ll = 1, nfrontleaves
1870  frontleaves(ll) = frontleavesnew(ll)
1871  end do
1872 
1873  end do treetraversalloop
1874  !
1875  ! Part 2: Loop over the selected bounding boxes and check if
1876  ! the corresponding elements contain the point.
1877  !
1878  !intInfo = nBB
1879  ! elementFound = .false.
1880 
1881  ! BBoxLoop: do mm=1,nBB
1882 
1883  ! ! Determine the element type stored in this bounding box.
1884 
1885  ! kk = BB(mm)
1886  ! select case (ADT%elementType(kk))
1887 
1888  ! case (adtQuadrilateral)
1889 
1890  ! !=========================================================
1891 
1892  ! case (adtTriangle)
1893 
1894  ! end select
1895 
1896  ! enddo BBoxLoop
1897 
1898  end subroutine intersectiontreesearchsinglepoint
1899 
1900  subroutine mind2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, d2, chi, iErr)
1901  !
1902  ! Subroutine to compute a fail-safe minimum distance computation
1903  ! between a point, P, and a hexahedral element. This subroutine
1904  ! can provide the minimum distance whether P is inside or
1905  ! outside of the hexahedral element. If P is inside the element
1906  ! the distance returned is zero, while if P is outside of the
1907  ! element, the distance returned is the minimum distance between
1908  ! P and any of the bounding faces of the hexahedral element.
1909  ! The basic idea of this subroutine is to perform a
1910  ! bound-constrained minimization of the distance between the
1911  ! point P and a point inside the element given by parametric
1912  ! coordinates chi(3) using a modified Newton step bound
1913  ! constrained optimizer.
1914  ! For more details of the optimization algorithms, a good
1915  ! reference is:
1916  ! Nocedal, J. and Wright, S. J., Numerical Optimization,
1917  ! Springer, 1999.
1918  ! Subroutine arguments:
1919  ! xp(3) - The Cartesian coordinates of point P.
1920  ! x1(3)-x8(3) - The Cartesian coordinates of the 8 nodes that
1921  ! make up the hexahedron, in the order specified
1922  ! in the CHIMPS standard.
1923  ! d2 - The squared of the minimum distance between the
1924  ! point P and the hexahedral element.
1925  ! chi(3) - Parametric coordinates of the point that
1926  ! belongs to the hexahedron where the minimum
1927  ! distance has been found. If P is inside the
1928  ! element, 0<(ksi,eta,zeta)<1, while if P is
1929  ! strictly outside of the element, then one or
1930  ! more of the values of (ksi,eta,zeta) will be
1931  ! exactly zero or one.
1932  ! iErr - Output status of this subroutine:
1933  ! iErr = 0, Proper minimum found.
1934  ! iErr = -1, Distance minimization failed.
1935  !
1936  use precision
1937 
1938  implicit none
1939  !
1940  ! Subroutine arguments.
1941  !
1942  real(kind=realtype), dimension(3), intent(in) :: xp
1943  real(kind=realtype), dimension(3), intent(in) :: x1, x2, x3, x4
1944  real(kind=realtype), dimension(3), intent(in) :: x5, x6, x7, x8
1945 
1946  integer(kind=intType), intent(out) :: iErr
1947 
1948  real(kind=realtype), intent(out) :: d2
1949  real(kind=realtype), dimension(3), intent(out) :: chi
1950  !
1951  ! Local variables.
1952  !
1953  integer(kind=intType), parameter :: maxIt = 30
1954  integer(kind=intType) :: i, itCount = 0
1955  integer(kind=intType), dimension(3) :: actSet, chiGradConv
1956 
1957  real(kind=realtype) :: inactgradnorm, normdeltachi, x0, y0, z0
1958  real(kind=realtype), parameter :: gradtol = 1.0e-14_realtype
1959  real(kind=realtype), parameter :: deltachitol = 1.0e-14_realtype
1960  real(kind=realtype), dimension(3) :: deltachi, actualdeltachi
1961  real(kind=realtype), dimension(3) :: lwrbnd, uppbnd
1962  real(kind=realtype), dimension(3) :: grad, oldchi
1963  real(kind=realtype), dimension(3, 3) :: hess
1964 
1965  logical :: convDeltaChi, convGradD2
1966  !
1967  ! Initialization section.
1968  !
1969 
1970  ! Setup initial values for the parametric coordinates of the
1971  ! minimum distance point. One may be able to do better than this
1972  ! in the future, but this is probably a good guess.
1973 
1974  chi(1) = 0.5_realtype
1975  chi(2) = 0.5_realtype
1976  chi(3) = 0.5_realtype
1977 
1978  ! Initialize the active set array, actSet(i), i=1,2,3
1979  ! actSet(i) = 0 if bound constraint i is inactive
1980  ! actSet(i) = +1 if bound constraint i is active at upper bound
1981  ! actSet(i) = -1 if bound constraint i is active at lower bound
1982 
1983  actset(:) = 0
1984 
1985  ! Initialize the upper and lower bounds for all chi variables
1986 
1987  lwrbnd(:) = 0.0_realtype
1988  uppbnd(:) = 1.0_realtype
1989 
1990  ! Initialize actualDeltaChi (step size) to a large value so that the
1991  ! convergence criteria for the step size is guaranteed not to pass
1992  ! on the first iteration.
1993 
1994  actualdeltachi(:) = 1.0_realtype
1995 
1996  !
1997  ! Main iteration loop.
1998  !
1999  itcount = 0
2000 
2001  iterloop: do while (itcount <= maxit)
2002 
2003  ! Increment iteration counter
2004 
2005  itcount = itcount + 1
2006 
2007  !
2008  ! Compute the gradient of d2
2009  !
2010 
2011  call gradd2hexa(xp, x1, x2, x3, x4, x5, x6, x7, x8, chi, x0, y0, z0, grad, ierr)
2012 
2013  !
2014  ! Convergence test
2015  !
2016 
2017  convdeltachi = .false.
2018  convgradd2 = .false.
2019 
2020  ! Convergence test for step size
2021 
2022  normdeltachi = sqrt(actualdeltachi(1)**2 + actualdeltachi(2)**2 + actualdeltachi(3)**2)
2023 
2024  if (normdeltachi < deltachitol) convdeltachi = .true.
2025 
2026  ! Convergence test for the gradient. Note that this gradient
2027  ! test is such that, in order to pass it, the following must
2028  ! be true:
2029  !
2030  ! 1) \frac{\partial d2}{\partial \chi_i} > 0 for i in the active
2031  ! set at the lower bound, actSet(i) = -1
2032  ! 2) \frac{\partial d2}{\partial \chi_i} < 0 for i in the active
2033  ! set at the upper bound, actSet(i) = +1
2034  ! 3) norm of components of the gradient in the inactive set must
2035  ! be smaller than gradTol
2036  !
2037 
2038  inactgradnorm = 0.0_realtype
2039  chigradconv(:) = 0
2040 
2041  do i = 1, 3
2042 
2043  ! If this component of chi is in the inactive set, accumulate
2044  ! the total value of the gradient and deal with it later.
2045 
2046  if (actset(i) == 0) inactgradnorm = inactgradnorm + grad(i)**2
2047 
2048  ! If this component of chi is active at a lower bound, check
2049  ! for the convergence criterion. Also deactivate the constraint
2050  ! if the gradient is pointing in the proper direction.
2051 
2052  if (actset(i) == -1) then
2053  if (grad(i) > 0.0_realtype) then
2054  chigradconv(i) = 1
2055  else
2056  actset(i) = 0
2057  end if
2058  end if
2059 
2060  ! If this component of chi is active at an upper bound, check
2061  ! for the convergence criterion.
2062 
2063  if (actset(i) == 1) then
2064  if (grad(i) < 0.0_realtype) then
2065  chigradconv(i) = 1
2066  else
2067  actset(i) = 0
2068  end if
2069  end if
2070 
2071  end do
2072 
2073  ! Check for convergence on the accumulated values of the inactive
2074  ! components of the gradient.
2075 
2076  if (inactgradnorm < gradtol) then
2077  do i = 1, 3
2078  if (actset(i) == 0) chigradconv(i) = 1
2079  end do
2080  end if
2081 
2082  if (sum(chigradconv(1:3)) == 3) convgradd2 = .true.
2083 
2084  ! Test for convergence using both criteria
2085 
2086  if (convdeltachi .and. convgradd2) exit iterloop
2087 
2088  !
2089  ! Compute the Hessian of d2
2090  !
2091 
2092  call hessd2hexa(xp, x1, x2, x3, x4, x5, x6, x7, x8, chi, hess, ierr)
2093 
2094  !
2095  ! Compute the Newton step
2096  !
2097 
2098  call newtonstep(hess, grad, deltachi, ierr)
2099 
2100  !
2101  ! Update the current guess (appropriately clipped at
2102  ! the bounds) and update the active set array as needed
2103  !
2104 
2105  ! Loop over the components of chi
2106 
2107  oldchi(:) = chi(:)
2108 
2109  do i = 1, 3
2110 
2111  ! Update only the components of chi that were not active
2112 
2113  if (actset(i) == 0) then
2114 
2115  chi(i) = chi(i) + deltachi(i)
2116 
2117  ! Check to see this degree of freedom has become
2118  ! actively constrained at either an upper or
2119  ! lower bound and update active set.
2120 
2121  if (chi(i) > uppbnd(i)) then
2122  chi(i) = uppbnd(i)
2123  actset(i) = 1
2124  else if (chi(i) < lwrbnd(i)) then
2125  chi(i) = lwrbnd(i)
2126  actset(i) = -1
2127  end if
2128  end if
2129  end do
2130 
2131  actualdeltachi(:) = chi(:) - oldchi(:)
2132  end do iterloop
2133 
2134  !
2135  ! Return the results to the calling function and print info
2136  ! for debugging purposes.
2137  !
2138 
2139  ! Compute the minimum distance (squared) that the algorithm found
2140 
2141  d2 = (xp(1) - x0)**2 + (xp(2) - y0)**2 + (xp(3) - z0)**2
2142 
2143  ! Print some stuff out to the screen for debugging purposes
2144 
2145  ! write(*,*)
2146  ! write(*,*) 'Results of minD2Hexa'
2147  ! write(*, coordinateFormat)'Point P = (',xP(1),xP(2),xP(3),' )'
2148  ! write(*, coordinateFormat)'Found point = (',x0,y0,z0,' )'
2149  ! write(*, coordinateFormat)'Parametric coordinates = (',chi(1),chi(2),chi(3),' )'
2150  ! write(*, distanceFormat)'Minimum distance =',sqrt(d2)
2151  ! write(*, iterationFormat)'Number of iterations =',itCount
2152 
2153  ! character(len=maxStringLen) :: iterationFormat = '(A, 1x, i3)'
2154  ! character(len=maxStringLen) :: coordinateFormat = '(A, 3f10.6, A)'
2155  ! character(len=maxStringLen) :: distanceFormat = '(A, f20.17, A)'
2156 
2157  return
2158 
2159  end subroutine mind2hexa
2160 
2161  subroutine newtonstep(hess, grad, step, iErr)
2162  !
2163  ! Compute the Newton step given by the Hessian matrix and the
2164  ! gradient vector of the distance squared function.
2165  !
2166  use precision
2167 
2168  implicit none
2169  !
2170  ! Subroutine arguments.
2171  !
2172  real(kind=realtype), dimension(3), intent(in) :: grad
2173  real(kind=realtype), dimension(3), intent(out) :: step
2174  real(kind=realtype), dimension(3, 3), intent(in) :: hess
2175 
2176  integer(kind=intType), intent(out) :: iErr
2177 
2178  !
2179  ! Local variables.
2180  !
2181  real(kind=realtype) :: determinant
2182 
2183  !
2184  ! Compute the Newton step as the solution of the problem
2185  ! Hessian . step = -grad
2186  ! using simple Cramer's rule
2187  !
2188  ! Compute the determinant of the Hessian Matrix
2189 
2190  determinant = hess(1, 1) * (hess(2, 2) * hess(3, 3) - hess(2, 3) * hess(3, 2)) &
2191  + hess(1, 2) * (hess(2, 3) * hess(3, 1) - hess(2, 1) * hess(3, 3)) &
2192  + hess(1, 3) * (hess(2, 1) * hess(3, 2) - hess(2, 2) * hess(3, 1))
2193 
2194  ! First component of the step
2195 
2196  step(1) = grad(1) * (hess(2, 2) * hess(3, 3) - hess(2, 3) * hess(3, 2)) &
2197  + hess(1, 2) * (hess(2, 3) * grad(3) - grad(2) * hess(3, 3)) &
2198  + hess(1, 3) * (grad(2) * hess(3, 2) - hess(2, 2) * grad(3))
2199  step(1) = -step(1) / determinant
2200 
2201  ! Second component of the step
2202 
2203  step(2) = hess(1, 1) * (grad(2) * hess(3, 3) - hess(2, 3) * grad(3)) &
2204  + grad(1) * (hess(2, 3) * hess(3, 1) - hess(2, 1) * hess(3, 3)) &
2205  + hess(1, 3) * (hess(2, 1) * grad(3) - grad(2) * hess(3, 1))
2206  step(2) = -step(2) / determinant
2207 
2208  ! First component of the step
2209 
2210  step(3) = hess(1, 1) * (hess(2, 2) * grad(3) - grad(2) * hess(3, 2)) &
2211  + hess(1, 2) * (grad(2) * hess(3, 1) - hess(2, 1) * grad(3)) &
2212  + grad(1) * (hess(2, 1) * hess(3, 2) - hess(2, 2) * hess(3, 1))
2213  step(3) = -step(3) / determinant
2214 
2215  ! Return iErr = 0 for the time being.
2216 
2217  ierr = 0
2218 
2219  return
2220 
2221  end subroutine newtonstep
2222 
2223  subroutine hessd2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, chi, hess, iErr)
2224  !
2225  ! Compute the Hessian of the square of the distance between the
2226  ! point xP, and the actual point in the hexahedron represented
2227  ! by chi(1:3).
2228  !
2229  use precision
2230 
2231  implicit none
2232  !
2233  ! Subroutine arguments.
2234  !
2235  real(kind=realtype), dimension(3), intent(in) :: xp, chi
2236  real(kind=realtype), dimension(3), intent(in) :: x1, x2, x3, x4
2237  real(kind=realtype), dimension(3), intent(in) :: x5, x6, x7, x8
2238  real(kind=realtype), dimension(3, 3), intent(out) :: hess
2239 
2240  integer(kind=intType), intent(out) :: iErr
2241 
2242  !
2243  ! Local variables.
2244  !
2245  integer(kind=intType) :: i
2246 
2247  real(kind=realtype) :: ksi, eta, zeta, x0, y0, z0
2248  real(kind=realtype), dimension(8, 3) :: alpha
2249  real(kind=realtype) :: dxdksi, dxdeta, dxdzeta
2250  real(kind=realtype) :: dydksi, dydeta, dydzeta
2251  real(kind=realtype) :: dzdksi, dzdeta, dzdzeta
2252  real(kind=realtype) :: d2xdksideta, d2xdksidzeta, d2xdetadzeta
2253  real(kind=realtype) :: d2ydksideta, d2ydksidzeta, d2ydetadzeta
2254  real(kind=realtype) :: d2zdksideta, d2zdksidzeta, d2zdetadzeta
2255 
2256  ! Initialize the parametric coordinates with a more recognizable
2257  ! name.
2258 
2259  ksi = chi(1)
2260  eta = chi(2)
2261  zeta = chi(3)
2262 
2263  ! Initialize the alpha array (obtained by regrouping terms in the
2264  ! parametric expansion of the (ksi,eta,zeta)->(x,y,z) mapping.
2265  ! The second index of the alpha array corresponds to the x, y, or
2266  ! z coordinate mapping.
2267 
2268  do i = 1, 3
2269  alpha(1, i) = x1(i)
2270  alpha(2, i) = x2(i) - x1(i)
2271  alpha(3, i) = x4(i) - x1(i)
2272  alpha(4, i) = x5(i) - x1(i)
2273  alpha(5, i) = x3(i) - x2(i) + x1(i) - x4(i)
2274  alpha(6, i) = x6(i) - x5(i) + x1(i) - x2(i)
2275  alpha(7, i) = x8(i) - x5(i) + x1(i) - x4(i)
2276  alpha(8, i) = x7(i) - x8(i) + x5(i) - x6(i) + x4(i) - x3(i) + x2(i) - x1(i)
2277  end do
2278 
2279  ! Compute the value of the partial derivatives of x, y, and z with
2280  ! respect to ksi, eta, and zeta
2281 
2282  dxdksi = alpha(2, 1) + alpha(5, 1) * eta + alpha(6, 1) * zeta + alpha(8, 1) * eta * zeta
2283  dydksi = alpha(2, 2) + alpha(5, 2) * eta + alpha(6, 2) * zeta + alpha(8, 2) * eta * zeta
2284  dzdksi = alpha(2, 3) + alpha(5, 3) * eta + alpha(6, 3) * zeta + alpha(8, 3) * eta * zeta
2285 
2286  dxdeta = alpha(3, 1) + alpha(5, 1) * ksi + alpha(7, 1) * zeta + alpha(8, 1) * ksi * zeta
2287  dydeta = alpha(3, 2) + alpha(5, 2) * ksi + alpha(7, 2) * zeta + alpha(8, 2) * ksi * zeta
2288  dzdeta = alpha(3, 3) + alpha(5, 3) * ksi + alpha(7, 3) * zeta + alpha(8, 3) * ksi * zeta
2289 
2290  dxdzeta = alpha(4, 1) + alpha(6, 1) * ksi + alpha(7, 1) * eta + alpha(8, 1) * ksi * eta
2291  dydzeta = alpha(4, 2) + alpha(6, 2) * ksi + alpha(7, 2) * eta + alpha(8, 2) * ksi * eta
2292  dzdzeta = alpha(4, 3) + alpha(6, 3) * ksi + alpha(7, 3) * eta + alpha(8, 3) * ksi * eta
2293 
2294  ! Compute the value of the non-zero partial derivatives of x, y, and z
2295  ! with respect to ksi, eta, and zeta. Note that the second derivatives
2296  ! of x, y, and z with respect to ksi2, eta2, and zeta2 are identically
2297  ! zero
2298 
2299  d2xdksideta = alpha(5, 1) + alpha(8, 1) * zeta
2300  d2ydksideta = alpha(5, 2) + alpha(8, 2) * zeta
2301  d2zdksideta = alpha(5, 3) + alpha(8, 3) * zeta
2302 
2303  d2xdksidzeta = alpha(6, 1) + alpha(8, 1) * eta
2304  d2ydksidzeta = alpha(6, 2) + alpha(8, 2) * eta
2305  d2zdksidzeta = alpha(6, 3) + alpha(8, 3) * eta
2306 
2307  d2xdetadzeta = alpha(7, 1) + alpha(8, 1) * ksi
2308  d2ydetadzeta = alpha(7, 2) + alpha(8, 2) * ksi
2309  d2zdetadzeta = alpha(7, 3) + alpha(8, 3) * ksi
2310 
2311  ! Compute the actual location of the point
2312 
2313  x0 = alpha(1, 1) + alpha(2, 1) * ksi + alpha(3, 1) * eta + alpha(4, 1) * zeta &
2314  + alpha(5, 1) * ksi * eta + alpha(6, 1) * ksi * zeta + alpha(7, 1) * eta * zeta &
2315  + alpha(8, 1) * ksi * eta * zeta
2316  y0 = alpha(1, 2) + alpha(2, 2) * ksi + alpha(3, 2) * eta + alpha(4, 2) * zeta &
2317  + alpha(5, 2) * ksi * eta + alpha(6, 2) * ksi * zeta + alpha(7, 2) * eta * zeta &
2318  + alpha(8, 2) * ksi * eta * zeta
2319  z0 = alpha(1, 3) + alpha(2, 3) * ksi + alpha(3, 3) * eta + alpha(4, 3) * zeta &
2320  + alpha(5, 3) * ksi * eta + alpha(6, 3) * ksi * zeta + alpha(7, 3) * eta * zeta &
2321  + alpha(8, 3) * ksi * eta * zeta
2322 
2323  ! Compute the elements of the Hessian matrix
2324 
2325  hess(1, 1) = 2.0_realtype * ((dxdksi * dxdksi) + (dydksi * dydksi) + (dzdksi * dzdksi))
2326  hess(2, 2) = 2.0_realtype * ((dxdeta * dxdeta) + (dydeta * dydeta) + (dzdeta * dzdeta))
2327  hess(3, 3) = 2.0_realtype * ((dxdzeta * dxdzeta) + (dydzeta * dydzeta) + (dzdzeta * dzdzeta))
2328 
2329  hess(1, 2) = 2.0_realtype * ((dxdksi * dxdeta) + (dydksi * dydeta) + (dzdksi * dzdeta) &
2330  - ((xp(1) - x0) * d2xdksideta) - &
2331  ((xp(2) - y0) * d2ydksideta) - ((xp(3) - z0) * d2zdksideta))
2332  hess(1, 3) = 2.0_realtype * ((dxdksi * dxdzeta) + (dydksi * dydzeta) + (dzdksi * dzdzeta) &
2333  - ((xp(1) - x0) * d2xdksidzeta) - &
2334  ((xp(2) - y0) * d2ydksidzeta) - ((xp(3) - z0) * d2zdksidzeta))
2335  hess(2, 3) = 2.0_realtype * ((dxdeta * dxdzeta) + (dydeta * dydzeta) + (dzdeta * dzdzeta) &
2336  - ((xp(1) - x0) * d2xdetadzeta) - &
2337  ((xp(2) - y0) * d2ydetadzeta) - ((xp(3) - z0) * d2zdetadzeta))
2338 
2339  hess(2, 1) = hess(1, 2)
2340  hess(3, 1) = hess(1, 3)
2341  hess(3, 2) = hess(2, 3)
2342 
2343  ! Return iErr = 0 for the time being.
2344 
2345  ierr = 0
2346 
2347  return
2348 
2349  end subroutine hessd2hexa
2350 
2351  subroutine gradd2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, chi, x0, y0, z0, grad, iErr)
2352  !
2353  ! Compute the gradient of the square of the distance between the
2354  ! point xP, and the actual point in the hexahedron represented
2355  ! by chi(1:3).
2356  !
2357  use precision
2358 
2359  implicit none
2360  !
2361  ! Subroutine arguments.
2362  !
2363  real(kind=realtype), intent(out) :: x0, y0, z0
2364  real(kind=realtype), dimension(3), intent(in) :: xp, chi
2365  real(kind=realtype), dimension(3), intent(in) :: x1, x2, x3, x4
2366  real(kind=realtype), dimension(3), intent(in) :: x5, x6, x7, x8
2367  real(kind=realtype), dimension(3), intent(out) :: grad
2368 
2369  integer(kind=intType), intent(out) :: iErr
2370 
2371  !
2372  ! Local variables.
2373  !
2374  integer(kind=intType) :: i
2375 
2376  real(kind=realtype) :: ksi, eta, zeta
2377  real(kind=realtype), dimension(8, 3) :: alpha
2378  real(kind=realtype) :: dxdksi, dxdeta, dxdzeta
2379  real(kind=realtype) :: dydksi, dydeta, dydzeta
2380  real(kind=realtype) :: dzdksi, dzdeta, dzdzeta
2381 
2382  ! Initialize the parametric coordinates with a more recognizable
2383  ! name.
2384 
2385  ksi = chi(1)
2386  eta = chi(2)
2387  zeta = chi(3)
2388 
2389  ! Initialize the alpha array (obtained by regrouping terms in the
2390  ! parametric expansion of the (ksi,eta,zeta)->(x,y,z) mapping.
2391  ! The second index of the alpha array corresponds to the x, y, or
2392  ! z coordinate mapping.
2393 
2394  do i = 1, 3
2395  alpha(1, i) = x1(i)
2396  alpha(2, i) = x2(i) - x1(i)
2397  alpha(3, i) = x4(i) - x1(i)
2398  alpha(4, i) = x5(i) - x1(i)
2399  alpha(5, i) = x3(i) - x2(i) + x1(i) - x4(i)
2400  alpha(6, i) = x6(i) - x5(i) + x1(i) - x2(i)
2401  alpha(7, i) = x8(i) - x5(i) + x1(i) - x4(i)
2402  alpha(8, i) = x7(i) - x8(i) + x5(i) - x6(i) + x4(i) - x3(i) + x2(i) - x1(i)
2403  end do
2404 
2405  ! Compute the value of the partial derivatives of x, y, and z with
2406  ! respect to ksi, eta, and zeta
2407 
2408  dxdksi = alpha(2, 1) + alpha(5, 1) * eta + alpha(6, 1) * zeta + alpha(8, 1) * eta * zeta
2409  dydksi = alpha(2, 2) + alpha(5, 2) * eta + alpha(6, 2) * zeta + alpha(8, 2) * eta * zeta
2410  dzdksi = alpha(2, 3) + alpha(5, 3) * eta + alpha(6, 3) * zeta + alpha(8, 3) * eta * zeta
2411 
2412  dxdeta = alpha(3, 1) + alpha(5, 1) * ksi + alpha(7, 1) * zeta + alpha(8, 1) * ksi * zeta
2413  dydeta = alpha(3, 2) + alpha(5, 2) * ksi + alpha(7, 2) * zeta + alpha(8, 2) * ksi * zeta
2414  dzdeta = alpha(3, 3) + alpha(5, 3) * ksi + alpha(7, 3) * zeta + alpha(8, 3) * ksi * zeta
2415 
2416  dxdzeta = alpha(4, 1) + alpha(6, 1) * ksi + alpha(7, 1) * eta + alpha(8, 1) * ksi * eta
2417  dydzeta = alpha(4, 2) + alpha(6, 2) * ksi + alpha(7, 2) * eta + alpha(8, 2) * ksi * eta
2418  dzdzeta = alpha(4, 3) + alpha(6, 3) * ksi + alpha(7, 3) * eta + alpha(8, 3) * ksi * eta
2419 
2420  ! Compute the actual location of the point
2421 
2422  x0 = alpha(1, 1) + alpha(2, 1) * ksi + alpha(3, 1) * eta + alpha(4, 1) * zeta &
2423  + alpha(5, 1) * ksi * eta + alpha(6, 1) * ksi * zeta + alpha(7, 1) * eta * zeta &
2424  + alpha(8, 1) * ksi * eta * zeta
2425  y0 = alpha(1, 2) + alpha(2, 2) * ksi + alpha(3, 2) * eta + alpha(4, 2) * zeta &
2426  + alpha(5, 2) * ksi * eta + alpha(6, 2) * ksi * zeta + alpha(7, 2) * eta * zeta &
2427  + alpha(8, 2) * ksi * eta * zeta
2428  z0 = alpha(1, 3) + alpha(2, 3) * ksi + alpha(3, 3) * eta + alpha(4, 3) * zeta &
2429  + alpha(5, 3) * ksi * eta + alpha(6, 3) * ksi * zeta + alpha(7, 3) * eta * zeta &
2430  + alpha(8, 3) * ksi * eta * zeta
2431 
2432  ! Compute the gradient components
2433 
2434  grad(1) = -2.0_realtype * (xp(1) - x0) * dxdksi &
2435  - 2.0_realtype * (xp(2) - y0) * dydksi &
2436  - 2.0_realtype * (xp(3) - z0) * dzdksi
2437  grad(2) = -2.0_realtype * (xp(1) - x0) * dxdeta &
2438  - 2.0_realtype * (xp(2) - y0) * dydeta &
2439  - 2.0_realtype * (xp(3) - z0) * dzdeta
2440  grad(3) = -2.0_realtype * (xp(1) - x0) * dxdzeta &
2441  - 2.0_realtype * (xp(2) - y0) * dydzeta &
2442  - 2.0_realtype * (xp(3) - z0) * dzdzeta
2443 
2444  ! Return iErr = 0 for the time being.
2445 
2446  ierr = 0
2447 
2448  return
2449 
2450  end subroutine gradd2hexa
2451 
2452 end module adtlocalsearch
adtlocalsearch::gradd2hexa
subroutine gradd2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, chi, x0, y0, z0, grad, iErr)
Definition: adtLocalSearch.F90:2352
adtlocalsearch::hessd2hexa
subroutine hessd2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, chi, hess, iErr)
Definition: adtLocalSearch.F90:2224
adtlocalsearch::containmenttreesearch
subroutine containmenttreesearch(ADT, coor, intInfo, uvw, arrDonor, nCoor, nInterpol)
Definition: adtLocalSearch.F90:23
adtlocalsearch::intersectiontreesearchsinglepoint
subroutine intersectiontreesearchsinglepoint(ADT, coor, intInfo, BB, frontLeaves, frontLeavesNew)
Definition: adtLocalSearch.F90:1725
adtlocalsearch::mindistancetreesearch
subroutine mindistancetreesearch(ADT, coor, intInfo, uvw, arrDonor, nCoor, nInterpol)
Definition: adtLocalSearch.F90:849
adtlocalsearch::containmenttreesearchsinglepoint
subroutine containmenttreesearchsinglepoint(ADT, coor, intInfo, uvw, arrDonor, nInterpol, BB, frontLeaves, frontLeavesNew, failed)
Definition: adtLocalSearch.F90:115
adtlocalsearch::newtonstep
subroutine newtonstep(hess, grad, step, iErr)
Definition: adtLocalSearch.F90:2162
adtlocalsearch::mindistancetreesearchsinglepoint
subroutine mindistancetreesearchsinglepoint(ADT, coor, intInfo, uvw, arrDonor, nInterpol, BB, frontLeaves, frontLeavesNew)
Definition: adtLocalSearch.F90:959
adtlocalsearch::mind2hexa
subroutine mind2hexa(xP, x1, x2, x3, x4, x5, x6, x7, x8, d2, chi, iErr)
Definition: adtLocalSearch.F90:1901
adtutils::qsortbboxtargets
subroutine qsortbboxtargets(arr, nn, ADT)
Definition: adtUtils.F90:515
adtutils::reallocbboxtargettypeplus
subroutine reallocbboxtargettypeplus(arr, nSize, nInc, ADT)
Definition: adtUtils.F90:796
adtutils::reallocplus
subroutine reallocplus(arr, nSize, nInc, ADT)
Definition: adtUtils.F90:860
adtutils::nstack
integer(kind=inttype) nstack
Definition: adtUtils.F90:19
adtutils::stack
integer(kind=inttype), dimension(:), pointer stack
Definition: adtUtils.F90:20
adtutils::adtterminate
subroutine adtterminate(ADT, routineName, errorMessage)
Definition: adtUtils.F90:29
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::adttetrahedron
integer(kind=adtelementtype), parameter adttetrahedron
Definition: constants.F90:251
constants::adtprism
integer(kind=adtelementtype), parameter adtprism
Definition: constants.F90:253
constants::adthexahedron
integer(kind=adtelementtype), parameter adthexahedron
Definition: constants.F90:254
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::adttriangle
integer(kind=adtelementtype), parameter adttriangle
Definition: constants.F90:249
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::adtquadrilateral
integer(kind=adtelementtype), parameter adtquadrilateral
Definition: constants.F90:250
constants::fourth
real(kind=realtype), parameter fourth
Definition: constants.F90:82
constants::adtpyramid
integer(kind=adtelementtype), parameter adtpyramid
Definition: constants.F90:252
precision::realtype
integer, parameter realtype
Definition: precision.F90:112
vf
Definition: vf.F90:1
Generated by doxygen 1.9.1