MFEM  v4.6.0
Finite element discretization library
gslib.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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),
37  fec_map_lin(NULL),
38  fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
39  dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
40  avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
41 {
42  mesh_split.SetSize(4);
43  ir_split.SetSize(4);
45  gf_rst_map.SetSize(4);
46  for (int i = 0; i < mesh_split.Size(); i++)
47  {
48  mesh_split[i] = NULL;
49  ir_split[i] = NULL;
50  fes_rst_map[i] = NULL;
51  gf_rst_map[i] = NULL;
52  }
53 
54  gsl_comm = new gslib::comm;
55  cr = new gslib::crystal;
56 #ifdef MFEM_USE_MPI
57  int initialized;
58  MPI_Initialized(&initialized);
59  if (!initialized) { MPI_Init(NULL, NULL); }
60  MPI_Comm comm = MPI_COMM_WORLD;
61  comm_init(gsl_comm, comm);
62 #else
63  comm_init(gsl_comm, 0);
64 #endif
65 }
66 
68 {
69  delete gsl_comm;
70  delete cr;
71  for (int i = 0; i < 4; i++)
72  {
73  if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = NULL; }
74  if (ir_split[i]) { delete ir_split[i]; ir_split[i] = NULL; }
75  if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = NULL; }
76  if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = NULL; }
77  }
78  if (fec_map_lin) { delete fec_map_lin; fec_map_lin = NULL; }
79 }
80 
81 #ifdef MFEM_USE_MPI
83  : mesh(NULL),
84  fec_map_lin(NULL),
85  fdata2D(NULL), fdata3D(NULL), cr(NULL), gsl_comm(NULL),
86  dim(-1), points_cnt(0), setupflag(false), default_interp_value(0),
87  avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
88 {
89  mesh_split.SetSize(4);
90  ir_split.SetSize(4);
92  gf_rst_map.SetSize(4);
93  for (int i = 0; i < mesh_split.Size(); i++)
94  {
95  mesh_split[i] = NULL;
96  ir_split[i] = NULL;
97  fes_rst_map[i] = NULL;
98  gf_rst_map[i] = NULL;
99  }
100 
101  gsl_comm = new gslib::comm;
102  cr = new gslib::crystal;
103  comm_init(gsl_comm, comm_);
104 }
105 #endif
106 
107 void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
108  const int npt_max)
109 {
110  MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
111  const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
112 
113  // call FreeData if FindPointsGSLIB::Setup has been called already
114  if (setupflag) { FreeData(); }
115 
116  crystal_init(cr, gsl_comm);
117  mesh = &m;
118  dim = mesh->Dimension();
119  unsigned dof1D = meshOrder + 1;
120 
122  if (dim == 2)
123  {
124  if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
125  ir_split[0] = new IntegrationRule(3*pow(dof1D, dim));
127 
128  if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
129  ir_split[1] = new IntegrationRule(pow(dof1D, dim));
131  }
132  else if (dim == 3)
133  {
134  if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
135  ir_split[0] = new IntegrationRule(pow(dof1D, dim));
137 
138  if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
139  ir_split[1] = new IntegrationRule(4*pow(dof1D, dim));
141 
142  if (ir_split[2]) { delete ir_split[2]; ir_split[2] = NULL; }
143  ir_split[2] = new IntegrationRule(3*pow(dof1D, dim));
145 
146  if (ir_split[3]) { delete ir_split[3]; ir_split[3] = NULL; }
147  ir_split[3] = new IntegrationRule(8*pow(dof1D, dim));
149  }
150 
152 
153  const int pts_cnt = gsl_mesh.Size()/dim,
154  NEtot = pts_cnt/(int)pow(dof1D, dim);
155 
156  if (dim == 2)
157  {
158  unsigned nr[2] = { dof1D, dof1D };
159  unsigned mr[2] = { 2*dof1D, 2*dof1D };
160  double * const elx[2] = { &gsl_mesh(0), &gsl_mesh(pts_cnt) };
161  fdata2D = findpts_setup_2(gsl_comm, elx, nr, NEtot, mr, bb_t,
162  pts_cnt, pts_cnt, npt_max, newt_tol);
163  }
164  else
165  {
166  unsigned nr[3] = { dof1D, dof1D, dof1D };
167  unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
168  double * const elx[3] =
169  { &gsl_mesh(0), &gsl_mesh(pts_cnt), &gsl_mesh(2*pts_cnt) };
170  fdata3D = findpts_setup_3(gsl_comm, elx, nr, NEtot, mr, bb_t,
171  pts_cnt, pts_cnt, npt_max, newt_tol);
172  }
173  setupflag = true;
174 }
175 
176 void FindPointsGSLIB::FindPoints(const Vector &point_pos,
177  int point_pos_ordering)
178 {
179  MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
180  points_cnt = point_pos.Size() / dim;
186 
187  auto xvFill = [&](const double *xv_base[], unsigned xv_stride[])
188  {
189  for (int d = 0; d < dim; d++)
190  {
191  if (point_pos_ordering == Ordering::byNODES)
192  {
193  xv_base[d] = point_pos.GetData() + d*points_cnt;
194  xv_stride[d] = sizeof(double);
195  }
196  else
197  {
198  xv_base[d] = point_pos.GetData() + d;
199  xv_stride[d] = dim*sizeof(double);
200  }
201  }
202  };
203  if (dim == 2)
204  {
205  const double *xv_base[2];
206  unsigned xv_stride[2];
207  xvFill(xv_base, xv_stride);
208  findpts_2(gsl_code.GetData(), sizeof(unsigned int),
209  gsl_proc.GetData(), sizeof(unsigned int),
210  gsl_elem.GetData(), sizeof(unsigned int),
211  gsl_ref.GetData(), sizeof(double) * dim,
212  gsl_dist.GetData(), sizeof(double),
213  xv_base, xv_stride, points_cnt, fdata2D);
214  }
215  else // dim == 3
216  {
217  const double *xv_base[3];
218  unsigned xv_stride[3];
219  xvFill(xv_base, xv_stride);
220  findpts_3(gsl_code.GetData(), sizeof(unsigned int),
221  gsl_proc.GetData(), sizeof(unsigned int),
222  gsl_elem.GetData(), sizeof(unsigned int),
223  gsl_ref.GetData(), sizeof(double) * dim,
224  gsl_dist.GetData(), sizeof(double),
225  xv_base, xv_stride, points_cnt, fdata3D);
226  }
227 
228  // Set the element number and reference position to 0 for points not found
229  for (int i = 0; i < points_cnt; i++)
230  {
231  if (gsl_code[i] == 2 ||
232  (gsl_code[i] == 1 && gsl_dist(i) > bdr_tol))
233  {
234  gsl_elem[i] = 0;
235  for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
236  gsl_code[i] = 2;
237  }
238  }
239 
240  // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for
241  // both simplices and quads.
243 }
244 
245 void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
246  int point_pos_ordering, const double bb_t,
247  const double newt_tol, const int npt_max)
248 {
249  if (!setupflag || (mesh != &m) )
250  {
251  Setup(m, bb_t, newt_tol, npt_max);
252  }
253  FindPoints(point_pos, point_pos_ordering);
254 }
255 
256 void FindPointsGSLIB::Interpolate(const Vector &point_pos,
257  const GridFunction &field_in, Vector &field_out,
258  int point_pos_ordering)
259 {
260  FindPoints(point_pos, point_pos_ordering);
261  Interpolate(field_in, field_out);
262 }
263 
264 void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
265  const GridFunction &field_in, Vector &field_out,
266  int point_pos_ordering)
267 {
268  FindPoints(m, point_pos, point_pos_ordering);
269  Interpolate(field_in, field_out);
270 }
271 
273 {
274  if (!setupflag) { return; }
275  crystal_free(cr);
276  if (dim == 2)
277  {
278  findpts_free_2(fdata2D);
279  }
280  else
281  {
282  findpts_free_3(fdata3D);
283  }
287  gsl_mesh.Destroy();
288  gsl_ref.Destroy();
289  gsl_dist.Destroy();
290  for (int i = 0; i < 4; i++)
291  {
292  if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = NULL; }
293  if (ir_split[i]) { delete ir_split[i]; ir_split[i] = NULL; }
294  if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = NULL; }
295  if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = NULL; }
296  }
297  if (fec_map_lin) { delete fec_map_lin; fec_map_lin = NULL; }
298  setupflag = false;
299 }
300 
302 {
303  fec_map_lin = new H1_FECollection(1, dim);
304  if (mesh->Dimension() == 2)
305  {
306  int Nvert = 7;
307  int NEsplit = 3;
308  mesh_split[0] = new Mesh(2, Nvert, NEsplit, 0, 2);
309 
310  const double quad_v[7][2] =
311  {
312  {0, 0}, {0.5, 0}, {1, 0}, {0, 0.5},
313  {1./3., 1./3.}, {0.5, 0.5}, {0, 1}
314  };
315  const int quad_e[3][4] =
316  {
317  {0, 1, 4, 3}, {1, 2, 5, 4}, {3, 4, 5, 6}
318  };
319 
320  for (int j = 0; j < Nvert; j++)
321  {
322  mesh_split[0]->AddVertex(quad_v[j]);
323  }
324  for (int j = 0; j < NEsplit; j++)
325  {
326  int attribute = j + 1;
327  mesh_split[0]->AddQuad(quad_e[j], attribute);
328  }
329  mesh_split[0]->FinalizeQuadMesh(1, 1, true);
330 
332  gf_rst_map[0] = new GridFunction(fes_rst_map[0]);
333  const int npt = gf_rst_map[0]->Size()/dim;
334  for (int k = 0; k < dim; k++)
335  {
336  for (int j = 0; j < npt; j++)
337  {
338  (*gf_rst_map[0])(j+k*npt) = quad_v[j][k];
339  }
340  }
341 
342  mesh_split[1] = new Mesh(Mesh::MakeCartesian2D(1, 1,
344  }
345  else if (mesh->Dimension() == 3)
346  {
347  mesh_split[0] = new Mesh(Mesh::MakeCartesian3D(1, 1, 1,
349  // Tetrahedron
350  {
351  int Nvert = 15;
352  int NEsplit = 4;
353  mesh_split[1] = new Mesh(3, Nvert, NEsplit, 0, 3);
354 
355  const double hex_v[15][3] =
356  {
357  {0, 0, 0.}, {1, 0., 0.}, {0., 1., 0.}, {0, 0., 1.},
358  {0.5, 0., 0.}, {0.5, 0.5, 0.}, {0., 0.5, 0.},
359  {0., 0., 0.5}, {0.5, 0., 0.5}, {0., 0.5, 0.5},
360  {1./3., 0., 1./3.}, {1./3., 1./3., 1./3.}, {0, 1./3., 1./3.},
361  {1./3., 1./3., 0}, {0.25, 0.25, 0.25}
362  };
363  const int hex_e[4][8] =
364  {
365  {7, 10, 4, 0, 12, 14, 13, 6},
366  {10, 8, 1, 4, 14, 11, 5, 13},
367  {14, 11, 5, 13, 12, 9, 2, 6},
368  {7, 3, 8, 10, 12, 9, 11, 14}
369  };
370 
371  for (int j = 0; j < Nvert; j++)
372  {
373  mesh_split[1]->AddVertex(hex_v[j]);
374  }
375  for (int j = 0; j < NEsplit; j++)
376  {
377  int attribute = j + 1;
378  mesh_split[1]->AddHex(hex_e[j], attribute);
379  }
380  mesh_split[1]->FinalizeHexMesh(1, 1, true);
381 
383  gf_rst_map[1] = new GridFunction(fes_rst_map[1]);
384  const int npt = gf_rst_map[1]->Size()/dim;
385  for (int k = 0; k < dim; k++)
386  {
387  for (int j = 0; j < npt; j++)
388  {
389  (*gf_rst_map[1])(j+k*npt) = hex_v[j][k];
390  }
391  }
392  }
393  // Prism
394  {
395  int Nvert = 14;
396  int NEsplit = 3;
397  mesh_split[2] = new Mesh(3, Nvert, NEsplit, 0, 3);
398 
399  const double hex_v[14][3] =
400  {
401  {0, 0, 0}, {0.5, 0, 0}, {1, 0, 0}, {0, 0.5, 0},
402  {1./3., 1./3., 0}, {0.5, 0.5, 0}, {0, 1, 0},
403  {0, 0, 1}, {0.5, 0, 1}, {1, 0, 1}, {0, 0.5, 1},
404  {1./3., 1./3., 1}, {0.5, 0.5, 1}, {0, 1, 1}
405  };
406  const int hex_e[3][8] =
407  {
408  {0, 1, 4, 3, 7, 8, 11, 10},
409  {1, 2, 5, 4, 8, 9, 12, 11},
410  {3, 4, 5, 6, 10, 11, 12, 13}
411  };
412 
413  for (int j = 0; j < Nvert; j++)
414  {
415  mesh_split[2]->AddVertex(hex_v[j]);
416  }
417  for (int j = 0; j < NEsplit; j++)
418  {
419  int attribute = j + 1;
420  mesh_split[2]->AddHex(hex_e[j], attribute);
421  }
422  mesh_split[2]->FinalizeHexMesh(1, 1, true);
423 
425  gf_rst_map[2] = new GridFunction(fes_rst_map[2]);
426  const int npt = gf_rst_map[2]->Size()/dim;
427  for (int k = 0; k < dim; k++)
428  {
429  for (int j = 0; j < npt; j++)
430  {
431  (*gf_rst_map[2])(j+k*npt) = hex_v[j][k];
432  }
433  }
434  }
435  // Pyramid
436  {
437  int Nvert = 23;
438  int NEsplit = 8;
439  mesh_split[3] = new Mesh(3, Nvert, NEsplit, 0, 3);
440 
441  const double hex_v[23][3] =
442  {
443  {0.0000, 0.0000, 0.0000}, {0.5000, 0.0000, 0.0000},
444  {0.0000, 0.0000, 0.5000}, {0.3333, 0.0000, 0.3333},
445  {0.0000, 0.5000, 0.0000}, {0.3333, 0.3333, 0.0000},
446  {0.0000, 0.3333, 0.3333}, {0.2500, 0.2500, 0.2500},
447  {1.0000, 0.0000, 0.0000}, {0.5000, 0.0000, 0.5000},
448  {0.5000, 0.5000, 0.0000}, {0.3333, 0.3333, 0.3333},
449  {0.0000, 1.0000, 0.0000}, {0.0000, 0.5000, 0.5000},
450  {0.0000, 0.0000, 1.0000}, {1.0000, 0.5000, 0.0000},
451  {0.6667, 0.3333, 0.3333}, {0.6667, 0.6667, 0.0000},
452  {0.5000, 0.5000, 0.2500}, {1.0000, 1.0000, 0.0000},
453  {0.5000, 0.5000, 0.5000}, {0.5000, 1.0000, 0.0000},
454  {0.3333, 0.6667, 0.3333}
455  };
456  const int hex_e[8][8] =
457  {
458  {2, 3, 1, 0, 6, 7, 5, 4}, {3, 9, 8, 1, 7, 11, 10, 5},
459  {7, 11, 10, 5, 6, 13, 12, 4}, {2, 14, 9, 3, 6, 13, 11, 7},
460  {9, 16, 15, 8, 11, 18, 17, 10}, {16, 20, 19, 15, 18, 22, 21, 17},
461  {18, 22, 21, 17, 11, 13, 12, 10}, {9, 14, 20, 16, 11, 13, 22, 18}
462  };
463 
464  for (int j = 0; j < Nvert; j++)
465  {
466  mesh_split[3]->AddVertex(hex_v[j]);
467  }
468  for (int j = 0; j < NEsplit; j++)
469  {
470  int attribute = j + 1;
471  mesh_split[3]->AddHex(hex_e[j], attribute);
472  }
473  mesh_split[3]->FinalizeHexMesh(1, 1, true);
474 
476  gf_rst_map[3] = new GridFunction(fes_rst_map[3]);
477  const int npt = gf_rst_map[3]->Size()/dim;
478  for (int k = 0; k < dim; k++)
479  {
480  for (int j = 0; j < npt; j++)
481  {
482  (*gf_rst_map[3])(j+k*npt) = hex_v[j][k];
483  }
484  }
485  }
486  }
487 
488  NE_split_total = 0;
491  int NEsplit = 0;
492  for (int e = 0; e < mesh->GetNE(); e++)
493  {
494  const Geometry::Type gt = mesh->GetElement(e)->GetGeometryType();
495  if (gt == Geometry::TRIANGLE || gt == Geometry::PRISM)
496  {
497  NEsplit = 3;
498  }
499  else if (gt == Geometry::TETRAHEDRON)
500  {
501  NEsplit = 4;
502  }
503  else if (gt == Geometry::PYRAMID)
504  {
505  NEsplit = 8;
506  }
507  else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
508  {
509  NEsplit = 1;
510  }
511  else
512  {
513  MFEM_ABORT("Unsupported geometry type.");
514  }
515  NE_split_total += NEsplit;
516  for (int i = 0; i < NEsplit; i++)
517  {
520  }
521  }
522 }
523 
525  IntegrationRule *irule,
526  int order)
527 {
528  H1_FECollection fec(order, dim);
529  FiniteElementSpace nodal_fes(meshin, &fec, dim);
530  meshin->SetNodalFESpace(&nodal_fes);
531  const int NEsplit = meshin->GetNE();
532 
533  const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
534  pts_cnt = NEsplit * dof_cnt;
535  Vector irlist(dim * pts_cnt);
536 
537  const TensorBasisElement *tbe =
538  dynamic_cast<const TensorBasisElement *>(nodal_fes.GetFE(0));
539  MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
540  const Array<int> &dof_map = tbe->GetDofMap();
541 
542  DenseMatrix pos(dof_cnt, dim);
543  Vector posV(pos.Data(), dof_cnt * dim);
544  Array<int> xdofs(dof_cnt * dim);
545 
546  // Create an IntegrationRule on the nodes of the reference submesh.
547  MFEM_ASSERT(irule->GetNPoints() == pts_cnt, "IntegrationRule does not have"
548  "the correct number of points.");
549  GridFunction *nodesplit = meshin->GetNodes();
550  int pt_id = 0;
551  for (int i = 0; i < NEsplit; i++)
552  {
553  nodal_fes.GetElementVDofs(i, xdofs);
554  nodesplit->GetSubVector(xdofs, posV);
555  for (int j = 0; j < dof_cnt; j++)
556  {
557  for (int d = 0; d < dim; d++)
558  {
559  irlist(pts_cnt * d + pt_id) = pos(dof_map[j], d);
560  }
561  irule->IntPoint(pt_id).x = irlist(pt_id);
562  irule->IntPoint(pt_id).y = irlist(pts_cnt + pt_id);
563  if (dim == 3)
564  {
565  irule->IntPoint(pt_id).z = irlist(2*pts_cnt + pt_id);
566  }
567  pt_id++;
568  }
569  }
570 }
571 
573  Vector &node_vals)
574 {
575  const GridFunction *nodes = gf_in;
576  const FiniteElementSpace *fes = nodes->FESpace();
577  const int NE = mesh->GetNE();
578  const int vdim = fes->GetVDim();
579 
580  IntegrationRule *ir_split_temp = NULL;
581 
582  const int maxOrder = fes->GetMaxElementOrder();
583  const int dof_1D = maxOrder+1;
584  const int pts_el = std::pow(dof_1D, dim);
585  const int pts_cnt = NE_split_total * pts_el;
586  node_vals.SetSize(vdim * pts_cnt);
587  node_vals = 0.0;
588 
589  int gsl_mesh_pt_index = 0;
590 
591  for (int e = 0; e < NE; e++)
592  {
593  const FiniteElement *fe = fes->GetFE(e);
594  const Geometry::Type gt = fe->GetGeomType();
595  bool el_to_split = true;
596  if (gt == Geometry::TRIANGLE)
597  {
598  ir_split_temp = ir_split[0];
599  }
600  else if (gt == Geometry::TETRAHEDRON)
601  {
602  ir_split_temp = ir_split[1];
603  }
604  else if (gt == Geometry::PRISM)
605  {
606  ir_split_temp = ir_split[2];
607  }
608  else if (gt == Geometry::PYRAMID)
609  {
610  ir_split_temp = ir_split[3];
611  }
612  else if (gt == Geometry::SQUARE)
613  {
614  ir_split_temp = ir_split[1];
615  el_to_split = gf_in->FESpace()->IsVariableOrder();
616  }
617  else if (gt == Geometry::CUBE)
618  {
619  ir_split_temp = ir_split[0];
620  el_to_split = gf_in->FESpace()->IsVariableOrder();
621  }
622  else
623  {
624  MFEM_ABORT("Unsupported geometry type.");
625  }
626 
627  if (el_to_split) // Triangle/Tet/Prism or Quads/Hex but variable order
628  {
629  // Fill gsl_mesh with location of split points.
630  Vector locval(vdim);
631  for (int i = 0; i < ir_split_temp->GetNPoints(); i++)
632  {
633  const IntegrationPoint &ip = ir_split_temp->IntPoint(i);
634  nodes->GetVectorValue(e, ip, locval);
635  for (int d = 0; d < vdim; d++)
636  {
637  node_vals(pts_cnt*d + gsl_mesh_pt_index) = locval(d);
638  }
639  gsl_mesh_pt_index++;
640  }
641  }
642  else // Quad/Hex and constant polynomial order
643  {
644  const int dof_cnt_split = fe->GetDof();
645 
646  const TensorBasisElement *tbe =
647  dynamic_cast<const TensorBasisElement *>(fes->GetFE(e));
648  MFEM_VERIFY(tbe != NULL, "TensorBasis FiniteElement expected.");
649  Array<int> dof_map(dof_cnt_split);
650  const Array<int> &dm = tbe->GetDofMap();
651  if (dm.Size() > 0) { dof_map = dm; }
652  else { for (int i = 0; i < dof_cnt_split; i++) { dof_map[i] = i; } }
653 
654  DenseMatrix pos(dof_cnt_split, vdim);
655  Vector posV(pos.Data(), dof_cnt_split * vdim);
656  Array<int> xdofs(dof_cnt_split * vdim);
657 
658  fes->GetElementVDofs(e, xdofs);
659  nodes->GetSubVector(xdofs, posV);
660  for (int j = 0; j < dof_cnt_split; j++)
661  {
662  for (int d = 0; d < vdim; d++)
663  {
664  node_vals(pts_cnt * d + gsl_mesh_pt_index) = pos(dof_map[j], d);
665  }
666  gsl_mesh_pt_index++;
667  }
668  }
669  }
670 }
671 
672 
674 {
677 
678  gsl_mfem_ref -= -1.; // map [-1, 1] to
679  gsl_mfem_ref *= 0.5; // [0, 1]
680 
681  int nptorig = points_cnt,
682  npt = points_cnt;
683 
684  GridFunction *gf_rst_map_temp = NULL;
685  int nptsend = 0;
686 
687  for (int index = 0; index < npt; index++)
688  {
689  if (gsl_code[index] != 2 && gsl_proc[index] != gsl_comm->id)
690  {
691  nptsend +=1;
692  }
693  }
694 
695  // Pack data to send via crystal router
696  struct gslib::array *outpt = new gslib::array;
697  struct out_pt { double r[3]; uint index, el, proc; };
698  struct out_pt *pt;
699  array_init(struct out_pt, outpt, nptsend);
700  outpt->n=nptsend;
701  pt = (struct out_pt *)outpt->ptr;
702  for (int index = 0; index < npt; index++)
703  {
704  if (gsl_code[index] == 2 || gsl_proc[index] == gsl_comm->id)
705  {
706  continue;
707  }
708  for (int d = 0; d < dim; ++d)
709  {
710  pt->r[d]= gsl_mfem_ref(index*dim + d);
711  }
712  pt->index = index;
713  pt->proc = gsl_proc[index];
714  pt->el = gsl_elem[index];
715  ++pt;
716  }
717 
718  // Transfer data to target MPI ranks
719  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
720 
721  // Map received points
722  npt = outpt->n;
723  pt = (struct out_pt *)outpt->ptr;
724  for (int index = 0; index < npt; index++)
725  {
726  IntegrationPoint ip;
727  ip.Set3(&pt->r[0]);
728  const int elem = pt->el;
729  const int mesh_elem = split_element_map[elem];
730  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(mesh_elem);
731  const Geometry::Type gt = fe->GetGeomType();
732  pt->el = mesh_elem;
733 
734  if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { ++pt; continue; }
735  else if (gt == Geometry::TRIANGLE)
736  {
737  gf_rst_map_temp = gf_rst_map[0];
738  }
739  else if (gt == Geometry::TETRAHEDRON)
740  {
741  gf_rst_map_temp = gf_rst_map[1];
742  }
743  else if (gt == Geometry::PRISM)
744  {
745  gf_rst_map_temp = gf_rst_map[2];
746  }
747  else if (gt == Geometry::PYRAMID)
748  {
749  gf_rst_map_temp = gf_rst_map[3];
750  }
751 
752  int local_elem = split_element_index[elem];
753  Vector mfem_ref(dim);
754  // map to rst of macro element
755  gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
756 
757  for (int d = 0; d < dim; d++)
758  {
759  pt->r[d] = mfem_ref(d);
760  }
761  ++pt;
762  }
763 
764  // Transfer data back to source MPI rank
765  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
766  npt = outpt->n;
767 
768  // First copy mapped information for points on other procs
769  pt = (struct out_pt *)outpt->ptr;
770  for (int index = 0; index < npt; index++)
771  {
772  gsl_mfem_elem[pt->index] = pt->el;
773  for (int d = 0; d < dim; d++)
774  {
775  gsl_mfem_ref(d + pt->index*dim) = pt->r[d];
776  }
777  ++pt;
778  }
779  array_free(outpt);
780  delete outpt;
781 
782  // Now map information for points on the same proc
783  for (int index = 0; index < nptorig; index++)
784  {
785  if (gsl_code[index] != 2 && gsl_proc[index] == gsl_comm->id)
786  {
787  const int elem = gsl_elem[index];
788  const int mesh_elem = split_element_map[elem];
789  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(mesh_elem);
790  const Geometry::Type gt = fe->GetGeomType();
791  gsl_mfem_elem[index] = mesh_elem;
792  if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { continue; }
793  else if (gt == Geometry::TRIANGLE)
794  {
795  gf_rst_map_temp = gf_rst_map[0];
796  }
797  else if (gt == Geometry::TETRAHEDRON)
798  {
799  gf_rst_map_temp = gf_rst_map[1];
800  }
801  else if (gt == Geometry::PRISM)
802  {
803  gf_rst_map_temp = gf_rst_map[2];
804  }
805  else if (gt == Geometry::PYRAMID)
806  {
807  gf_rst_map_temp = gf_rst_map[3];
808  }
809 
810  int local_elem = split_element_index[elem];
811  IntegrationPoint ip;
812  Vector mfem_ref(gsl_mfem_ref.GetData()+index*dim, dim);
813  ip.Set2(mfem_ref.GetData());
814  if (dim == 3) { ip.z = mfem_ref(2); }
815  gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
816  }
817  }
818 }
819 
821  Vector &field_out)
822 {
823  const int gf_order = field_in.FESpace()->GetMaxElementOrder(),
824  mesh_order = mesh->GetNodalFESpace()->GetMaxElementOrder();
825 
826  const FiniteElementCollection *fec_in = field_in.FESpace()->FEColl();
827  const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>(fec_in);
828  const L2_FECollection *fec_l2 = dynamic_cast<const L2_FECollection *>(fec_in);
829 
830  if (fec_h1 && gf_order == mesh_order &&
831  fec_h1->GetBasisType() == BasisType::GaussLobatto &&
832  field_in.FESpace()->IsVariableOrder() ==
834  {
835  InterpolateH1(field_in, field_out);
836  return;
837  }
838  else
839  {
840  InterpolateGeneral(field_in, field_out);
841  if (!fec_l2 || avgtype == AvgType::NONE) { return; }
842  }
843 
844  // For points on element borders, project the L2 GridFunction to H1 and
845  // re-interpolate.
846  if (fec_l2)
847  {
848  Array<int> indl2;
849  for (int i = 0; i < points_cnt; i++)
850  {
851  if (gsl_code[i] == 1) { indl2.Append(i); }
852  }
853  int borderPts = indl2.Size();
854 #ifdef MFEM_USE_MPI
855  MPI_Allreduce(MPI_IN_PLACE, &borderPts, 1, MPI_INT, MPI_SUM, gsl_comm->c);
856 #endif
857  if (borderPts == 0) { return; } // no points on element borders
858 
859  Vector field_out_l2(field_out.Size());
860  VectorGridFunctionCoefficient field_in_dg(&field_in);
861  int gf_order_h1 = std::max(gf_order, 1); // H1 should be at least order 1
862  H1_FECollection fec(gf_order_h1, dim);
863  const int ncomp = field_in.FESpace()->GetVDim();
864  FiniteElementSpace fes(mesh, &fec, ncomp);
865  GridFunction field_in_h1(&fes);
866 
867  if (avgtype == AvgType::ARITHMETIC)
868  {
869  field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
870  }
871  else if (avgtype == AvgType::HARMONIC)
872  {
873  field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
874  }
875  else
876  {
877  MFEM_ABORT("Invalid averaging type.");
878  }
879 
880  if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
881  {
882  InterpolateH1(field_in_h1, field_out_l2);
883  }
884  else
885  {
886  InterpolateGeneral(field_in_h1, field_out_l2);
887  }
888 
889  // Copy interpolated values for the points on element border
890  for (int j = 0; j < ncomp; j++)
891  {
892  for (int i = 0; i < indl2.Size(); i++)
893  {
894  int idx = field_in.FESpace()->GetOrdering() == Ordering::byNODES ?
895  indl2[i] + j*points_cnt:
896  indl2[i]*ncomp + j;
897  field_out(idx) = field_out_l2(idx);
898  }
899  }
900  }
901 }
902 
904  Vector &field_out)
905 {
906  FiniteElementSpace ind_fes(mesh, field_in.FESpace()->FEColl());
907  if (field_in.FESpace()->IsVariableOrder())
908  {
909  for (int e = 0; e < ind_fes.GetMesh()->GetNE(); e++)
910  {
911  ind_fes.SetElementOrder(e, field_in.FESpace()->GetElementOrder(e));
912  }
913  ind_fes.Update(false);
914  }
915  GridFunction field_in_scalar(&ind_fes);
916  Vector node_vals;
917 
918  const int ncomp = field_in.FESpace()->GetVDim(),
919  points_fld = field_in.Size() / ncomp;
920  MFEM_VERIFY(points_cnt == gsl_code.Size(),
921  "FindPointsGSLIB::InterpolateH1: Inconsistent size of gsl_code");
922 
923  field_out.SetSize(points_cnt*ncomp);
924  field_out = default_interp_value;
925 
926  for (int i = 0; i < ncomp; i++)
927  {
928  const int dataptrin = i*points_fld,
929  dataptrout = i*points_cnt;
930  if (field_in.FESpace()->GetOrdering() == Ordering::byNODES)
931  {
932  field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin, points_fld);
933  }
934  else
935  {
936  for (int j = 0; j < points_fld; j++)
937  {
938  field_in_scalar(j) = field_in(i + j*ncomp);
939  }
940  }
941  GetNodalValues(&field_in_scalar, node_vals);
942 
943  if (dim==2)
944  {
945  findpts_eval_2(field_out.GetData()+dataptrout, sizeof(double),
946  gsl_code.GetData(), sizeof(unsigned int),
947  gsl_proc.GetData(), sizeof(unsigned int),
948  gsl_elem.GetData(), sizeof(unsigned int),
949  gsl_ref.GetData(), sizeof(double) * dim,
950  points_cnt, node_vals.GetData(), fdata2D);
951  }
952  else
953  {
954  findpts_eval_3(field_out.GetData()+dataptrout, sizeof(double),
955  gsl_code.GetData(), sizeof(unsigned int),
956  gsl_proc.GetData(), sizeof(unsigned int),
957  gsl_elem.GetData(), sizeof(unsigned int),
958  gsl_ref.GetData(), sizeof(double) * dim,
959  points_cnt, node_vals.GetData(), fdata3D);
960  }
961  }
962  if (field_in.FESpace()->GetOrdering() == Ordering::byVDIM)
963  {
964  Vector field_out_temp = field_out;
965  for (int i = 0; i < ncomp; i++)
966  {
967  for (int j = 0; j < points_cnt; j++)
968  {
969  field_out(i + j*ncomp) = field_out_temp(j + i*points_cnt);
970  }
971  }
972  }
973 }
974 
976  Vector &field_out)
977 {
978  int ncomp = field_in.VectorDim(),
979  nptorig = points_cnt,
980  npt = points_cnt;
981 
982  field_out.SetSize(points_cnt*ncomp);
983  field_out = default_interp_value;
984 
985  if (gsl_comm->np == 1) // serial
986  {
987  for (int index = 0; index < npt; index++)
988  {
989  if (gsl_code[index] == 2) { continue; }
990  IntegrationPoint ip;
992  if (dim == 3) { ip.z = gsl_mfem_ref(index*dim + 2); }
993  Vector localval(ncomp);
994  field_in.GetVectorValue(gsl_mfem_elem[index], ip, localval);
995  if (field_in.FESpace()->GetOrdering() == Ordering::byNODES)
996  {
997  for (int i = 0; i < ncomp; i++)
998  {
999  field_out(index + i*npt) = localval(i);
1000  }
1001  }
1002  else //byVDIM
1003  {
1004  for (int i = 0; i < ncomp; i++)
1005  {
1006  field_out(index*ncomp + i) = localval(i);
1007  }
1008  }
1009  }
1010  }
1011  else // parallel
1012  {
1013  // Determine number of points to be sent
1014  int nptsend = 0;
1015  for (int index = 0; index < npt; index++)
1016  {
1017  if (gsl_code[index] != 2) { nptsend +=1; }
1018  }
1019 
1020  // Pack data to send via crystal router
1021  struct gslib::array *outpt = new gslib::array;
1022  struct out_pt { double r[3], ival; uint index, el, proc; };
1023  struct out_pt *pt;
1024  array_init(struct out_pt, outpt, nptsend);
1025  outpt->n=nptsend;
1026  pt = (struct out_pt *)outpt->ptr;
1027  for (int index = 0; index < npt; index++)
1028  {
1029  if (gsl_code[index] == 2) { continue; }
1030  for (int d = 0; d < dim; ++d) { pt->r[d]= gsl_mfem_ref(index*dim + d); }
1031  pt->index = index;
1032  pt->proc = gsl_proc[index];
1033  pt->el = gsl_mfem_elem[index];
1034  ++pt;
1035  }
1036 
1037  // Transfer data to target MPI ranks
1038  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
1039 
1040  if (ncomp == 1)
1041  {
1042  // Interpolate the grid function
1043  npt = outpt->n;
1044  pt = (struct out_pt *)outpt->ptr;
1045  for (int index = 0; index < npt; index++)
1046  {
1047  IntegrationPoint ip;
1048  ip.Set3(&pt->r[0]);
1049  pt->ival = field_in.GetValue(pt->el, ip, 1);
1050  ++pt;
1051  }
1052 
1053  // Transfer data back to source MPI rank
1054  sarray_transfer(struct out_pt, outpt, proc, 1, cr);
1055  npt = outpt->n;
1056  pt = (struct out_pt *)outpt->ptr;
1057  for (int index = 0; index < npt; index++)
1058  {
1059  field_out(pt->index) = pt->ival;
1060  ++pt;
1061  }
1062  array_free(outpt);
1063  delete outpt;
1064  }
1065  else // ncomp > 1
1066  {
1067  // Interpolate data and store in a Vector
1068  npt = outpt->n;
1069  pt = (struct out_pt *)outpt->ptr;
1070  Vector vec_int_vals(npt*ncomp);
1071  for (int index = 0; index < npt; index++)
1072  {
1073  IntegrationPoint ip;
1074  ip.Set3(&pt->r[0]);
1075  Vector localval(vec_int_vals.GetData()+index*ncomp, ncomp);
1076  field_in.GetVectorValue(pt->el, ip, localval);
1077  ++pt;
1078  }
1079 
1080  // Save index and proc data in a struct
1081  struct gslib::array *savpt = new gslib::array;
1082  struct sav_pt { uint index, proc; };
1083  struct sav_pt *spt;
1084  array_init(struct sav_pt, savpt, npt);
1085  savpt->n=npt;
1086  spt = (struct sav_pt *)savpt->ptr;
1087  pt = (struct out_pt *)outpt->ptr;
1088  for (int index = 0; index < npt; index++)
1089  {
1090  spt->index = pt->index;
1091  spt->proc = pt->proc;
1092  ++pt; ++spt;
1093  }
1094 
1095  array_free(outpt);
1096  delete outpt;
1097 
1098  // Copy data from save struct to send struct and send component wise
1099  struct gslib::array *sendpt = new gslib::array;
1100  struct send_pt { double ival; uint index, proc; };
1101  struct send_pt *sdpt;
1102  for (int j = 0; j < ncomp; j++)
1103  {
1104  array_init(struct send_pt, sendpt, npt);
1105  sendpt->n=npt;
1106  spt = (struct sav_pt *)savpt->ptr;
1107  sdpt = (struct send_pt *)sendpt->ptr;
1108  for (int index = 0; index < npt; index++)
1109  {
1110  sdpt->index = spt->index;
1111  sdpt->proc = spt->proc;
1112  sdpt->ival = vec_int_vals(j + index*ncomp);
1113  ++sdpt; ++spt;
1114  }
1115 
1116  sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
1117  sdpt = (struct send_pt *)sendpt->ptr;
1118  for (int index = 0; index < static_cast<int>(sendpt->n); index++)
1119  {
1120  int idx = field_in.FESpace()->GetOrdering() == Ordering::byNODES ?
1121  sdpt->index + j*nptorig :
1122  sdpt->index*ncomp + j;
1123  field_out(idx) = sdpt->ival;
1124  ++sdpt;
1125  }
1126  array_free(sendpt);
1127  }
1128  array_free(savpt);
1129  delete sendpt;
1130  delete savpt;
1131  } // ncomp > 1
1132  } // parallel
1133 }
1134 
1135 void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
1136  GridFunction *gfmax,
1137  const double bb_t, const double newt_tol,
1138  const int npt_max)
1139 {
1140  MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
1141  const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
1142 
1143  // FreeData if OversetFindPointsGSLIB::Setup has been called already
1144  if (setupflag) { FreeData(); }
1145 
1146  crystal_init(cr, gsl_comm);
1147  mesh = &m;
1148  dim = mesh->Dimension();
1149  const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(0);
1150  unsigned dof1D = fe->GetOrder() + 1;
1151 
1152  SetupSplitMeshes();
1153  if (dim == 2)
1154  {
1155  if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
1156  ir_split[0] = new IntegrationRule(3*pow(dof1D, dim));
1158 
1159  if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
1160  ir_split[1] = new IntegrationRule(pow(dof1D, dim));
1162  }
1163  else if (dim == 3)
1164  {
1165  if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
1166  ir_split[0] = new IntegrationRule(pow(dof1D, dim));
1168 
1169  if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
1170  ir_split[1] = new IntegrationRule(4*pow(dof1D, dim));
1172 
1173  if (ir_split[2]) { delete ir_split[2]; ir_split[2] = NULL; }
1174  ir_split[2] = new IntegrationRule(3*pow(dof1D, dim));
1176 
1177  if (ir_split[3]) { delete ir_split[3]; ir_split[3] = NULL; }
1178  ir_split[3] = new IntegrationRule(8*pow(dof1D, dim));
1180  }
1181 
1183 
1184  MFEM_ASSERT(meshid>=0, " The ID should be greater than or equal to 0.");
1185 
1186  const int pts_cnt = gsl_mesh.Size()/dim,
1187  NEtot = pts_cnt/(int)pow(dof1D, dim);
1188 
1189  distfint.SetSize(pts_cnt);
1190  if (!gfmax)
1191  {
1192  distfint = 0.0;
1193  }
1194  else
1195  {
1196  GetNodalValues(gfmax, distfint);
1197  }
1198  u_meshid = (unsigned int)meshid;
1199 
1200  if (dim == 2)
1201  {
1202  unsigned nr[2] = { dof1D, dof1D };
1203  unsigned mr[2] = { 2*dof1D, 2*dof1D };
1204  double * const elx[2] = { &gsl_mesh(0), &gsl_mesh(pts_cnt) };
1205  fdata2D = findptsms_setup_2(gsl_comm, elx, nr, NEtot, mr, bb_t,
1206  pts_cnt, pts_cnt, npt_max, newt_tol,
1207  &u_meshid, &distfint(0));
1208  }
1209  else
1210  {
1211  unsigned nr[3] = { dof1D, dof1D, dof1D };
1212  unsigned mr[3] = { 2*dof1D, 2*dof1D, 2*dof1D };
1213  double * const elx[3] =
1214  { &gsl_mesh(0), &gsl_mesh(pts_cnt), &gsl_mesh(2*pts_cnt) };
1215  fdata3D = findptsms_setup_3(gsl_comm, elx, nr, NEtot, mr, bb_t,
1216  pts_cnt, pts_cnt, npt_max, newt_tol,
1217  &u_meshid, &distfint(0));
1218  }
1219  setupflag = true;
1220  overset = true;
1221 }
1222 
1224  Array<unsigned int> &point_id,
1225  int point_pos_ordering)
1226 {
1227  MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
1228  "finding points.");
1229  MFEM_VERIFY(overset, " Please setup FindPoints for overlapping grids.");
1230  points_cnt = point_pos.Size() / dim;
1231  unsigned int match = 0; // Don't find points in the mesh if point_id = mesh_id
1232 
1238 
1239  auto xvFill = [&](const double *xv_base[], unsigned xv_stride[])
1240  {
1241  for (int d = 0; d < dim; d++)
1242  {
1243  if (point_pos_ordering == Ordering::byNODES)
1244  {
1245  xv_base[d] = point_pos.GetData() + d*points_cnt;
1246  xv_stride[d] = sizeof(double);
1247  }
1248  else
1249  {
1250  xv_base[d] = point_pos.GetData() + d;
1251  xv_stride[d] = dim*sizeof(double);
1252  }
1253  }
1254  };
1255  if (dim == 2)
1256  {
1257  const double *xv_base[2];
1258  unsigned xv_stride[2];
1259  xvFill(xv_base, xv_stride);
1260  findptsms_2(gsl_code.GetData(), sizeof(unsigned int),
1261  gsl_proc.GetData(), sizeof(unsigned int),
1262  gsl_elem.GetData(), sizeof(unsigned int),
1263  gsl_ref.GetData(), sizeof(double) * dim,
1264  gsl_dist.GetData(), sizeof(double),
1265  xv_base, xv_stride,
1266  point_id.GetData(), sizeof(unsigned int), &match,
1267  points_cnt, fdata2D);
1268  }
1269  else // dim == 3
1270  {
1271  const double *xv_base[3];
1272  unsigned xv_stride[3];
1273  xvFill(xv_base, xv_stride);
1274  findptsms_3(gsl_code.GetData(), sizeof(unsigned int),
1275  gsl_proc.GetData(), sizeof(unsigned int),
1276  gsl_elem.GetData(), sizeof(unsigned int),
1277  gsl_ref.GetData(), sizeof(double) * dim,
1278  gsl_dist.GetData(), sizeof(double),
1279  xv_base, xv_stride,
1280  point_id.GetData(), sizeof(unsigned int), &match,
1281  points_cnt, fdata3D);
1282  }
1283 
1284  // Set the element number and reference position to 0 for points not found
1285  for (int i = 0; i < points_cnt; i++)
1286  {
1287  if (gsl_code[i] == 2 ||
1288  (gsl_code[i] == 1 && gsl_dist(i) > bdr_tol))
1289  {
1290  gsl_elem[i] = 0;
1291  for (int d = 0; d < dim; d++) { gsl_ref(i*dim + d) = -1.; }
1292  gsl_code[i] = 2;
1293  }
1294  }
1295 
1296  // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for both
1297  // simplices and quads.
1299 }
1300 
1302  Array<unsigned int> &point_id,
1303  const GridFunction &field_in,
1304  Vector &field_out,
1305  int point_pos_ordering)
1306 {
1307  FindPoints(point_pos, point_id, point_pos_ordering);
1308  Interpolate(field_in, field_out);
1309 }
1310 
1311 
1312 } // namespace mfem
1313 
1314 #endif // MFEM_USE_GSLIB
Abstract class for all finite elements.
Definition: fe_base.hpp:233
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:820
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Array< unsigned int > gsl_elem
Definition: gslib.hpp:72
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:5630
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:577
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:160
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:476
virtual ~FindPointsGSLIB()
Definition: gslib.cpp:67
void FindPoints(const Vector &point_pos, int point_pos_ordering=Ordering::byNODES)
Definition: gslib.cpp:176
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Array< GridFunction * > gf_rst_map
Definition: gslib.hpp:65
void SetElementOrder(int i, int p)
Sets the order of the i&#39;th finite element.
Definition: fespace.cpp:151
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, double sx=1.0, double sy=1.0, double sz=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3783
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2589
Array< unsigned int > gsl_mfem_elem
Definition: gslib.hpp:72
T * GetData()
Returns the data.
Definition: array.hpp:115
unsigned int uint
Definition: gecko.hpp:204
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void SetupSplitMeshes()
Definition: gslib.cpp:301
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void Interpolate(const Vector &point_pos, Array< unsigned int > &point_id, const GridFunction &field_in, Vector &field_out, int point_pos_ordering=Ordering::byNODES)
Definition: gslib.cpp:1301
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
int GetBasisType() const
Definition: fe_coll.hpp:282
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
Array< unsigned int > gsl_proc
Definition: gslib.hpp:72
Array< Mesh * > mesh_split
Definition: gslib.hpp:59
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:449
void DeleteAll()
Delete the whole array.
Definition: array.hpp:854
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:2841
const Array< int > & GetDofMap() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:1210
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
virtual void MapRefPosAndElemIndices()
Definition: gslib.cpp:673
struct gslib::findpts_data_3 * fdata3D
Definition: gslib.hpp:68
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
virtual void GetNodalValues(const GridFunction *gf_in, Vector &node_vals)
Get GridFunction value at the points expected by GSLIB.
Definition: gslib.cpp:572
Array< IntegrationRule * > ir_split
Definition: gslib.hpp:62
int GetMaxElementOrder() const
Return the maximum polynomial order.
Definition: fespace.hpp:573
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
virtual void FreeData()
Definition: gslib.cpp:272
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
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:107
double default_interp_value
Definition: gslib.hpp:75
struct gslib::crystal * cr
Definition: gslib.hpp:69
void Set2(const double x1, const double x2)
Definition: intrules.hpp:86
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5577
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
Array< int > split_element_index
Definition: gslib.hpp:78
Array< FiniteElementSpace * > fes_rst_map
Definition: gslib.hpp:64
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:76
Array< unsigned int > gsl_code
Definition: gslib.hpp:72
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
virtual void InterpolateGeneral(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:975
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
virtual void SetupIntegrationRuleForSplitMesh(Mesh *mesh, IntegrationRule *irule, int order)
Definition: gslib.cpp:524
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3773
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
struct gslib::findpts_data_2 * fdata2D
Definition: gslib.hpp:67
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Class for integration point with weight.
Definition: intrules.hpp:31
int VectorDim() const
Definition: gridfunc.cpp:323
int GetElementOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:178
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
Array< int > split_element_map
Definition: gslib.hpp:77
int dim
Definition: ex24.cpp:53
void Destroy()
Destroy a vector.
Definition: vector.hpp:594
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
struct gslib::comm * gsl_comm
Definition: gslib.hpp:70
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
FiniteElementCollection * fec_map_lin
Definition: gslib.hpp:66
Vector data type.
Definition: vector.hpp:58
Vector coefficient defined by a vector GridFunction.
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
Definition: gslib.cpp:23
void FindPoints(const Vector &point_pos, Array< unsigned int > &point_id, int point_pos_ordering=Ordering::byNODES)
Definition: gslib.cpp:1223
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:1135
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
virtual void InterpolateH1(const GridFunction &field_in, Vector &field_out)
Use GSLIB for communication and interpolation.
Definition: gslib.cpp:903
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327