MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
gslib.cpp
Go to the documentation of this file.
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 
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 
55 {
56  delete gsl_comm;
57  delete cr;
58  delete ir_simplex;
59  delete meshsplit;
60 }
61 
62 #ifdef MFEM_USE_MPI
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 
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  {
96  }
97  else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
98  {
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 
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;
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.
185 }
186 
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 
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 
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 
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  }
227  gsl_mesh.Destroy();
228  gsl_ref.Destroy();
229  gsl_dist.Destroy();
230  setupflag = false;
231 }
232 
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 
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 
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 
493 {
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 
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 &&
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 
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,
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 
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;
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 
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  {
870  }
871  else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
872  {
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 
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 
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,
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,
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.
984 }
985 
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
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe.hpp:243
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:582
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
Array< unsigned int > gsl_elem
Definition: gslib.hpp:60
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1324
virtual ~FindPointsGSLIB()
Definition: gslib.cpp:54
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:425
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5428
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:262
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2433
Array< unsigned int > gsl_mfem_elem
Definition: gslib.hpp:60
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:525
T * GetData()
Returns the data.
Definition: array.hpp:108
unsigned int uint
Definition: gecko.hpp:204
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:327
int VectorDim() const
Definition: gridfunc.cpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
Array< unsigned int > gsl_proc
Definition: gslib.hpp:60
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
Definition: mesh.cpp:2470
virtual void GetNodeValues(const GridFunction &gf_in, Vector &node_vals)
Get GridFunction from MFEM format to GSLIB format.
Definition: gslib.cpp:233
void DeleteAll()
Delete the whole array.
Definition: array.hpp:841
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:320
virtual void MapRefPosAndElemIndices()
Definition: gslib.cpp:492
struct gslib::findpts_data_3 * fdata3D
Definition: gslib.hpp:56
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1263
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:746
virtual void FreeData()
Definition: gslib.cpp:212
int GetBasisType() const
Definition: fe_coll.hpp:240
void Setup(Mesh &m, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
Definition: gslib.cpp:75
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.hpp:390
double default_interp_value
Definition: gslib.hpp:63
struct gslib::crystal * cr
Definition: gslib.hpp:57
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:466
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4843
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:534
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
int Dimension() const
Definition: mesh.hpp:911
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
void FindPoints(const Vector &point_pos, Array< unsigned int > &point_id)
Definition: gslib.cpp:919
Array< unsigned int > gsl_code
Definition: gslib.hpp:60
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
virtual void GetSimplexNodalCoordinates()
Definition: gslib.cpp:328
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:323
void FindPoints(const Vector &point_pos)
Definition: gslib.cpp:129
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Definition: mesh.cpp:1373
virtual void InterpolateGeneral(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:699
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:674
virtual void GetQuadHexNodalCoordinates()
Definition: gslib.cpp:290
struct gslib::findpts_data_2 * fdata2D
Definition: gslib.hpp:55
const Array< int > & GetDofMap() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe.hpp:2154
Class for integration point with weight.
Definition: intrules.hpp:25
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:4877
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1585
int dim
Definition: ex24.cpp:53
void Destroy()
Destroy a vector.
Definition: vector.hpp:590
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2388
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
struct gslib::comm * gsl_comm
Definition: gslib.hpp:58
void Interpolate(const Vector &point_pos, Array< unsigned int > &point_id, const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:986
IntegrationRule * ir_simplex
Definition: gslib.hpp:54
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:554
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
Vector coefficient defined by a vector GridFunction.
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
void Setup(Mesh &m, const int meshid, GridFunction *gfmax=NULL, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
Definition: gslib.cpp:847
virtual void InterpolateH1(const GridFunction &field_in, Vector &field_out)
Use GSLIB for communication and interpolation.
Definition: gslib.cpp:657
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:402
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285