1 // Copyright (c) 2010-2021, 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
14 #ifdef MFEM_USE_GSLIB
15
16 // Ignore warnings from the gslib header (GCC version)
17 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
18 #pragma GCC diagnostic push
19 #pragma GCC diagnostic ignored "-Wunused-function"
20 #endif
21
22 // External GSLIB header (the MFEM header is gslib.hpp)
23 namespace gslib
24 {
25 #include "gslib.h"
26 }
27
28 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
29 #pragma GCC diagnostic pop
30 #endif
31
32 namespace mfem
33 {
34
FindPointsGSLIB()35 FindPointsGSLIB::FindPointsGSLIB()
36 : mesh(NULL), meshsplit(NULL), ir_simplex(NULL),
37 fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
38 dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
39 avgtype(AvgType::ARITHMETIC)
40 {
41 gsl_comm = new gslib::comm;
42 cr = new gslib::crystal;
43 #ifdef MFEM_USE_MPI
44 int initialized;
45 MPI_Initialized(&initialized);
46 if (!initialized) { MPI_Init(NULL, NULL); }
47 MPI_Comm comm = MPI_COMM_WORLD;
48 comm_init(gsl_comm, comm);
49 #else
50 comm_init(gsl_comm, 0);
51 #endif
52 }
53
~FindPointsGSLIB()54 FindPointsGSLIB::~FindPointsGSLIB()
55 {
56 delete gsl_comm;
57 delete cr;
58 delete ir_simplex;
59 delete meshsplit;
60 }
61
62 #ifdef MFEM_USE_MPI
FindPointsGSLIB(MPI_Comm comm_)63 FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
64 : mesh(NULL), meshsplit(NULL), ir_simplex(NULL),
65 fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
66 dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
67 avgtype(AvgType::ARITHMETIC)
68 {
69 gsl_comm = new gslib::comm;
70 cr = new gslib::crystal;
71 comm_init(gsl_comm, comm_);
72 }
73 #endif
74
Setup(Mesh & m,const double bb_t,const double newt_tol,const int npt_max)75 void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
76 const int npt_max)
77 {
78 MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
79 MFEM_VERIFY(m.GetNumGeometries(m.Dimension()) == 1,
80 "Mixed meshes are not currently supported in FindPointsGSLIB.");
81
82 // call FreeData if FindPointsGSLIB::Setup has been called already
83 if (setupflag) { FreeData(); }
84
85 crystal_init(cr, gsl_comm);
86 mesh = &m;
87 dim = mesh->Dimension();
88 const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
89 unsigned dof1D = fe->GetOrder() + 1;
90 const int gt = fe->GetGeomType();
91
92 if (gt == Geometry::TRIANGLE || gt == Geometry::TETRAHEDRON ||
93 gt == Geometry::PRISM)
94 {
95 GetSimplexNodalCoordinates();
96 }
97 else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
98 {
99 GetQuadHexNodalCoordinates();
100 }
101 else
102 {
103 MFEM_ABORT("Element type not currently supported in FindPointsGSLIB.");
104 }
105
106 const int pts_cnt = gsl_mesh.Size()/dim,
107 NEtot = pts_cnt/(int)pow(dof1D, dim);
108
109 if (dim == 2)
110 {
111 unsigned nr[2] = { dof1D, dof1D };
112 unsigned mr[2] = { 2*dof1D, 2*dof1D };
113 double * const elx[2] = { &gsl_mesh(0), &gsl_mesh(pts_cnt) };
114 fdata2D = findpts_setup_2(gsl_comm, elx, nr, NEtot, mr, bb_t,
115 pts_cnt, pts_cnt, npt_max, newt_tol);
116 }
117 else
118 {
119 unsigned nr[3] = { dof1D, dof1D, dof1D };
120 unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
121 double * const elx[3] =
122 { &gsl_mesh(0), &gsl_mesh(pts_cnt), &gsl_mesh(2*pts_cnt) };
123 fdata3D = findpts_setup_3(gsl_comm, elx, nr, NEtot, mr, bb_t,
124 pts_cnt, pts_cnt, npt_max, newt_tol);
125 }
126 setupflag = true;
127 }
128
FindPoints(const Vector & point_pos)129 void FindPointsGSLIB::FindPoints(const Vector &point_pos)
130 {
131 MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
132 points_cnt = point_pos.Size() / dim;
133 gsl_code.SetSize(points_cnt);
134 gsl_proc.SetSize(points_cnt);
135 gsl_elem.SetSize(points_cnt);
136 gsl_ref.SetSize(points_cnt * dim);
137 gsl_dist.SetSize(points_cnt);
138
139 if (dim == 2)
140 {
141 const double *xv_base[2];
142 xv_base[0] = point_pos.GetData();
143 xv_base[1] = point_pos.GetData() + points_cnt;
144 unsigned xv_stride[2];
145 xv_stride[0] = sizeof(double);
146 xv_stride[1] = sizeof(double);
147 findpts_2(gsl_code.GetData(), sizeof(unsigned int),
148 gsl_proc.GetData(), sizeof(unsigned int),
149 gsl_elem.GetData(), sizeof(unsigned int),
150 gsl_ref.GetData(), sizeof(double) * dim,
151 gsl_dist.GetData(), sizeof(double),
152 xv_base, xv_stride, points_cnt, fdata2D);
153 }
154 else
155 {
156 const double *xv_base[3];
157 xv_base[0] = point_pos.GetData();
158 xv_base[1] = point_pos.GetData() + points_cnt;
159 xv_base[2] = point_pos.GetData() + 2*points_cnt;
160 unsigned xv_stride[3];
161 xv_stride[0] = sizeof(double);
162 xv_stride[1] = sizeof(double);
163 xv_stride[2] = sizeof(double);
164 findpts_3(gsl_code.GetData(), sizeof(unsigned int),
165 gsl_proc.GetData(), sizeof(unsigned int),
166 gsl_elem.GetData(), sizeof(unsigned int),
167 gsl_ref.GetData(), sizeof(double) * dim,
168 gsl_dist.GetData(), sizeof(double),
169 xv_base, xv_stride, points_cnt, fdata3D);
170 }
171
172 // Set the element number and reference position to 0 for points not found
173 for (int i = 0; i < points_cnt; i++)
174 {
175 if (gsl_code[i] == 2)
176 {
177 gsl_elem[i] = 0;
178 for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
179 }
180 }
181
182 // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for
183 // both simplices and quads.
184 MapRefPosAndElemIndices();
185 }
186
FindPoints(Mesh & m,const Vector & point_pos,const double bb_t,const double newt_tol,const int npt_max)187 void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
188 const double bb_t, const double newt_tol,
189 const int npt_max)
190 {
191 if (!setupflag || (mesh != &m) )
192 {
193 Setup(m, bb_t, newt_tol, npt_max);
194 }
195 FindPoints(point_pos);
196 }
197
Interpolate(const Vector & point_pos,const GridFunction & field_in,Vector & field_out)198 void FindPointsGSLIB::Interpolate(const Vector &point_pos,
199 const GridFunction &field_in, Vector &field_out)
200 {
201 FindPoints(point_pos);
202 Interpolate(field_in, field_out);
203 }
204
Interpolate(Mesh & m,const Vector & point_pos,const GridFunction & field_in,Vector & field_out)205 void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
206 const GridFunction &field_in, Vector &field_out)
207 {
208 FindPoints(m, point_pos);
209 Interpolate(field_in, field_out);
210 }
211
FreeData()212 void FindPointsGSLIB::FreeData()
213 {
214 if (!setupflag) { return; }
215 crystal_free(cr);
216 if (dim == 2)
217 {
218 findpts_free_2(fdata2D);
219 }
220 else
221 {
222 findpts_free_3(fdata3D);
223 }
224 gsl_code.DeleteAll();
225 gsl_proc.DeleteAll();
226 gsl_elem.DeleteAll();
227 gsl_mesh.Destroy();
228 gsl_ref.Destroy();
229 gsl_dist.Destroy();
230 setupflag = false;
231 }
232
GetNodeValues(const GridFunction & gf_in,Vector & node_vals)233 void FindPointsGSLIB::GetNodeValues(const GridFunction &gf_in,
234 Vector &node_vals)
235 {
236 MFEM_ASSERT(gf_in.FESpace()->GetVDim() == 1, "Scalar function expected.");
237
238 const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
239 const Geometry::Type gt = fe->GetGeomType();
240 const int NE = mesh->GetNE();
241
242 if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
243 {
244 const GridFunction *nodes = mesh->GetNodes();
245 const FiniteElementSpace *fes = nodes->FESpace();
246 const IntegrationRule &ir = fes->GetFE(0)->GetNodes();
247 const int dof_cnt = ir.GetNPoints();
248
249 node_vals.SetSize(NE * dof_cnt);
250
251 const TensorBasisElement *tbe =
252 dynamic_cast<const TensorBasisElement *>(fes->GetFE(0));
253 MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
254 const Array<int> &dof_map = tbe->GetDofMap();
255
256 int pt_id = 0;
257 Vector vals_el;
258 for (int i = 0; i < NE; i++)
259 {
260 gf_in.GetValues(i, ir, vals_el);
261 for (int j = 0; j < dof_cnt; j++)
262 {
263 node_vals(pt_id++) = vals_el(dof_map[j]);
264 }
265 }
266 }
267 else if (gt == Geometry::TRIANGLE || gt == Geometry::TETRAHEDRON ||
268 gt == Geometry::PRISM)
269 {
270 const int dof_cnt = ir_simplex->GetNPoints();
271 node_vals.SetSize(NE * dof_cnt);
272
273 int pt_id = 0;
274 Vector vals_el;
275 for (int j = 0; j < NE; j++)
276 {
277 gf_in.GetValues(j, *ir_simplex, vals_el);
278 for (int i = 0; i < dof_cnt; i++)
279 {
280 node_vals(pt_id++) = vals_el(i);
281 }
282 }
283 }
284 else
285 {
286 MFEM_ABORT("Element type not currently supported.");
287 }
288 }
289
GetQuadHexNodalCoordinates()290 void FindPointsGSLIB::GetQuadHexNodalCoordinates()
291 {
292 const GridFunction *nodes = mesh->GetNodes();
293 const FiniteElementSpace *fes = nodes->FESpace();
294
295 const int NE = mesh->GetNE(),
296 dof_cnt = fes->GetFE(0)->GetDof(),
297 pts_cnt = NE * dof_cnt;
298 gsl_mesh.SetSize(dim * pts_cnt);
299
300 const TensorBasisElement *tbe =
301 dynamic_cast<const TensorBasisElement *>(fes->GetFE(0));
302 MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
303 Array<int> dof_map(dof_cnt);
304 const Array<int> &dm = tbe->GetDofMap();
305 if (dm.Size() > 0) { dof_map = dm; }
306 else { for (int i = 0; i < dof_cnt; i++) { dof_map[i] = i; } }
307
308 DenseMatrix pos(dof_cnt, dim);
309 Vector posV(pos.Data(), dof_cnt * dim);
310 Array<int> xdofs(dof_cnt * dim);
311
312 int pt_id = 0;
313 for (int i = 0; i < NE; i++)
314 {
315 fes->GetElementVDofs(i, xdofs);
316 nodes->GetSubVector(xdofs, posV);
317 for (int j = 0; j < dof_cnt; j++)
318 {
319 for (int d = 0; d < dim; d++)
320 {
321 gsl_mesh(pts_cnt * d + pt_id) = pos(dof_map[j], d);
322 }
323 pt_id++;
324 }
325 }
326 }
327
GetSimplexNodalCoordinates()328 void FindPointsGSLIB::GetSimplexNodalCoordinates()
329 {
330 const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
331 const Geometry::Type gt = fe->GetGeomType();
332 const GridFunction *nodes = mesh->GetNodes();
333 const int NE = mesh->GetNE();
334 int NEsplit = 0;
335
336 // Split the reference element into a reference submesh of quads or hexes.
337 if (gt == Geometry::TRIANGLE)
338 {
339 int Nvert = 7;
340 NEsplit = 3;
341 meshsplit = new Mesh(2, Nvert, NEsplit, 0, 2);
342
343 const double quad_v[7][2] =
344 {
345 {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
346 {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
347 };
348 const int quad_e[3][4] =
349 {
350 {3, 4, 1, 0}, {4, 5, 2, 1}, {6, 5, 4, 3}
351 };
352
353 for (int j = 0; j < Nvert; j++)
354 {
355 meshsplit->AddVertex(quad_v[j]);
356 }
357 for (int j = 0; j < NEsplit; j++)
358 {
359 int attribute = j + 1;
360 meshsplit->AddQuad(quad_e[j], attribute);
361 }
362 meshsplit->FinalizeQuadMesh(1, 1, true);
363 }
364 else if (gt == Geometry::TETRAHEDRON)
365 {
366 int Nvert = 15;
367 NEsplit = 4;
368 meshsplit = new Mesh(3, Nvert, NEsplit, 0, 3);
369
370 const double hex_v[15][3] =
371 {
372 {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
373 {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
374 {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
375 {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
376 {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
377 };
378 const int hex_e[4][8] =
379 {
380 {0, 4, 10, 7, 6, 13, 14, 12},
381 {4, 1, 8, 10, 13, 5, 11, 14},
382 {13, 5, 11, 14, 6, 2, 9, 12},
383 {10, 8, 3, 7, 14, 11, 9, 12}
384 };
385
386 for (int j = 0; j < Nvert; j++)
387 {
388 meshsplit->AddVertex(hex_v[j]);
389 }
390 for (int j = 0; j < NEsplit; j++)
391 {
392 int attribute = j + 1;
393 meshsplit->AddHex(hex_e[j], attribute);
394 }
395 meshsplit->FinalizeHexMesh(1, 1, true);
396 }
397 else if (gt == Geometry::PRISM)
398 {
399 int Nvert = 14;
400 NEsplit = 3;
401 meshsplit = new Mesh(3, Nvert, NEsplit, 0, 3);
402
403 const double hex_v[14][3] =
404 {
405 {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
406 {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
407 {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
408 {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
409 };
410 const int hex_e[3][8] =
411 {
412 {3, 4, 1, 0, 10, 11, 8, 7},
413 {4, 5, 2, 1, 11, 12, 9, 8},
414 {6, 5, 4, 3, 13, 12, 11, 10}
415 };
416
417 for (int j = 0; j < Nvert; j++)
418 {
419 meshsplit->AddVertex(hex_v[j]);
420 }
421 for (int j = 0; j < NEsplit; j++)
422 {
423 int attribute = j + 1;
424 meshsplit->AddHex(hex_e[j], attribute);
425 }
426 meshsplit->FinalizeHexMesh(1, 1, true);
427 }
428 else { MFEM_ABORT("Unsupported geometry type."); }
429
430 // Curve the reference submesh.
431 H1_FECollection fec(fe->GetOrder(), dim);
432 FiniteElementSpace nodal_fes(meshsplit, &fec, dim);
433 meshsplit->SetNodalFESpace(&nodal_fes);
434
435 const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
436 pts_cnt = NEsplit * dof_cnt;
437 Vector irlist(dim * pts_cnt);
438
439 const TensorBasisElement *tbe =
440 dynamic_cast<const TensorBasisElement *>(nodal_fes.GetFE(0));
441 MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
442 const Array<int> &dof_map = tbe->GetDofMap();
443
444 DenseMatrix pos(dof_cnt, dim);
445 Vector posV(pos.Data(), dof_cnt * dim);
446 Array<int> xdofs(dof_cnt * dim);
447
448 // Create an IntegrationRule on the nodes of the reference submesh.
449 ir_simplex = new IntegrationRule(pts_cnt);
450 GridFunction *nodesplit = meshsplit->GetNodes();
451 int pt_id = 0;
452 for (int i = 0; i < NEsplit; i++)
453 {
454 nodal_fes.GetElementVDofs(i, xdofs);
455 nodesplit->GetSubVector(xdofs, posV);
456 for (int j = 0; j < dof_cnt; j++)
457 {
458 for (int d = 0; d < dim; d++)
459 {
460 irlist(pts_cnt * d + pt_id) = pos(dof_map[j], d);
461 }
462 ir_simplex->IntPoint(pt_id).x = irlist(pt_id);
463 ir_simplex->IntPoint(pt_id).y = irlist(pts_cnt + pt_id);
464 if (dim == 3)
465 {
466 ir_simplex->IntPoint(pt_id).z = irlist(2*pts_cnt + pt_id);
467 }
468 pt_id++;
469 }
470 }
471
472 // Initialize gsl_mesh with the positions of the split physical elements.
473 pt_id = 0;
474 Vector locval(dim);
475 const int tot_pts_cnt = pts_cnt*NE;
476 gsl_mesh.SetSize(tot_pts_cnt*dim);
477 for (int j = 0; j < NE; j++)
478 {
479 for (int i = 0; i < dof_cnt*NEsplit; i++)
480 {
481 const IntegrationPoint &ip = ir_simplex->IntPoint(i);
482 nodes->GetVectorValue(j, ip, locval);
483 for (int d = 0; d < dim; d++)
484 {
485 gsl_mesh(tot_pts_cnt*d + pt_id) = locval(d);
486 }
487 pt_id++;
488 }
489 }
490 }
491
MapRefPosAndElemIndices()492 void FindPointsGSLIB::MapRefPosAndElemIndices()
493 {
494 gsl_mfem_ref = gsl_ref;
495 gsl_mfem_elem = gsl_elem;
496 const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
497 const Geometry::Type gt = fe->GetGeomType();
498 int NEsplit = 0;
499
500 gsl_mfem_ref -= -1.; // map [-1, 1] to
501 gsl_mfem_ref *= 0.5; // [0, 1]
502 if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { return; }
503
504 H1_FECollection feclin(1, dim);
505 FiniteElementSpace nodal_fes_lin(meshsplit, &feclin, dim);
506 GridFunction gf_lin(&nodal_fes_lin);
507
508 if (gt == Geometry::TRIANGLE)
509 {
510 const double quad_v[7][2] =
511 {
512 {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
513 {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
514 };
515 for (int k = 0; k < dim; k++)
516 {
517 for (int j = 0; j < gf_lin.Size()/dim; j++)
518 {
519 gf_lin(j+k*gf_lin.Size()/dim) = quad_v[j][k];
520 }
521 }
522 NEsplit = 3;
523 }
524 else if (gt == Geometry::TETRAHEDRON)
525 {
526 const double hex_v[15][3] =
527 {
528 {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
529 {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
530 {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
531 {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
532 {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
533 };
534 for (int k = 0; k < dim; k++)
535 {
536 for (int j = 0; j < gf_lin.Size()/dim; j++)
537 {
538 gf_lin(j+k*gf_lin.Size()/dim) = hex_v[j][k];
539 }
540 }
541 NEsplit = 4;
542 }
543 else if (gt == Geometry::PRISM)
544 {
545 const double hex_v[14][3] =
546 {
547 {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
548 {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
549 {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
550 {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
551 };
552 for (int k = 0; k < dim; k++)
553 {
554 for (int j = 0; j < gf_lin.Size()/dim; j++)
555 {
556 gf_lin(j+k*gf_lin.Size()/dim) = hex_v[j][k];
557 }
558 }
559 NEsplit = 3;
560 }
561 else
562 {
563 MFEM_ABORT("Element type not currently supported.");
564 }
565
566 // Simplices are split into quads/hexes for GSLIB. For MFEM, we need to find
567 // the original element number and map the rst from micro to macro element.
568 for (int i = 0; i < points_cnt; i++)
569 {
570 if (gsl_code[i] == 2) { continue; }
571 int local_elem = gsl_elem[i]%NEsplit;
572 gsl_mfem_elem[i] = (gsl_elem[i] - local_elem)/NEsplit; // macro element number
573
574 IntegrationPoint ip;
575 Vector mfem_ref(gsl_mfem_ref.GetData()+i*dim, dim);
576 ip.Set2(mfem_ref.GetData());
577 if (dim == 3) { ip.z = mfem_ref(2); }
578 gf_lin.GetVectorValue(local_elem, ip, mfem_ref); // map to rst of macro element
579 }
580 }
581
Interpolate(const GridFunction & field_in,Vector & field_out)582 void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
583 Vector &field_out)
584 {
585 const int gf_order = field_in.FESpace()->GetFE(0)->GetOrder(),
586 mesh_order = mesh->GetNodalFESpace()->GetFE(0)->GetOrder();
587
588 const FiniteElementCollection *fec_in = field_in.FESpace()->FEColl();
589 const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>(fec_in);
590 const L2_FECollection *fec_l2 = dynamic_cast<const L2_FECollection *>(fec_in);
591
592 if (fec_h1 && gf_order == mesh_order &&
593 fec_h1->GetBasisType() == BasisType::GaussLobatto)
594 {
595 InterpolateH1(field_in, field_out);
596 return;
597 }
598 else
599 {
600 InterpolateGeneral(field_in, field_out);
601 if (!fec_l2 || avgtype == AvgType::NONE) { return; }
602 }
603
604 // For points on element borders, project the L2 GridFunction to H1 and
605 // re-interpolate.
606 if (fec_l2)
607 {
608 Array<int> indl2;
609 for (int i = 0; i < points_cnt; i++)
610 {
611 if (gsl_code[i] == 1) { indl2.Append(i); }
612 }
613 if (indl2.Size() == 0) { return; } // no points on element borders
614
615 Vector field_out_l2(field_out.Size());
616 VectorGridFunctionCoefficient field_in_dg(&field_in);
617 int gf_order_h1 = std::max(gf_order, 1); // H1 should be at least order 1
618 H1_FECollection fec(gf_order_h1, dim);
619 const int ncomp = field_in.FESpace()->GetVDim();
620 FiniteElementSpace fes(mesh, &fec, ncomp);
621 GridFunction field_in_h1(&fes);
622
623 if (avgtype == AvgType::ARITHMETIC)
624 {
625 field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
626 }
627 else if (avgtype == AvgType::HARMONIC)
628 {
629 field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
630 }
631 else
632 {
633 MFEM_ABORT("Invalid averaging type.");
634 }
635
636 if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
637 {
638 InterpolateH1(field_in_h1, field_out_l2);
639 }
640 else
641 {
642 InterpolateGeneral(field_in_h1, field_out_l2);
643 }
644
645 // Copy interpolated values for the points on element border
646 for (int j = 0; j < ncomp; j++)
647 {
648 for (int i = 0; i < indl2.Size(); i++)
649 {
650 int idx = indl2[i] + j*points_cnt;
651 field_out(idx) = field_out_l2(idx);
652 }
653 }
654 }
655 }
656
InterpolateH1(const GridFunction & field_in,Vector & field_out)657 void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
658 Vector &field_out)
659 {
660 FiniteElementSpace ind_fes(mesh, field_in.FESpace()->FEColl());
661 GridFunction field_in_scalar(&ind_fes);
662 Vector node_vals;
663
664 const int ncomp = field_in.FESpace()->GetVDim(),
665 points_fld = field_in.Size() / ncomp,
666 points_cnt = gsl_code.Size();
667
668 field_out.SetSize(points_cnt*ncomp);
669 field_out = default_interp_value;
670
671 for (int i = 0; i < ncomp; i++)
672 {
673 const int dataptrin = i*points_fld,
674 dataptrout = i*points_cnt;
675 field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin, points_fld);
676 GetNodeValues(field_in_scalar, node_vals);
677
678 if (dim==2)
679 {
680 findpts_eval_2(field_out.GetData()+dataptrout, sizeof(double),
681 gsl_code.GetData(), sizeof(unsigned int),
682 gsl_proc.GetData(), sizeof(unsigned int),
683 gsl_elem.GetData(), sizeof(unsigned int),
684 gsl_ref.GetData(), sizeof(double) * dim,
685 points_cnt, node_vals.GetData(), fdata2D);
686 }
687 else
688 {
689 findpts_eval_3(field_out.GetData()+dataptrout, sizeof(double),
690 gsl_code.GetData(), sizeof(unsigned int),
691 gsl_proc.GetData(), sizeof(unsigned int),
692 gsl_elem.GetData(), sizeof(unsigned int),
693 gsl_ref.GetData(), sizeof(double) * dim,
694 points_cnt, node_vals.GetData(), fdata3D);
695 }
696 }
697 }
698
InterpolateGeneral(const GridFunction & field_in,Vector & field_out)699 void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
700 Vector &field_out)
701 {
702 int ncomp = field_in.VectorDim(),
703 nptorig = points_cnt,
704 npt = points_cnt;
705
706 field_out.SetSize(points_cnt*ncomp);
707 field_out = default_interp_value;
708
709 if (gsl_comm->np == 1) // serial
710 {
711 for (int index = 0; index < npt; index++)
712 {
713 if (gsl_code[index] == 2) { continue; }
714 IntegrationPoint ip;
715 ip.Set2(gsl_mfem_ref.GetData()+index*dim);
716 if (dim == 3) { ip.z = gsl_mfem_ref(index*dim + 2); }
717 Vector localval(ncomp);
718 field_in.GetVectorValue(gsl_mfem_elem[index], ip, localval);
719 for (int i = 0; i < ncomp; i++)
720 {
721 field_out(index + i*npt) = localval(i);
722 }
723 }
724 }
725 else // parallel
726 {
727 // Determine number of points to be sent
728 int nptsend = 0;
729 for (int index = 0; index < npt; index++)
730 {
731 if (gsl_code[index] != 2) { nptsend +=1; }
732 }
733
734 // Pack data to send via crystal router
735 struct gslib::array *outpt = new gslib::array;
736 struct out_pt { double r[3], ival; uint index, el, proc; };
737 struct out_pt *pt;
738 array_init(struct out_pt, outpt, nptsend);
739 outpt->n=nptsend;
740 pt = (struct out_pt *)outpt->ptr;
741 for (int index = 0; index < npt; index++)
742 {
743 if (gsl_code[index] == 2) { continue; }
744 for (int d = 0; d < dim; ++d) { pt->r[d]= gsl_mfem_ref(index*dim + d); }
745 pt->index = index;
746 pt->proc = gsl_proc[index];
747 pt->el = gsl_mfem_elem[index];
748 ++pt;
749 }
750
751 // Transfer data to target MPI ranks
752 sarray_transfer(struct out_pt, outpt, proc, 1, cr);
753
754 if (ncomp == 1)
755 {
756 // Interpolate the grid function
757 npt = outpt->n;
758 pt = (struct out_pt *)outpt->ptr;
759 for (int index = 0; index < npt; index++)
760 {
761 IntegrationPoint ip;
762 ip.Set3(&pt->r[0]);
763 pt->ival = field_in.GetValue(pt->el, ip, 1);
764 ++pt;
765 }
766
767 // Transfer data back to source MPI rank
768 sarray_transfer(struct out_pt, outpt, proc, 1, cr);
769 npt = outpt->n;
770 pt = (struct out_pt *)outpt->ptr;
771 for (int index = 0; index < npt; index++)
772 {
773 field_out(pt->index) = pt->ival;
774 ++pt;
775 }
776 array_free(outpt);
777 delete outpt;
778 }
779 else // ncomp > 1
780 {
781 // Interpolate data and store in a Vector
782 npt = outpt->n;
783 pt = (struct out_pt *)outpt->ptr;
784 Vector vec_int_vals(npt*ncomp);
785 for (int index = 0; index < npt; index++)
786 {
787 IntegrationPoint ip;
788 ip.Set3(&pt->r[0]);
789 Vector localval(vec_int_vals.GetData()+index*ncomp, ncomp);
790 field_in.GetVectorValue(pt->el, ip, localval);
791 ++pt;
792 }
793
794 // Save index and proc data in a struct
795 struct gslib::array *savpt = new gslib::array;
796 struct sav_pt { uint index, proc; };
797 struct sav_pt *spt;
798 array_init(struct sav_pt, savpt, npt);
799 savpt->n=npt;
800 spt = (struct sav_pt *)savpt->ptr;
801 pt = (struct out_pt *)outpt->ptr;
802 for (int index = 0; index < npt; index++)
803 {
804 spt->index = pt->index;
805 spt->proc = pt->proc;
806 ++pt; ++spt;
807 }
808
809 array_free(outpt);
810 delete outpt;
811
812 // Copy data from save struct to send struct and send component wise
813 struct gslib::array *sendpt = new gslib::array;
814 struct send_pt { double ival; uint index, proc; };
815 struct send_pt *sdpt;
816 for (int j = 0; j < ncomp; j++)
817 {
818 array_init(struct send_pt, sendpt, npt);
819 sendpt->n=npt;
820 spt = (struct sav_pt *)savpt->ptr;
821 sdpt = (struct send_pt *)sendpt->ptr;
822 for (int index = 0; index < npt; index++)
823 {
824 sdpt->index = spt->index;
825 sdpt->proc = spt->proc;
826 sdpt->ival = vec_int_vals(j + index*ncomp);
827 ++sdpt; ++spt;
828 }
829
830 sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
831 sdpt = (struct send_pt *)sendpt->ptr;
832 for (int index = 0; index < nptorig; index++)
833 {
834 int idx = sdpt->index + j*nptorig;
835 field_out(idx) = sdpt->ival;
836 ++sdpt;
837 }
838 array_free(sendpt);
839 }
840 array_free(savpt);
841 delete sendpt;
842 delete savpt;
843 } // ncomp > 1
844 } // parallel
845 }
846
Setup(Mesh & m,const int meshid,GridFunction * gfmax,const double bb_t,const double newt_tol,const int npt_max)847 void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
848 GridFunction *gfmax,
849 const double bb_t, const double newt_tol,
850 const int npt_max)
851 {
852 MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
853 MFEM_VERIFY(m.GetNumGeometries(m.Dimension()) == 1,
854 "Mixed meshes are not currently supported in FindPointsGSLIB.");
855
856 // FreeData if OversetFindPointsGSLIB::Setup has been called already
857 if (setupflag) { FreeData(); }
858
859 crystal_init(cr, gsl_comm);
860 mesh = &m;
861 dim = mesh->Dimension();
862 const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
863 unsigned dof1D = fe->GetOrder() + 1;
864 int gt = fe->GetGeomType();
865
866 if (gt == Geometry::TRIANGLE || gt == Geometry::TETRAHEDRON ||
867 gt == Geometry::PRISM)
868 {
869 GetSimplexNodalCoordinates();
870 }
871 else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
872 {
873 GetQuadHexNodalCoordinates();
874 }
875 else
876 {
877 MFEM_ABORT("Element type not currently supported in FindPointsGSLIB.");
878 }
879
880 MFEM_ASSERT(meshid>=0, " The ID should be greater than or equal to 0.");
881
882 const int pts_cnt = gsl_mesh.Size()/dim,
883 NEtot = pts_cnt/(int)pow(dof1D, dim);
884
885 distfint.SetSize(pts_cnt);
886 if (!gfmax)
887 {
888 distfint = 0.;
889 }
890 else
891 {
892 GetNodeValues(*gfmax, distfint);
893 }
894 u_meshid = (unsigned int)meshid;
895
896 if (dim == 2)
897 {
898 unsigned nr[2] = { dof1D, dof1D };
899 unsigned mr[2] = { 2*dof1D, 2*dof1D };
900 double * const elx[2] = { &gsl_mesh(0), &gsl_mesh(pts_cnt) };
901 fdata2D = findptsms_setup_2(gsl_comm, elx, nr, NEtot, mr, bb_t,
902 pts_cnt, pts_cnt, npt_max, newt_tol,
903 &u_meshid, &distfint(0));
904 }
905 else
906 {
907 unsigned nr[3] = { dof1D, dof1D, dof1D };
908 unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
909 double * const elx[3] =
910 { &gsl_mesh(0), &gsl_mesh(pts_cnt), &gsl_mesh(2*pts_cnt) };
911 fdata3D = findptsms_setup_3(gsl_comm, elx, nr, NEtot, mr, bb_t,
912 pts_cnt, pts_cnt, npt_max, newt_tol,
913 &u_meshid, &distfint(0));
914 }
915 setupflag = true;
916 overset = true;
917 }
918
FindPoints(const Vector & point_pos,Array<unsigned int> & point_id)919 void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
920 Array<unsigned int> &point_id)
921 {
922 MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
923 "finding points.");
924 MFEM_VERIFY(overset, " Please setup FindPoints for overlapping grids.");
925 points_cnt = point_pos.Size() / dim;
926 unsigned int match = 0; // Don't find points in the mesh if point_id = mesh_id
927
928 gsl_code.SetSize(points_cnt);
929 gsl_proc.SetSize(points_cnt);
930 gsl_elem.SetSize(points_cnt);
931 gsl_ref.SetSize(points_cnt * dim);
932 gsl_dist.SetSize(points_cnt);
933
934 if (dim == 2)
935 {
936 const double *xv_base[2];
937 xv_base[0] = point_pos.GetData();
938 xv_base[1] = point_pos.GetData() + points_cnt;
939 unsigned xv_stride[2];
940 xv_stride[0] = sizeof(double);
941 xv_stride[1] = sizeof(double);
942 findptsms_2(gsl_code.GetData(), sizeof(unsigned int),
943 gsl_proc.GetData(), sizeof(unsigned int),
944 gsl_elem.GetData(), sizeof(unsigned int),
945 gsl_ref.GetData(), sizeof(double) * dim,
946 gsl_dist.GetData(), sizeof(double),
947 xv_base, xv_stride,
948 point_id.GetData(), sizeof(unsigned int), &match,
949 points_cnt, fdata2D);
950 }
951 else
952 {
953 const double *xv_base[3];
954 xv_base[0] = point_pos.GetData();
955 xv_base[1] = point_pos.GetData() + points_cnt;
956 xv_base[2] = point_pos.GetData() + 2*points_cnt;
957 unsigned xv_stride[3];
958 xv_stride[0] = sizeof(double);
959 xv_stride[1] = sizeof(double);
960 xv_stride[2] = sizeof(double);
961 findptsms_3(gsl_code.GetData(), sizeof(unsigned int),
962 gsl_proc.GetData(), sizeof(unsigned int),
963 gsl_elem.GetData(), sizeof(unsigned int),
964 gsl_ref.GetData(), sizeof(double) * dim,
965 gsl_dist.GetData(), sizeof(double),
966 xv_base, xv_stride,
967 point_id.GetData(), sizeof(unsigned int), &match,
968 points_cnt, fdata3D);
969 }
970
971 // Set the element number and reference position to 0 for points not found
972 for (int i = 0; i < points_cnt; i++)
973 {
974 if (gsl_code[i] == 2)
975 {
976 gsl_elem[i] = 0;
977 for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
978 }
979 }
980
981 // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for both
982 // simplices and quads.
983 MapRefPosAndElemIndices();
984 }
985
Interpolate(const Vector & point_pos,Array<unsigned int> & point_id,const GridFunction & field_in,Vector & field_out)986 void OversetFindPointsGSLIB::Interpolate(const Vector &point_pos,
987 Array<unsigned int> &point_id,
988 const GridFunction &field_in,
989 Vector &field_out)
990 {
991 FindPoints(point_pos, point_id);
992 Interpolate(field_in, field_out);
993 }
994
995
996 } // namespace mfem
997
998 #endif // MFEM_USE_GSLIB
999