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