MFEM  v4.2.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-2020, 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 #include "gslib.h"
23 
24 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
25 #pragma GCC diagnostic pop
26 #endif
27 
28 namespace mfem
29 {
30 
32  : mesh(NULL), meshsplit(NULL), ir_simplex(NULL),
33  fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
34  dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
35  avgtype(AvgType::ARITHMETIC)
36 {
37  gsl_comm = new comm;
38  cr = new crystal;
39 #ifdef MFEM_USE_MPI
40  int initialized;
41  MPI_Initialized(&initialized);
42  if (!initialized) { MPI_Init(NULL, NULL); }
43  MPI_Comm comm = MPI_COMM_WORLD;;
44  comm_init(gsl_comm, comm);
45 #else
46  comm_init(gsl_comm, 0);
47 #endif
48 }
49 
51 {
52  delete gsl_comm;
53  delete cr;
54  delete ir_simplex;
55  delete meshsplit;
56 }
57 
58 #ifdef MFEM_USE_MPI
60  : mesh(NULL), meshsplit(NULL), ir_simplex(NULL),
61  fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
62  dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
63  avgtype(AvgType::ARITHMETIC)
64 {
65  gsl_comm = new comm;
66  cr = new crystal;
67  comm_init(gsl_comm, _comm);
68 }
69 #endif
70 
71 void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
72  const int npt_max)
73 {
74  MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
75  MFEM_VERIFY(m.GetNumGeometries(m.Dimension()) == 1,
76  "Mixed meshes are not currently supported in FindPointsGSLIB.");
77 
78  // call FreeData if FindPointsGSLIB::Setup has been called already
79  if (setupflag) { FreeData(); }
80 
81  crystal_init(cr, gsl_comm);
82  mesh = &m;
83  dim = mesh->Dimension();
84  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
85  unsigned dof1D = fe->GetOrder() + 1;
86  const int gt = fe->GetGeomType();
87 
88  if (gt == Geometry::TRIANGLE || gt == Geometry::TETRAHEDRON ||
89  gt == Geometry::PRISM)
90  {
92  }
93  else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
94  {
96  }
97  else
98  {
99  MFEM_ABORT("Element type not currently supported in FindPointsGSLIB.");
100  }
101 
102  const int pts_cnt = gsl_mesh.Size()/dim,
103  NEtot = pts_cnt/(int)pow(dof1D, dim);
104 
105  if (dim == 2)
106  {
107  unsigned nr[2] = { dof1D, dof1D };
108  unsigned mr[2] = { 2*dof1D, 2*dof1D };
109  double * const elx[2] = { &gsl_mesh(0), &gsl_mesh(pts_cnt) };
110  fdata2D = findpts_setup_2(gsl_comm, elx, nr, NEtot, mr, bb_t,
111  pts_cnt, pts_cnt, npt_max, newt_tol);
112  }
113  else
114  {
115  unsigned nr[3] = { dof1D, dof1D, dof1D };
116  unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
117  double * const elx[3] =
118  { &gsl_mesh(0), &gsl_mesh(pts_cnt), &gsl_mesh(2*pts_cnt) };
119  fdata3D = findpts_setup_3(gsl_comm, elx, nr, NEtot, mr, bb_t,
120  pts_cnt, pts_cnt, npt_max, newt_tol);
121  }
122  setupflag = true;
123 }
124 
125 void FindPointsGSLIB::FindPoints(const Vector &point_pos)
126 {
127  MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
128  points_cnt = point_pos.Size() / dim;
134 
135  if (dim == 2)
136  {
137  const double *xv_base[2];
138  xv_base[0] = point_pos.GetData();
139  xv_base[1] = point_pos.GetData() + points_cnt;
140  unsigned xv_stride[2];
141  xv_stride[0] = sizeof(double);
142  xv_stride[1] = sizeof(double);
143  findpts_2(gsl_code.GetData(), sizeof(unsigned int),
144  gsl_proc.GetData(), sizeof(unsigned int),
145  gsl_elem.GetData(), sizeof(unsigned int),
146  gsl_ref.GetData(), sizeof(double) * dim,
147  gsl_dist.GetData(), sizeof(double),
148  xv_base, xv_stride, points_cnt, fdata2D);
149  }
150  else
151  {
152  const double *xv_base[3];
153  xv_base[0] = point_pos.GetData();
154  xv_base[1] = point_pos.GetData() + points_cnt;
155  xv_base[2] = point_pos.GetData() + 2*points_cnt;
156  unsigned xv_stride[3];
157  xv_stride[0] = sizeof(double);
158  xv_stride[1] = sizeof(double);
159  xv_stride[2] = sizeof(double);
160  findpts_3(gsl_code.GetData(), sizeof(unsigned int),
161  gsl_proc.GetData(), sizeof(unsigned int),
162  gsl_elem.GetData(), sizeof(unsigned int),
163  gsl_ref.GetData(), sizeof(double) * dim,
164  gsl_dist.GetData(), sizeof(double),
165  xv_base, xv_stride, points_cnt, fdata3D);
166  }
167 
168  // Set the element number and reference position to 0 for points not found
169  for (int i = 0; i < points_cnt; i++)
170  {
171  if (gsl_code[i] == 2)
172  {
173  gsl_elem[i] = 0;
174  for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
175  }
176  }
177 
178  // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for
179  // both simplices and quads.
181 }
182 
183 void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
184  const double bb_t, const double newt_tol,
185  const int npt_max)
186 {
187  if (!setupflag || (mesh != &m) )
188  {
189  Setup(m, bb_t, newt_tol, npt_max);
190  }
191  FindPoints(point_pos);
192 }
193 
194 void FindPointsGSLIB::Interpolate(const Vector &point_pos,
195  const GridFunction &field_in, Vector &field_out)
196 {
197  FindPoints(point_pos);
198  Interpolate(field_in, field_out);
199 }
200 
201 void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
202  const GridFunction &field_in, Vector &field_out)
203 {
204  FindPoints(m, point_pos);
205  Interpolate(field_in, field_out);
206 }
207 
209 {
210  if (!setupflag) { return; }
211  crystal_free(cr);
212  if (dim == 2)
213  {
214  findpts_free_2(fdata2D);
215  }
216  else
217  {
218  findpts_free_3(fdata3D);
219  }
223  gsl_mesh.Destroy();
224  gsl_ref.Destroy();
225  gsl_dist.Destroy();
226  setupflag = false;
227 }
228 
230  Vector &node_vals)
231 {
232  MFEM_ASSERT(gf_in.FESpace()->GetVDim() == 1, "Scalar function expected.");
233 
234  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
235  const Geometry::Type gt = fe->GetGeomType();
236  const int NE = mesh->GetNE();
237 
238  if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
239  {
240  const GridFunction *nodes = mesh->GetNodes();
241  const FiniteElementSpace *fes = nodes->FESpace();
242  const IntegrationRule &ir = fes->GetFE(0)->GetNodes();
243  const int dof_cnt = ir.GetNPoints();
244 
245  node_vals.SetSize(NE * dof_cnt);
246 
247  const TensorBasisElement *tbe =
248  dynamic_cast<const TensorBasisElement *>(fes->GetFE(0));
249  MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
250  const Array<int> &dof_map = tbe->GetDofMap();
251 
252  int pt_id = 0;
253  Vector vals_el;
254  for (int i = 0; i < NE; i++)
255  {
256  gf_in.GetValues(i, ir, vals_el);
257  for (int j = 0; j < dof_cnt; j++)
258  {
259  node_vals(pt_id++) = vals_el(dof_map[j]);
260  }
261  }
262  }
263  else if (gt == Geometry::TRIANGLE || gt == Geometry::TETRAHEDRON ||
264  gt == Geometry::PRISM)
265  {
266  const int dof_cnt = ir_simplex->GetNPoints();
267  node_vals.SetSize(NE * dof_cnt);
268 
269  int pt_id = 0;
270  Vector vals_el;
271  for (int j = 0; j < NE; j++)
272  {
273  gf_in.GetValues(j, *ir_simplex, vals_el);
274  for (int i = 0; i < dof_cnt; i++)
275  {
276  node_vals(pt_id++) = vals_el(i);
277  }
278  }
279  }
280  else
281  {
282  MFEM_ABORT("Element type not currently supported.");
283  }
284 }
285 
287 {
288  const GridFunction *nodes = mesh->GetNodes();
289  const FiniteElementSpace *fes = nodes->FESpace();
290 
291  const int NE = mesh->GetNE(),
292  dof_cnt = fes->GetFE(0)->GetDof(),
293  pts_cnt = NE * dof_cnt;
294  gsl_mesh.SetSize(dim * pts_cnt);
295 
296  const TensorBasisElement *tbe =
297  dynamic_cast<const TensorBasisElement *>(fes->GetFE(0));
298  MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
299  const Array<int> &dof_map = tbe->GetDofMap();
300 
301  DenseMatrix pos(dof_cnt, dim);
302  Vector posV(pos.Data(), dof_cnt * dim);
303  Array<int> xdofs(dof_cnt * dim);
304 
305  int pt_id = 0;
306  for (int i = 0; i < NE; i++)
307  {
308  fes->GetElementVDofs(i, xdofs);
309  nodes->GetSubVector(xdofs, posV);
310  for (int j = 0; j < dof_cnt; j++)
311  {
312  for (int d = 0; d < dim; d++)
313  {
314  gsl_mesh(pts_cnt * d + pt_id) = pos(dof_map[j], d);
315  }
316  pt_id++;
317  }
318  }
319 }
320 
322 {
323  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
324  const Geometry::Type gt = fe->GetGeomType();
325  const GridFunction *nodes = mesh->GetNodes();
326  const int NE = mesh->GetNE();
327  int NEsplit = 0;
328 
329  // Split the reference element into a reference submesh of quads or hexes.
330  if (gt == Geometry::TRIANGLE)
331  {
332  int Nvert = 7;
333  NEsplit = 3;
334  meshsplit = new Mesh(2, Nvert, NEsplit, 0, 2);
335 
336  const double quad_v[7][2] =
337  {
338  {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
339  {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
340  };
341  const int quad_e[3][4] =
342  {
343  {3, 4, 1, 0}, {4, 5, 2, 1}, {6, 5, 4, 3}
344  };
345 
346  for (int j = 0; j < Nvert; j++)
347  {
348  meshsplit->AddVertex(quad_v[j]);
349  }
350  for (int j = 0; j < NEsplit; j++)
351  {
352  int attribute = j + 1;
353  meshsplit->AddQuad(quad_e[j], attribute);
354  }
355  meshsplit->FinalizeQuadMesh(1, 1, true);
356  }
357  else if (gt == Geometry::TETRAHEDRON)
358  {
359  int Nvert = 15;
360  NEsplit = 4;
361  meshsplit = new Mesh(3, Nvert, NEsplit, 0, 3);
362 
363  const double hex_v[15][3] =
364  {
365  {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
366  {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
367  {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
368  {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
369  {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
370  };
371  const int hex_e[4][8] =
372  {
373  {0, 4, 10, 7, 6, 13, 14, 12},
374  {4, 1, 8, 10, 13, 5, 11, 14},
375  {13, 5, 11, 14, 6, 2, 9, 12},
376  {10, 8, 3, 7, 14, 11, 9, 12}
377  };
378 
379  for (int j = 0; j < Nvert; j++)
380  {
381  meshsplit->AddVertex(hex_v[j]);
382  }
383  for (int j = 0; j < NEsplit; j++)
384  {
385  int attribute = j + 1;
386  meshsplit->AddHex(hex_e[j], attribute);
387  }
388  meshsplit->FinalizeHexMesh(1, 1, true);
389  }
390  else if (gt == Geometry::PRISM)
391  {
392  int Nvert = 14;
393  NEsplit = 3;
394  meshsplit = new Mesh(3, Nvert, NEsplit, 0, 3);
395 
396  const double hex_v[14][3] =
397  {
398  {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
399  {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
400  {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
401  {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
402  };
403  const int hex_e[3][8] =
404  {
405  {3, 4, 1, 0, 10, 11, 8, 7},
406  {4, 5, 2, 1, 11, 12, 9, 8},
407  {6, 5, 4, 3, 13, 12, 11, 10}
408  };
409 
410  for (int j = 0; j < Nvert; j++)
411  {
412  meshsplit->AddVertex(hex_v[j]);
413  }
414  for (int j = 0; j < NEsplit; j++)
415  {
416  int attribute = j + 1;
417  meshsplit->AddHex(hex_e[j], attribute);
418  }
419  meshsplit->FinalizeHexMesh(1, 1, true);
420  }
421  else { MFEM_ABORT("Unsupported geometry type."); }
422 
423  // Curve the reference submesh.
424  H1_FECollection fec(fe->GetOrder(), dim);
425  FiniteElementSpace nodal_fes(meshsplit, &fec, dim);
426  meshsplit->SetNodalFESpace(&nodal_fes);
427 
428  const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
429  pts_cnt = NEsplit * dof_cnt;
430  Vector irlist(dim * pts_cnt);
431 
432  const TensorBasisElement *tbe =
433  dynamic_cast<const TensorBasisElement *>(nodal_fes.GetFE(0));
434  MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
435  const Array<int> &dof_map = tbe->GetDofMap();
436 
437  DenseMatrix pos(dof_cnt, dim);
438  Vector posV(pos.Data(), dof_cnt * dim);
439  Array<int> xdofs(dof_cnt * dim);
440 
441  // Create an IntegrationRule on the nodes of the reference submesh.
442  ir_simplex = new IntegrationRule(pts_cnt);
443  GridFunction *nodesplit = meshsplit->GetNodes();
444  int pt_id = 0;
445  for (int i = 0; i < NEsplit; i++)
446  {
447  nodal_fes.GetElementVDofs(i, xdofs);
448  nodesplit->GetSubVector(xdofs, posV);
449  for (int j = 0; j < dof_cnt; j++)
450  {
451  for (int d = 0; d < dim; d++)
452  {
453  irlist(pts_cnt * d + pt_id) = pos(dof_map[j], d);
454  }
455  ir_simplex->IntPoint(pt_id).x = irlist(pt_id);
456  ir_simplex->IntPoint(pt_id).y = irlist(pts_cnt + pt_id);
457  if (dim == 3)
458  {
459  ir_simplex->IntPoint(pt_id).z = irlist(2*pts_cnt + pt_id);
460  }
461  pt_id++;
462  }
463  }
464 
465  // Initialize gsl_mesh with the positions of the split physical elements.
466  pt_id = 0;
467  Vector locval(dim);
468  const int tot_pts_cnt = pts_cnt*NE;
469  gsl_mesh.SetSize(tot_pts_cnt*dim);
470  for (int j = 0; j < NE; j++)
471  {
472  for (int i = 0; i < dof_cnt*NEsplit; i++)
473  {
474  const IntegrationPoint &ip = ir_simplex->IntPoint(i);
475  nodes->GetVectorValue(j, ip, locval);
476  for (int d = 0; d < dim; d++)
477  {
478  gsl_mesh(tot_pts_cnt*d + pt_id) = locval(d);
479  }
480  pt_id++;
481  }
482  }
483 }
484 
486 {
489  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
490  const Geometry::Type gt = fe->GetGeomType();
491  int NEsplit = 0;
492 
493  gsl_mfem_ref -= -1.; // map [-1, 1] to
494  gsl_mfem_ref *= 0.5; // [0, 1]
495  if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { return; }
496 
497  H1_FECollection feclin(1, dim);
498  FiniteElementSpace nodal_fes_lin(meshsplit, &feclin, dim);
499  GridFunction gf_lin(&nodal_fes_lin);
500 
501  if (gt == Geometry::TRIANGLE)
502  {
503  const double quad_v[7][2] =
504  {
505  {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
506  {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
507  };
508  for (int k = 0; k < dim; k++)
509  {
510  for (int j = 0; j < gf_lin.Size()/dim; j++)
511  {
512  gf_lin(j+k*gf_lin.Size()/dim) = quad_v[j][k];
513  }
514  }
515  NEsplit = 3;
516  }
517  else if (gt == Geometry::TETRAHEDRON)
518  {
519  const double hex_v[15][3] =
520  {
521  {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
522  {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
523  {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
524  {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
525  {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
526  };
527  for (int k = 0; k < dim; k++)
528  {
529  for (int j = 0; j < gf_lin.Size()/dim; j++)
530  {
531  gf_lin(j+k*gf_lin.Size()/dim) = hex_v[j][k];
532  }
533  }
534  NEsplit = 4;
535  }
536  else if (gt == Geometry::PRISM)
537  {
538  const double hex_v[14][3] =
539  {
540  {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
541  {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
542  {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
543  {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
544  };
545  for (int k = 0; k < dim; k++)
546  {
547  for (int j = 0; j < gf_lin.Size()/dim; j++)
548  {
549  gf_lin(j+k*gf_lin.Size()/dim) = hex_v[j][k];
550  }
551  }
552  NEsplit = 3;
553  }
554  else
555  {
556  MFEM_ABORT("Element type not currently supported.");
557  }
558 
559  // Simplices are split into quads/hexes for GSLIB. For MFEM, we need to find
560  // the original element number and map the rst from micro to macro element.
561  for (int i = 0; i < points_cnt; i++)
562  {
563  if (gsl_code[i] == 2) { continue; }
564  int local_elem = gsl_elem[i]%NEsplit;
565  gsl_mfem_elem[i] = (gsl_elem[i] - local_elem)/NEsplit; // macro element number
566 
567  IntegrationPoint ip;
568  Vector mfem_ref(gsl_mfem_ref.GetData()+i*dim, dim);
569  ip.Set2(mfem_ref.GetData());
570  if (dim == 3) { ip.z = mfem_ref(2); }
571  gf_lin.GetVectorValue(local_elem, ip, mfem_ref); // map to rst of macro element
572  }
573 }
574 
576  Vector &field_out)
577 {
578  const int gf_order = field_in.FESpace()->GetFE(0)->GetOrder(),
579  mesh_order = mesh->GetNodalFESpace()->GetFE(0)->GetOrder();
580 
581  const FiniteElementCollection *fec_in = field_in.FESpace()->FEColl();
582  const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>(fec_in);
583  const L2_FECollection *fec_l2 = dynamic_cast<const L2_FECollection *>(fec_in);
584 
585  if (fec_h1 && gf_order == mesh_order &&
587  {
588  InterpolateH1(field_in, field_out);
589  return;
590  }
591  else
592  {
593  InterpolateGeneral(field_in, field_out);
594  if (!fec_l2 || avgtype == AvgType::NONE) { return; }
595  }
596 
597  // For points on element borders, project the L2 GridFunction to H1 and
598  // re-interpolate.
599  if (fec_l2)
600  {
601  Array<int> indl2;
602  for (int i = 0; i < points_cnt; i++)
603  {
604  if (gsl_code[i] == 1) { indl2.Append(i); }
605  }
606  if (indl2.Size() == 0) { return; } // no points on element borders
607 
608  Vector field_out_l2(field_out.Size());
609  VectorGridFunctionCoefficient field_in_dg(&field_in);
610  int gf_order_h1 = std::max(gf_order, 1); // H1 should be at least order 1
611  H1_FECollection fec(gf_order_h1, dim);
612  const int ncomp = field_in.FESpace()->GetVDim();
613  FiniteElementSpace fes(mesh, &fec, ncomp);
614  GridFunction field_in_h1(&fes);
615 
616  if (avgtype == AvgType::ARITHMETIC)
617  {
618  field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
619  }
620  else if (avgtype == AvgType::HARMONIC)
621  {
622  field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
623  }
624  else
625  {
626  MFEM_ABORT("Invalid averaging type.");
627  }
628 
629  if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
630  {
631  InterpolateH1(field_in_h1, field_out_l2);
632  }
633  else
634  {
635  InterpolateGeneral(field_in_h1, field_out_l2);
636  }
637 
638  // Copy interpolated values for the points on element border
639  for (int j = 0; j < ncomp; j++)
640  {
641  for (int i = 0; i < indl2.Size(); i++)
642  {
643  int idx = indl2[i] + j*points_cnt;
644  field_out(idx) = field_out_l2(idx);
645  }
646  }
647  }
648 }
649 
651  Vector &field_out)
652 {
653  FiniteElementSpace ind_fes(mesh, field_in.FESpace()->FEColl());
654  GridFunction field_in_scalar(&ind_fes);
655  Vector node_vals;
656 
657  const int ncomp = field_in.FESpace()->GetVDim(),
658  points_fld = field_in.Size() / ncomp,
660 
661  field_out.SetSize(points_cnt*ncomp);
662  field_out = default_interp_value;
663 
664  for (int i = 0; i < ncomp; i++)
665  {
666  const int dataptrin = i*points_fld,
667  dataptrout = i*points_cnt;
668  field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin, points_fld);
669  GetNodeValues(field_in_scalar, node_vals);
670 
671  if (dim==2)
672  {
673  findpts_eval_2(field_out.GetData()+dataptrout, sizeof(double),
674  gsl_code.GetData(), sizeof(unsigned int),
675  gsl_proc.GetData(), sizeof(unsigned int),
676  gsl_elem.GetData(), sizeof(unsigned int),
677  gsl_ref.GetData(), sizeof(double) * dim,
678  points_cnt, node_vals.GetData(), fdata2D);
679  }
680  else
681  {
682  findpts_eval_3(field_out.GetData()+dataptrout, sizeof(double),
683  gsl_code.GetData(), sizeof(unsigned int),
684  gsl_proc.GetData(), sizeof(unsigned int),
685  gsl_elem.GetData(), sizeof(unsigned int),
686  gsl_ref.GetData(), sizeof(double) * dim,
687  points_cnt, node_vals.GetData(), fdata3D);
688  }
689  }
690 }
691 
693  Vector &field_out)
694 {
695  int ncomp = field_in.VectorDim(),
696  nptorig = points_cnt,
697  npt = points_cnt;
698 
699  field_out.SetSize(points_cnt*ncomp);
700  field_out = default_interp_value;
701 
702  if (gsl_comm->np == 1) // serial
703  {
704  for (int index = 0; index < npt; index++)
705  {
706  if (gsl_code[index] == 2) { continue; }
707  IntegrationPoint ip;
709  if (dim == 3) { ip.z = gsl_mfem_ref(index*dim + 2); }
710  Vector localval(ncomp);
711  field_in.GetVectorValue(gsl_mfem_elem[index], ip, localval);
712  for (int i = 0; i < ncomp; i++)
713  {
714  field_out(index + i*npt) = localval(i);
715  }
716  }
717  }
718  else // parallel
719  {
720  // Determine number of points to be sent
721  int nptsend = 0;
722  for (int index = 0; index < npt; index++)
723  {
724  if (gsl_code[index] != 2) { nptsend +=1; }
725  }
726 
727  // Pack data to send via crystal router
728  struct array *outpt = new array;
729  struct out_pt { double r[3], ival; uint index, el, proc; };
730  struct out_pt *pt;
731  array_init(struct out_pt, outpt, nptsend);
732  outpt->n=nptsend;
733  pt = (struct out_pt *)outpt->ptr;
734  for (int index = 0; index < npt; index++)
735  {
736  if (gsl_code[index] == 2) { continue; }
737  for (int d = 0; d < dim; ++d) { pt->r[d]= gsl_mfem_ref(index*dim + d); }
738  pt->index = index;
739  pt->proc = gsl_proc[index];
740  pt->el = gsl_mfem_elem[index];
741  ++pt;
742  }
743 
744  // Transfer data to target MPI ranks
745  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
746 
747  if (ncomp == 1)
748  {
749  // Interpolate the grid function
750  npt = outpt->n;
751  pt = (struct out_pt *)outpt->ptr;
752  for (int index = 0; index < npt; index++)
753  {
754  IntegrationPoint ip;
755  ip.Set3(&pt->r[0]);
756  pt->ival = field_in.GetValue(pt->el, ip, 1);
757  ++pt;
758  }
759 
760  // Transfer data back to source MPI rank
761  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
762  npt = outpt->n;
763  pt = (struct out_pt *)outpt->ptr;
764  for (int index = 0; index < npt; index++)
765  {
766  field_out(pt->index) = pt->ival;
767  ++pt;
768  }
769  array_free(outpt);
770  delete outpt;
771  }
772  else // ncomp > 1
773  {
774  // Interpolate data and store in a Vector
775  npt = outpt->n;
776  pt = (struct out_pt *)outpt->ptr;
777  Vector vec_int_vals(npt*ncomp);
778  for (int index = 0; index < npt; index++)
779  {
780  IntegrationPoint ip;
781  ip.Set3(&pt->r[0]);
782  Vector localval(vec_int_vals.GetData()+index*ncomp, ncomp);
783  field_in.GetVectorValue(pt->el, ip, localval);
784  ++pt;
785  }
786 
787  // Save index and proc data in a struct
788  struct array *savpt = new array;
789  struct sav_pt { uint index, proc; };
790  struct sav_pt *spt;
791  array_init(struct sav_pt, savpt, npt);
792  savpt->n=npt;
793  spt = (struct sav_pt *)savpt->ptr;
794  pt = (struct out_pt *)outpt->ptr;
795  for (int index = 0; index < npt; index++)
796  {
797  spt->index = pt->index;
798  spt->proc = pt->proc;
799  ++pt; ++spt;
800  }
801 
802  array_free(outpt);
803  delete outpt;
804 
805  // Copy data from save struct to send struct and send component wise
806  struct array *sendpt = new array;
807  struct send_pt { double ival; uint index, proc; };
808  struct send_pt *sdpt;
809  for (int j = 0; j < ncomp; j++)
810  {
811  array_init(struct send_pt, sendpt, npt);
812  sendpt->n=npt;
813  spt = (struct sav_pt *)savpt->ptr;
814  sdpt = (struct send_pt *)sendpt->ptr;
815  for (int index = 0; index < npt; index++)
816  {
817  sdpt->index = spt->index;
818  sdpt->proc = spt->proc;
819  sdpt->ival = vec_int_vals(j + index*ncomp);
820  ++sdpt; ++spt;
821  }
822 
823  sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
824  sdpt = (struct send_pt *)sendpt->ptr;
825  for (int index = 0; index < nptorig; index++)
826  {
827  int idx = sdpt->index + j*nptorig;
828  field_out(idx) = sdpt->ival;
829  ++sdpt;
830  }
831  array_free(sendpt);
832  }
833  array_free(savpt);
834  delete sendpt;
835  delete savpt;
836  } // ncomp > 1
837  } // parallel
838 }
839 
840 } // namespace mfem
841 
842 #endif // MFEM_USE_GSLIB
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for all finite elements.
Definition: fe.hpp:235
void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:575
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Array< unsigned int > gsl_elem
Definition: gslib.hpp:58
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:1293
struct findpts_data_3 * fdata3D
Definition: gslib.hpp:54
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:414
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:4697
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
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:173
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2387
Array< unsigned int > gsl_mfem_elem
Definition: gslib.hpp:58
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:495
T * GetData()
Returns the data.
Definition: array.hpp:98
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:319
int VectorDim() const
Definition: gridfunc.cpp:300
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
Array< unsigned int > gsl_proc
Definition: gslib.hpp:58
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
Definition: mesh.cpp:2432
void GetNodeValues(const GridFunction &gf_in, Vector &node_vals)
Get GridFunction from MFEM format to GSLIB format.
Definition: gslib.cpp:229
void DeleteAll()
Delete the whole array.
Definition: array.hpp:821
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:312
void MapRefPosAndElemIndices()
Definition: gslib.cpp:485
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1232
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:726
int GetBasisType() const
Definition: fe_coll.hpp:184
struct comm * gsl_comm
Definition: gslib.hpp:56
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:71
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.hpp:382
double default_interp_value
Definition: gslib.hpp:61
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:455
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4126
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:389
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:595
int Dimension() const
Definition: mesh.hpp:788
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
Array< unsigned int > gsl_code
Definition: gslib.hpp:58
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:96
void GetSimplexNodalCoordinates()
Definition: gslib.cpp:321
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:315
void FindPoints(const Vector &point_pos)
Definition: gslib.cpp:125
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:1342
void InterpolateGeneral(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:692
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:654
void GetQuadHexNodalCoordinates()
Definition: gslib.cpp:286
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:2029
Class for integration point with weight.
Definition: intrules.hpp:25
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:4160
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition: mesh.cpp:1547
int dim
Definition: ex24.cpp:53
struct crystal * cr
Definition: gslib.hpp:55
void Destroy()
Destroy a vector.
Definition: vector.hpp:530
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:1798
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
IntegrationRule * ir_simplex
Definition: gslib.hpp:52
Vector data type.
Definition: vector.hpp:51
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:414
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
Vector coefficient defined by a vector GridFunction.
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
struct findpts_data_2 * fdata2D
Definition: gslib.hpp:53
void InterpolateH1(const GridFunction &field_in, Vector &field_out)
Use GSLIB for communication and interpolation.
Definition: gslib.cpp:650
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:391
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:221