MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
findptsedge_local_2.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "../gslib.hpp"
13#include "../../general/forall.hpp"
14#include "gslib_kernel_helpers.hpp"
15
16#ifdef MFEM_USE_GSLIB
17
18#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
19#pragma GCC diagnostic push
20#pragma GCC diagnostic ignored "-Wunused-function"
21#endif
22#include "gslib.h"
23#ifndef GSLIB_RELEASE_VERSION //gslib v1.0.7
24#define GSLIB_RELEASE_VERSION 10007
25#endif
26#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
27#pragma GCC diagnostic pop
28#endif
29
30namespace mfem
31{
32#if GSLIB_RELEASE_VERSION >= 10009
33#define CODE_INTERNAL 0
34#define CODE_BORDER 1
35#define CODE_NOT_FOUND 2
36#define sDIM 2
37#define sDIM2 4
38#define rDIM 1
39
40struct findptsElementPoint_t
41{
42 double x[sDIM], r, oldr, dist2, dist2p, tr;
43 int flags;
44};
45
46struct findptsElementGEdge_t
47{
48 double *x[sDIM];
49};
50
51struct findptsElementGPT_t
52{
53 double x[sDIM], jac[sDIM*rDIM], hes[sDIM*rDIM];
54};
55
56using dbl_range_t = gslib::dbl_range_t;
57using obbox_t = gslib::obbox_t<sDIM>;
58using findptsLocalHashData_t = gslib::findptsLocalHashData_t<sDIM>;
59using gslib::AABB_test;
60using gslib::bbox_test;
61using gslib::hash_index;
62using gslib::l2norm2;
63using gslib::lag_eval_second_der;
64
65/* the bit structure of flags is CRR
66 the C bit --- 1<<2 --- is set when the point is converged
67 RR is 0 = 00b if r is unconstrained,
68 1 = 01b if r is constrained at -1, i.e., rmin
69 2 = 10b if r is constrained at +1, i.e., rmax
70*/
71
72#define CONVERGED_FLAG (1u<<2)
73#define FLAG_MASK 0x07u // = 111b
74
75/* returns 1 if r direction (the only free direction in 2D) is constrained.
76 returns 1 if either 1st or 2nd bit of flags is set.
77*/
78static MFEM_HOST_DEVICE inline int num_constrained(const int flags)
79{
80 return ((flags | flags>>1) & 1u);
81}
82
83/* pi=0, r=-1; pi=1, r=+1 */
84static MFEM_HOST_DEVICE inline int point_index(const int x)
85{
86 return ((x>>1) & 1u);
87}
88
89/* check reduction in objective against prediction, and adjust
90 trust region radius (p->tr) accordingly;
91 may reject the prior step, returning 1; otherwise returns 0
92 sets out_pt->dist2, out_pt->index, out_pt->x, out_pt->oldr in any event,
93 leaving out_pt->r, out_pt->dr, out_pt->flags to be set when returning 0 */
94static MFEM_HOST_DEVICE bool reject_prior_step_q(findptsElementPoint_t *out_pt,
95 const double resid[2],
96 const findptsElementPoint_t *p,
97 const double tol)
98{
99 const double dist2 = l2norm2<2>(resid);
100 const double decr = p->dist2 - dist2;
101 const double pred = p->dist2p;
102 out_pt->x[0] = p->x[0];
103 out_pt->x[1] = p->x[1];
104 out_pt->oldr = p->r;
105 out_pt->dist2 = dist2;
106 if (decr >= 0.01*pred)
107 {
108 if (decr >= 0.9*pred) // very good iteration
109 {
110 out_pt->tr = p->tr*2;
111 }
112 else // somewhat good iteration
113 {
114 out_pt->tr = p->tr;
115 }
116 return false;
117 }
118 else
119 {
120 /* reject step; note: the point will pass through this routine
121 again, and we set things up here so it gets classed as a
122 "very good iteration" --- this doubles the trust radius,
123 which is why we divide by 4 below */
124 double v0 = fabs(p->r - p->oldr);
125 out_pt->tr = v0/4.0;
126 out_pt->dist2 = p->dist2;
127 out_pt->r = p->oldr;
128 out_pt->flags = p->flags>>3;
129 out_pt->dist2p = -HUGE_VAL;
130 if (pred < dist2*tol)
131 {
132 out_pt->flags |= CONVERGED_FLAG;
133 }
134 return true;
135 }
136}
137
138static MFEM_HOST_DEVICE inline void newton_edge( findptsElementPoint_t *const
139 out_pt,
140 const double jac[2],
141 const double rhess,
142 const double resid[2],
143 int flags,
144 const findptsElementPoint_t *const p,
145 const double tol )
146{
147 const double tr = p->tr;
148 const double A = jac[0] * jac[0] + jac[1] * jac[1] -
149 rhess; // A = J^T J - resid_d H_d
150 const double y = jac[0]*resid[0] + jac[1]*resid[1]; // y = J^T resid
151
152 const double oldr = p->r;
153 double dr, newr, tdr, tnewr, v, tv;
154 int new_flags=0, tnew_flags=0;
155
156#define EVAL(dr) ( (dr*A - 2*y) * dr )
157 if (A>0)
158 {
159 dr = y/A;
160 if (fabs(dr)<tol)
161 {
162 dr=0.0;
163 newr = oldr;
164 }
165 else
166 {
167 newr = oldr+dr;
168 }
169
170 if (fabs(dr)<tr && fabs(newr)<1)
171 {
172 v = EVAL(dr);
173 goto newton_edge_fin;
174 }
175 }
176
177 if ((newr=oldr-tr) > -1)
178 {
179 dr = -tr;
180 }
181 else
182 {
183 newr = -1, dr = -1-oldr, new_flags = flags|1u;
184 }
185 v = EVAL(dr);
186
187 if ((tnewr=oldr+tr) < 1)
188 {
189 tdr = tr;
190 }
191 else
192 {
193 tnewr = 1, tdr = 1-oldr, tnew_flags = flags|2u;
194 }
195 tv = EVAL(tdr);
196
197 if (tv<v)
198 {
199 newr = tnewr, dr = tdr, v = tv, new_flags = tnew_flags;
200 }
201#undef EVAL
202
203newton_edge_fin:
204 // check convergence by testing if change in r is less than tol
205 if (fabs(dr)<tol)
206 {
207 new_flags |= CONVERGED_FLAG;
208 }
209 out_pt->r = newr;
210 out_pt->dist2p = -v;
211 out_pt->flags = flags | new_flags | ((p->flags & FLAG_MASK)<<3);
212}
213
214static MFEM_HOST_DEVICE void seed_j( const double *elx[sDIM],
215 const double x[sDIM],
216 const double *z,
217 double *dist2,
218 double *r,
219 const int ir,
220 const int pN )
221{
222 double dx[sDIM];
223 for (int d=0; d<sDIM; ++d)
224 {
225 dx[d] = x[d] - elx[d][ir];
226 }
227 dist2[ir] = HUGE_VAL;
228 const double dist2_rs = l2norm2(dx);
229 if (dist2[ir]>dist2_rs)
230 {
231 dist2[ir] = dist2_rs;
232 r[ir] = z[ir];
233 }
234}
235
236template<int T_D1D = 0>
237static void FindPointsEdgeLocal2DKernel( const int npt,
238 const double tol,
239 const double dist2tol,
240 const double *x,
241 const int point_pos_ordering,
242 const double *xElemCoord,
243 const int nel,
244 const double *wtend,
245 const double *boxinfo,
246 const bool obb_check,
247 const int hash_n,
248 const double *hashMin,
249 const double *hashFac,
250 unsigned int *hashOffset,
251 unsigned int *const code_base,
252 unsigned int *const el_base,
253 double *const r_base,
254 double *const dist2_base,
255 const double *gll1D,
256 const double *lagcoeff,
257 const int pN = 0 )
258{
259 const int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
260 const int D1D = T_D1D ? T_D1D : pN;
261 const int p_NEL = nel*D1D;
262 MFEM_VERIFY(MD1<=DofQuadLimits::MAX_D1D,
263 "Increase Max allowable polynomial order.");
264 MFEM_VERIFY(pN<=DofQuadLimits::MAX_D1D,
265 "Increase Max allowable polynomial order.");
266 MFEM_VERIFY(D1D!=0, "Polynomial order not specified.");
267 const int nThreads = D1D*sDIM;
268
269 mfem::forall_2D(npt, nThreads, 1, [=] MFEM_HOST_DEVICE (int i)
270 {
271 // 2D1D for seed, 3D1D + 7 for edge
272 constexpr int size1 = 3*MD1 + 7;
273 // edge coordinates = D1D*2
274 constexpr int size2 = 2*MD1;
275 // local element coordinates in shared memory
276 constexpr int size3 = MD1*sDIM;
277
278 MFEM_SHARED findptsElementPoint_t el_pts[2];
279 MFEM_SHARED double r_workspace[size1];
280
281 MFEM_SHARED double constraint_workspace[size2];
282
283 MFEM_SHARED double elem_coords[MD1 <= 6 ? size3 : 1];
284
285 double *r_workspace_ptr = r_workspace;
286 findptsElementPoint_t *fpt, *tmp;
287 fpt = el_pts + 0;
288 tmp = el_pts + 1;
289
290 // x and y coord index within point_pos for point i
291 int id_x = point_pos_ordering == 0 ? i : i*sDIM;
292 int id_y = point_pos_ordering == 0 ? i+npt : i*sDIM+1;
293 double x_i[2] = {x[id_x], x[id_y]};
294
295 unsigned int *code_i = code_base + i;
296 double *dist2_i = dist2_base + i;
297
298 //---------------- map_points_to_els --------------------
299 findptsLocalHashData_t hash;
300 for (int d=0; d<sDIM; ++d)
301 {
302 hash.bnd[d].min = hashMin[d];
303 hash.fac[d] = hashFac[d];
304 }
305 hash.hash_n = hash_n;
306 hash.offset = hashOffset;
307
308 const int hi = hash_index(&hash, x_i);
309 const unsigned int *elp = hash.offset + hash.offset[hi];
310 const unsigned int *const ele = hash.offset + hash.offset[hi+1];
311 *code_i = CODE_NOT_FOUND;
312 *dist2_i = HUGE_VAL;
313
314 for (; elp!=ele; ++elp)
315 {
316 const unsigned int el = *elp;
317
318 const int n_box_ents = obb_check ? (3*sDIM + sDIM2) : (2*sDIM);
319 bool pass_bb = true;
320 obbox_t box;
321 if (obb_check)
322 {
323 for (int idx = 0; idx < sDIM; ++idx)
324 {
325 box.c0[idx] = boxinfo[n_box_ents*el + idx];
326 box.x[idx].min = boxinfo[n_box_ents*el + sDIM + idx];
327 box.x[idx].max = boxinfo[n_box_ents*el + 2*sDIM + idx];
328 }
329 for (int idx = 0; idx < sDIM2; ++idx)
330 {
331 box.A[idx] = boxinfo[n_box_ents*el + 3*sDIM + idx];
332 }
333 pass_bb = (bbox_test(&box, x_i) >= 0);
334 }
335 else
336 {
337 for (int d = 0; d < sDIM; ++d)
338 {
339 box.x[d].min = boxinfo[n_box_ents*el + d];
340 box.x[d].max = boxinfo[n_box_ents*el + sDIM + d];
341 }
342 pass_bb = (AABB_test(&box, x_i) >= 0);
343 }
344
345 if (pass_bb)
346 {
347 //------------ findpts_local ------------------
348 {
349 if (MD1 <= 6)
350 {
351 MFEM_FOREACH_THREAD(j,x,D1D*sDIM)
352 {
353 const int qp = j % D1D;
354 const int d = j / D1D;
355 elem_coords[qp + d*D1D] =
356 xElemCoord[qp + el*D1D + d*p_NEL];
357 }
358 MFEM_SYNC_THREAD;
359 }
360
361 const double *elx[sDIM];
362 for (int d=0; d<sDIM; d++)
363 {
364 elx[d] = MD1<= 6 ? &elem_coords[d*D1D] :
365 xElemCoord + d*p_NEL + el*D1D;
366 }
367 MFEM_SYNC_THREAD;
368 //// findpts_el ////
369 {
370 MFEM_FOREACH_THREAD(j,x,1)
371 {
372 fpt->dist2 = HUGE_VAL;
373 fpt->dist2p = 0;
374 fpt->tr = 1;
375 }
376 MFEM_FOREACH_THREAD(j,x,sDIM)
377 {
378 fpt->x[j] = x_i[j];
379 }
380 MFEM_SYNC_THREAD;
381
382 {
383 double *dist2_temp = r_workspace_ptr;
384 double *r_temp = dist2_temp + D1D;
385 MFEM_FOREACH_THREAD(j,x,D1D)
386 {
387 seed_j(elx, x_i, gll1D, dist2_temp, r_temp, j, D1D);
388 }
389 MFEM_SYNC_THREAD;
390
391 MFEM_FOREACH_THREAD(j,x,1)
392 {
393 for (int ir=0; ir<D1D; ++ir)
394 {
395 if (dist2_temp[ir]<fpt->dist2)
396 {
397 fpt->dist2 = dist2_temp[ir];
398 fpt->r = r_temp[ir];
399 }
400 }
401 }
402 MFEM_SYNC_THREAD;
403 } //seed done
404
405 // Initialize tmp struct with fpt values before starting Newton iterations
406 MFEM_FOREACH_THREAD(j,x,1)
407 {
408 tmp->dist2 = HUGE_VAL;
409 tmp->dist2p = 0;
410 tmp->tr = 1;
411 tmp->flags = 0;
412 tmp->r = fpt->r;
413 }
414 MFEM_FOREACH_THREAD(j,x,sDIM)
415 {
416 tmp->x[j] = fpt->x[j];
417 }
418 MFEM_SYNC_THREAD;
419
420
421 for (int step=0; step<50; step++)
422 {
423 int nc = num_constrained(tmp->flags & FLAG_MASK);
424 switch (nc)
425 {
426 case 0:
427 {
428 double *wt = r_workspace_ptr;
429 double *resid = wt + 3*D1D;
430 double *jac = resid + sDIM;
431 double *hess = jac + sDIM*rDIM;
432
433 findptsElementGEdge_t edge;
434 for (int d=0; d<sDIM; ++d)
435 {
436 edge.x[d] = constraint_workspace + d*D1D;
437 }
438 MFEM_FOREACH_THREAD(j,x,D1D)
439 {
440 for (int d=0; d<sDIM; ++d)
441 {
442 edge.x[d][j] = elx[d][j];
443 }
444 }
445 MFEM_SYNC_THREAD;
446
447 // compute basis function info upto 2nd derivative
448 MFEM_FOREACH_THREAD(j,x,D1D)
449 {
450 lag_eval_second_der(wt, tmp->r, j, gll1D,
451 lagcoeff, D1D);
452 }
453 MFEM_SYNC_THREAD;
454
455 MFEM_FOREACH_THREAD(j,x,sDIM)
456 {
457 resid[j] = tmp->x[j];
458 jac[j] = 0.0;
459 hess[j] = 0.0;
460 for (int k=0; k<D1D; ++k)
461 {
462 resid[j] -= wt[ k]*edge.x[j][k];
463 jac[j] += wt[D1D+k]*edge.x[j][k];
464 hess[j] += wt[2*D1D+k]*edge.x[j][k];
465 }
466 }
467 MFEM_SYNC_THREAD;
468
469 MFEM_FOREACH_THREAD(j,x,1)
470 {
471 hess[2] = resid[0]*hess[0] + resid[1]*hess[1];
472 }
473 MFEM_SYNC_THREAD;
474
475 MFEM_FOREACH_THREAD(j,x,1)
476 {
477 if (!reject_prior_step_q(fpt, resid, tmp, tol))
478 {
479 newton_edge(fpt, jac, hess[2], resid,
480 tmp->flags & FLAG_MASK, tmp, tol);
481 }
482 }
483 MFEM_SYNC_THREAD;
484 break;
485 }
486 case 1: // r is constrained to either -1 or 1
487 {
488 MFEM_FOREACH_THREAD(j,x,1)
489 {
490 const int pi = point_index(tmp->flags &
491 FLAG_MASK);
492 const double *wt = wtend + pi*3*D1D;
493 findptsElementGPT_t gpt;
494 for (int d=0; d<sDIM; ++d)
495 {
496 gpt.x[d] = elx[d][pi*(D1D-1)];
497 gpt.jac[d] = 0.0;
498 gpt.hes[d] = 0.0;
499 for (int k=0; k<D1D; ++k)
500 {
501 gpt.jac[d] += wt[D1D +k]*elx[d][k];
502 gpt.hes[d] += wt[2*D1D+k]*elx[d][k];
503 }
504 }
505
506 const double *const pt_x = gpt.x;
507 const double *const jac = gpt.jac;
508 const double *const hes = gpt.hes;
509 double resid[sDIM], steep, sr;
510 resid[0] = fpt->x[0] - pt_x[0];
511 resid[1] = fpt->x[1] - pt_x[1];
512 steep = jac[0]*resid[0] + jac[1]*resid[1];
513 sr = steep*tmp->r;
514 if ( !reject_prior_step_q(fpt, resid, tmp, tol) )
515 {
516 if (sr<0)
517 {
518 const double rhess = resid[0]*hes[0] +
519 resid[1]*hes[1];
520 newton_edge(fpt, jac, rhess,
521 resid, 0, tmp, tol);
522 }
523 else // sr==0
524 {
525 fpt->r = tmp->r;
526 fpt->dist2p = 0;
527 fpt->flags = tmp->flags | CONVERGED_FLAG;
528 }
529 }
530 }
531 MFEM_SYNC_THREAD;
532 break;
533 } // case 1
534 } //switch
535 if (fpt->flags & CONVERGED_FLAG)
536 {
537 break;
538 }
539 MFEM_SYNC_THREAD;
540 MFEM_FOREACH_THREAD(j,x,1)
541 {
542 *tmp = *fpt;
543 }
544 MFEM_SYNC_THREAD;
545 } //for int step<50
546 } //findpts_el
547
548 bool converged_internal =
549 ((fpt->flags&FLAG_MASK) == CONVERGED_FLAG) &&
550 (fpt->dist2<dist2tol);
551
552 if (*code_i == CODE_NOT_FOUND || converged_internal ||
553 fpt->dist2 < *dist2_i)
554 {
555 MFEM_FOREACH_THREAD(j,x,1)
556 {
557 *(el_base+i) = el;
558 *code_i = converged_internal ? CODE_INTERNAL : CODE_BORDER;
559 *dist2_i = fpt->dist2;
560 *(r_base+i) = fpt->r;
561 }
562 MFEM_SYNC_THREAD;
563 if (converged_internal)
564 {
565 break;
566 }
567 }
568 } //findpts_local
569 } //obbox_test
570 } //elp
571 });
572}
573
574void FindPointsGSLIB::FindPointsEdgeLocal2( const Vector &point_pos,
575 int point_pos_ordering,
576 Array<unsigned int> &code,
577 Array<unsigned int> &elem,
578 Vector &ref,
579 Vector &dist,
580 int npt )
581{
582 if (npt==0)
583 {
584 return;
585 }
586 MFEM_VERIFY(dim==1 && spacedim==2,"Function for 2D edges only");
587 bool use_dev = point_pos.UseDevice();
588 auto pp = point_pos.Read(use_dev);
589 auto pgslm = gsl_mesh.Read(use_dev);
590 auto pwt = DEV.wtend.Read(use_dev);
591 auto pbb = DEV.bb.Read(use_dev);
592 auto plhm = DEV.lh_min.Read(use_dev);
593 auto plhf = DEV.lh_fac.Read(use_dev);
594 auto plho = DEV.lh_offset.ReadWrite(use_dev);
595 auto pcode = code.Write(use_dev);
596 auto pelem = elem.Write(use_dev);
597 auto pref = ref.Write(use_dev);
598 auto pdist = dist.Write(use_dev);
599 auto pgll1d = DEV.gll1d.ReadWrite(use_dev);
600 auto plc = DEV.lagcoeff.Read(use_dev);
601 double dist2tol = DEV.surf_dist_tol;
602 const bool obb_chk = obb_check;
603 switch (DEV.dof1d)
604 {
605 case 2:
606 FindPointsEdgeLocal2DKernel<2>(npt, DEV.newt_tol, dist2tol,
607 pp, point_pos_ordering, pgslm,
608 NE_split_total, pwt, pbb, obb_chk,
609 DEV.lh_nx, plhm, plhf, plho,
610 pcode, pelem, pref, pdist,
611 pgll1d, plc);
612 break;
613 case 3:
614 FindPointsEdgeLocal2DKernel<3>(npt, DEV.newt_tol, dist2tol,
615 pp, point_pos_ordering, pgslm,
616 NE_split_total, pwt, pbb, obb_chk,
617 DEV.lh_nx, plhm, plhf, plho,
618 pcode, pelem, pref, pdist,
619 pgll1d, plc);
620 break;
621 case 4:
622 FindPointsEdgeLocal2DKernel<4>(npt, DEV.newt_tol, dist2tol,
623 pp, point_pos_ordering, pgslm,
624 NE_split_total, pwt, pbb, obb_chk,
625 DEV.lh_nx, plhm, plhf, plho,
626 pcode, pelem, pref, pdist,
627 pgll1d, plc);
628 break;
629 default:
630 FindPointsEdgeLocal2DKernel(npt, DEV.newt_tol, dist2tol, pp,
631 point_pos_ordering, pgslm,
632 NE_split_total, pwt, pbb, obb_chk,
633 DEV.lh_nx, plhm, plhf, plho,
634 pcode, pelem, pref, pdist,
635 pgll1d, plc, DEV.dof1d);
636 break;
637 }
638}
639#undef sDIM
640#undef rDIM
641#undef sDIM2
642#undef CODE_INTERNAL
643#undef CODE_BORDER
644#undef CODE_NOT_FOUND
645#else
646void FindPointsGSLIB::FindPointsEdgeLocal2( const Vector &point_pos,
647 int point_pos_ordering,
648 Array<unsigned int> &code,
649 Array<unsigned int> &elem,
650 Vector &ref,
651 Vector &dist,
652 int npt ) {} ;
653#endif
654} // namespace mfem
655
656#endif //ifdef MFEM_USE_GSLIB
mfem::Array::ReadWrite
T * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:426
mfem::Array::Write
T * Write(bool on_dev=true)
Shortcut for mfem::Write(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:418
mfem::FindPointsGSLIB::FindPointsEdgeLocal2
void FindPointsEdgeLocal2(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 2D edge elements.
Definition findptsedge_local_2.cpp:574
mfem::FindPointsGSLIB::DEV
struct mfem::FindPointsGSLIB::DevStruct DEV
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Read
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:520
mfem::Vector::ReadWrite
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:536
mfem::Vector::UseDevice
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:145
mfem::Vector::Write
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:528
gslib_kernel_helpers.hpp
mfem::future::tr
MFEM_HOST_DEVICE T tr(const tensor< T, n, n > &A)
Returns the trace of a square matrix.
Definition tensor.hpp:1317
mfem::gslib::bbox_test
MFEM_HOST_DEVICE double bbox_test(const obbox_t< SDIM > *const b, const double(&x)[SDIM])
Definition gslib_kernel_helpers.hpp:128
mfem::gslib::l2norm2
MFEM_HOST_DEVICE double l2norm2(const double(&x)[SDIM])
Definition gslib_kernel_helpers.hpp:176
mfem::gslib::hash_index
MFEM_HOST_DEVICE int hash_index(const findptsLocalHashData_t< SDIM > *const p, const double(&x)[SDIM])
Definition gslib_kernel_helpers.hpp:159
mfem::gslib::lag_eval_second_der
MFEM_HOST_DEVICE void lag_eval_second_der(double *p0, double x, int i, const double *z, const double *lCoeff, int pN)
Definition gslib_kernel_helpers.hpp:81
mfem::gslib::AABB_test
MFEM_HOST_DEVICE double AABB_test(const obbox_t< SDIM > *const b, const double(&x)[SDIM])
Definition gslib_kernel_helpers.hpp:114
mfem::obbox_t
gslib::obbox_t< DIM > obbox_t
Definition findpts_local_2.cpp:56
mfem::findptsLocalHashData_t
gslib::findptsLocalHashData_t< DIM > findptsLocalHashData_t
Definition findpts_local_2.cpp:57
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:1220
mfem::dbl_range_t
gslib::dbl_range_t dbl_range_t
Definition findpts_local_2.cpp:55
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
mfem::FindPointsGSLIB::DevStruct::lh_offset
Array< unsigned int > lh_offset
Definition gslib.hpp:172
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0