19#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
20#pragma GCC diagnostic push
21#pragma GCC diagnostic ignored "-Wunused-function"
24#ifndef GSLIB_RELEASE_VERSION
25#define GSLIB_RELEASE_VERSION 10007
27#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
28#pragma GCC diagnostic pop
33#if GSLIB_RELEASE_VERSION >= 10009
34#define CODE_INTERNAL 0
36#define CODE_NOT_FOUND 2
39#define sDIM2 (sDIM*sDIM)
40#define rDIM2 (rDIM*rDIM)
42struct findptsElementPoint_t
44 double x[sDIM], r, oldr, dist2, dist2p, tr;
48struct findptsElementGEdge_t
50 double *x[sDIM], *dxdn[sDIM], *d2xdn[sDIM];
53struct findptsElementGPT_t
55 double x[sDIM], jac[sDIM], hes[sDIM*(1+1)];
59using obbox_t = gslib::obbox_t<sDIM>;
73#define CONVERGED_FLAG (1u<<2)
74#define FLAG_MASK 0x07u
78static MFEM_HOST_DEVICE
inline int num_constrained(
const int flags)
80 return ((flags | flags>>1) & 1u);
83static MFEM_HOST_DEVICE
inline int point_index(
const int x)
93static MFEM_HOST_DEVICE
bool reject_prior_step_q(findptsElementPoint_t *out_pt,
94 const double resid[3],
95 const findptsElementPoint_t *
p,
98 const double dist2 = l2norm2<sDIM>(resid);
99 const double decr =
p->dist2 - dist2;
100 const double pred =
p->dist2p;
101 for (
int d=0; d<sDIM; ++d)
103 out_pt->x[d] =
p->x[d];
106 out_pt->dist2 = dist2;
111 out_pt->tr = 2*
p->tr;
125 double v0 = fabs(
p->r -
p->oldr);
127 out_pt->dist2 =
p->dist2;
129 out_pt->flags =
p->flags>>3;
130 out_pt->dist2p = -HUGE_VAL;
133 out_pt->flags |= CONVERGED_FLAG;
139static MFEM_HOST_DEVICE
inline void newton_edge(findptsElementPoint_t *
const
141 const double jac[sDIM*rDIM],
143 const double resid[sDIM],
145 const findptsElementPoint_t *
const p,
148 const double tr =
p->tr;
150 const double A = jac[0]*jac[0]+ jac[1] * jac[1] + jac[2] * jac[2]
153 const double y = jac[0]*resid[0] + jac[1]*resid[1] + jac[0+2]*resid[2];
155 const double oldr =
p->r;
156 double dr, nr, tdr, tnr;
158 int new_flags = 0, tnew_flags = 0;
160#define EVAL(dr) (dr*A - 2*y)*dr
176 if ( fabs(dr)<tr && fabs(nr)<1 )
179 goto newton_edge_fin;
183 if ( (nr=oldr-tr)>-1 )
189 nr = -1, dr = -1-oldr, new_flags = flags | 1u;
193 if ( (tnr = oldr+tr)<1 )
199 tnr = 1, tdr = 1-oldr, tnew_flags = flags | 2u;
205 nr = tnr, dr = tdr, v = tv, new_flags = tnew_flags;
212 new_flags |= CONVERGED_FLAG;
216 out_pt->flags = flags | new_flags | ((
p->flags & FLAG_MASK)<<3);
220static MFEM_HOST_DEVICE
void seed_j(
const double *elx[sDIM],
221 const double x[sDIM],
234 for (
int d=0; d<sDIM; ++d)
236 dx[d] = x[d] - elx[d][ir];
238 dist2[ir] = l2norm2(dx);
242template<
int T_D1D = 0>
243static void FindPointsEdgeLocal3DKernel(
const int npt,
245 const double dist2tol,
247 const int point_pos_ordering,
248 const double *xElemCoord,
251 const double *boxinfo,
252 const bool obb_check,
254 const double *hashMin,
255 const double *hashFac,
256 unsigned int *hashOffset,
257 unsigned int *
const code_base,
258 unsigned int *
const el_base,
259 double *
const r_base,
260 double *
const dist2_base,
262 const double *lagcoeff,
265 const int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
266 const int D1D = T_D1D ? T_D1D : pN;
267 const int p_NEL = nel*D1D;
268 MFEM_VERIFY(MD1<=DofQuadLimits::MAX_D1D,
269 "Increase Max allowable polynomial order.");
270 MFEM_VERIFY(pN<=DofQuadLimits::MAX_D1D,
271 "Increase Max allowable polynomial order.");
272 MFEM_VERIFY(D1D!=0,
"Polynomial order not specified.");
273 const int nThreads = D1D*sDIM;
277 constexpr int size1 = 3*MD1 + 13;
278 constexpr int size2 = 3*MD1;
279 constexpr int size3 = MD1*sDIM;
281 MFEM_SHARED findptsElementPoint_t el_pts[2];
282 MFEM_SHARED
double r_workspace[size1];
284 MFEM_SHARED
double constraint_workspace[size2];
286 MFEM_SHARED
double elem_coords[MD1 <= 6 ? size3 : 1];
288 double *r_workspace_ptr = r_workspace;
289 findptsElementPoint_t *fpt, *tmp;
293 int id_x = point_pos_ordering==0 ? i : i*sDIM;
294 int id_y = point_pos_ordering==0 ? npt+i : 1+i*sDIM;
295 int id_z = point_pos_ordering==0 ? 2*npt+i : 2+i*sDIM;
296 double x_i[3] = {x[id_x], x[id_y], x[id_z]};
298 unsigned int *code_i = code_base + i;
299 double *dist2_i = dist2_base + i;
303 for (
int d=0; d<sDIM; ++d)
305 hash.bnd[d].min = hashMin[d];
306 hash.fac[d] = hashFac[d];
308 hash.hash_n = hash_n;
309 hash.offset = hashOffset;
311 const unsigned int hi = hash_index(&hash, x_i);
312 const unsigned int *elp = hash.offset + hash.offset[hi];
313 const unsigned int *
const ele = hash.offset + hash.offset[hi+1];
314 *code_i = CODE_NOT_FOUND;
317 for (; elp!=ele; ++elp)
319 const unsigned int el = *elp;
321 const int n_box_ents = obb_check ? (3*sDIM + sDIM2) : (2*sDIM);
326 for (
int idx = 0; idx < sDIM; ++idx)
328 box.c0[idx] = boxinfo[n_box_ents*el + idx];
329 box.x[idx].min = boxinfo[n_box_ents*el + sDIM + idx];
330 box.x[idx].max = boxinfo[n_box_ents*el + 2*sDIM + idx];
332 for (
int idx = 0; idx < sDIM2; ++idx)
334 box.A[idx] = boxinfo[n_box_ents*el + 3*sDIM + idx];
336 pass_bb = (bbox_test(&box, x_i) >= 0);
340 for (
int d = 0; d < sDIM; ++d)
342 box.x[d].min = boxinfo[n_box_ents*el + d];
343 box.x[d].max = boxinfo[n_box_ents*el + sDIM + d];
345 pass_bb = (AABB_test(&box, x_i) >= 0);
354 MFEM_FOREACH_THREAD(j,x,D1D*sDIM)
356 const int qp = j % D1D;
357 const int d = j / D1D;
358 elem_coords[qp + d*D1D] =
359 xElemCoord[qp + el*D1D + d*p_NEL];
364 const double *elx[sDIM];
365 for (
int d=0; d<sDIM; d++)
367 elx[d] = MD1<= 6 ? &elem_coords[d*D1D] :
368 xElemCoord + d*p_NEL + el*D1D;
373 MFEM_FOREACH_THREAD(j,x,1)
375 fpt->dist2 = HUGE_VAL;
379 MFEM_FOREACH_THREAD(j,x,sDIM)
387 double *dist2_temp = r_workspace_ptr;
388 double *r_temp = dist2_temp + D1D;
389 MFEM_FOREACH_THREAD(j,x,nThreads)
391 seed_j(elx, x_i, gll1D, dist2_temp, r_temp, j, D1D);
395 MFEM_FOREACH_THREAD(j,x,1)
397 fpt->dist2 = HUGE_VAL;
398 for (
int ir=0; ir<D1D; ++ir)
400 if (dist2_temp[ir] < fpt->dist2)
402 fpt->dist2 = dist2_temp[ir];
410 MFEM_FOREACH_THREAD(j,x,1)
412 tmp->dist2 = HUGE_VAL;
418 MFEM_FOREACH_THREAD(j,x,sDIM)
420 tmp->x[j] = fpt->x[j];
424 for (
int step=0; step<50; step++)
426 switch (num_constrained(tmp->flags & FLAG_MASK))
430 double *wt = r_workspace_ptr;
431 double *resid = wt + 3*D1D;
432 double *jac = resid + sDIM;
433 double *hess = jac + sDIM*rDIM;
435 findptsElementGEdge_t edge;
436 for (
int d=0; d<sDIM; ++d)
438 edge.x[d] = constraint_workspace + d*D1D;
440 MFEM_FOREACH_THREAD(j,x,D1D)
442 for (
int d=0; d<sDIM; ++d)
444 edge.x[d][j] = elx[d][j];
449 MFEM_FOREACH_THREAD(j,x,D1D)
451 lag_eval_second_der(wt, tmp->r, j, gll1D,
456 MFEM_FOREACH_THREAD(j,x,sDIM)
458 resid[j] = tmp->x[j];
461 for (
int k=0; k<D1D; ++k)
463 resid[j] -= wt[ k]*edge.x[j][k];
464 jac[j] += wt[D1D+k]*edge.x[j][k];
465 hess[j] += wt[2*D1D+k]*edge.x[j][k];
470 MFEM_FOREACH_THREAD(j,x,1)
472 hess[3] = resid[0]*hess[0] + resid[1]*hess[1] +
476 MFEM_FOREACH_THREAD(l,x,1)
478 if (!reject_prior_step_q(fpt,resid,tmp,tol))
480 newton_edge(fpt,jac,hess[3],resid,
481 tmp->flags&FLAG_MASK,tmp,tol);
489 MFEM_FOREACH_THREAD(j,x,1)
491 const int pi = point_index(tmp->flags &
493 const double *wt = wtend + pi*3*D1D;
494 findptsElementGPT_t gpt;
495 for (
int d=0; d<sDIM; ++d)
497 gpt.x[d] = elx[d][pi*(D1D-1)];
500 for (
int k=0; k<D1D; ++k)
502 gpt.jac[d] += wt[D1D +k]*elx[d][k];
503 gpt.hes[d] += wt[2*D1D+k]*elx[d][k];
507 const double *
const pt_x = gpt.x;
508 const double *
const jac = gpt.jac;
509 const double *
const hes = gpt.hes;
510 double resid[sDIM], steep, sr;
511 resid[0] = fpt->x[0] - pt_x[0];
512 resid[1] = fpt->x[1] - pt_x[1];
513 resid[2] = fpt->x[2] - pt_x[2];
514 steep = jac[0]*resid[0] + jac[1]*resid[1] +
517 if (!reject_prior_step_q(fpt, resid, tmp, tol))
521 const double rhess = resid[0]*hes[0] +
524 newton_edge(fpt, jac, rhess,
531 fpt->flags = tmp->flags | CONVERGED_FLAG;
539 if (fpt->flags & CONVERGED_FLAG)
545 MFEM_FOREACH_THREAD(j,x,1)
553 bool converged_internal =
554 ((fpt->flags&FLAG_MASK) == CONVERGED_FLAG) &&
555 (fpt->dist2<dist2tol);
556 if (*code_i==CODE_NOT_FOUND || converged_internal ||
559 MFEM_FOREACH_THREAD(j,x,1)
562 *code_i = converged_internal?CODE_INTERNAL:CODE_BORDER;
563 *dist2_i = fpt->dist2;
564 *(r_base+i) = fpt->r;
567 if (converged_internal)
579 int point_pos_ordering,
590 MFEM_VERIFY(
spacedim==3 &&
dim == 1,
"Function for 3D edges only");
592 auto pp = point_pos.
Read(use_dev);
599 auto pcode = code.
Write(use_dev);
600 auto pelem = elem.
Write(use_dev);
601 auto pref = ref.
Write(use_dev);
602 auto pdist = dist.
Write(use_dev);
610 FindPointsEdgeLocal3DKernel<2>(npt,
DEV.
newt_tol, dist2tol,
611 pp, point_pos_ordering, pgslm,
614 pcode, pelem, pref, pdist,
618 FindPointsEdgeLocal3DKernel<3>(npt,
DEV.
newt_tol, dist2tol,
619 pp, point_pos_ordering, pgslm,
622 pcode, pelem, pref, pdist,
626 FindPointsEdgeLocal3DKernel<4>(npt,
DEV.
newt_tol, dist2tol,
627 pp, point_pos_ordering, pgslm,
630 pcode, pelem, pref, pdist,
634 FindPointsEdgeLocal3DKernel(npt,
DEV.
newt_tol, dist2tol, pp,
635 point_pos_ordering, pgslm,
638 pcode, pelem, pref, pdist,
652 int point_pos_ordering,
T * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), on_dev).
T * Write(bool on_dev=true)
Shortcut for mfem::Write(a.GetMemory(), a.Size(), on_dev).
void FindPointsEdgeLocal3(const Vector &point_pos, int point_pos_ordering, Array< unsigned int > &gsl_code_dev_l, Array< unsigned int > &gsl_elem_dev_l, Vector &gsl_ref_l, Vector &gsl_dist_l, int npt)
FindPoints locally on device for 3D edge elements.
struct mfem::FindPointsGSLIB::DevStruct DEV
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
MFEM_HOST_DEVICE T tr(const tensor< T, n, n > &A)
Returns the trace of a square matrix.
MFEM_HOST_DEVICE double bbox_test(const obbox_t< SDIM > *const b, const double(&x)[SDIM])
MFEM_HOST_DEVICE double l2norm2(const double(&x)[SDIM])
MFEM_HOST_DEVICE int hash_index(const findptsLocalHashData_t< SDIM > *const p, const double(&x)[SDIM])
MFEM_HOST_DEVICE void lag_eval_second_der(double *p0, double x, int i, const double *z, const double *lCoeff, int pN)
MFEM_HOST_DEVICE double AABB_test(const obbox_t< SDIM > *const b, const double(&x)[SDIM])
gslib::obbox_t< DIM > obbox_t
gslib::findptsLocalHashData_t< DIM > findptsLocalHashData_t
void forall_2D(int N, int X, int Y, lambda &&body)
gslib::dbl_range_t dbl_range_t
real_t p(const Vector &x, real_t t)
Array< unsigned int > lh_offset