MFEM  v4.5.2
Finite element discretization library
pfindpts.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 // ----------------------------------------------------------------------
13 // Parallel Find Points Miniapp: Evaluate grid function in physical space
14 // ----------------------------------------------------------------------
15 //
16 // This miniapp demonstrates the interpolation of a high-order grid function on
17 // a set of points in physical-space. The miniapp is based on GSLIB-FindPoints,
18 // which provides two key functionalities. First, for a given set of points in
19 // the physical-space, it determines the computational coordinates (element
20 // number, reference-space coordinates inside the element, and processor number
21 // [in parallel]) for each point. Second, based on computational coordinates, it
22 // interpolates a grid function in the given points. Inside GSLIB, computation
23 // of the coordinates requires use of a Hash Table to identify the candidate
24 // processor and element for each point, followed by the Newton's method to
25 // determine the reference-space coordinates inside the candidate element.
26 //
27 // Compile with: make pfindpts
28 //
29 // Sample runs:
30 // mpirun -np 2 pfindpts -m ../../data/rt-2d-p4-tri.mesh -o 4
31 // mpirun -np 2 pfindpts -m ../../data/inline-tri.mesh -o 3
32 // mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -o 3
33 // mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -o 3 -po 1
34 // mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -o 3 -po 1 -gfo 1 -nc 2
35 // mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -o 3 -hr
36 // mpirun -np 2 pfindpts -m ../../data/inline-tet.mesh -o 3
37 // mpirun -np 2 pfindpts -m ../../data/inline-hex.mesh -o 3
38 // mpirun -np 2 pfindpts -m ../../data/inline-wedge.mesh -o 3
39 // mpirun -np 2 pfindpts -m ../../data/amr-quad.mesh -o 2
40 // mpirun -np 2 pfindpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
41 // mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -ft 1 -no-vis -sr0
42 // mpirun -np 2 pfindpts -m ../../data/square-mixed.mesh -o 2 -mo 2
43 // mpirun -np 2 pfindpts -m ../../data/square-mixed.mesh -o 2 -mo 2 -hr
44 // mpirun -np 2 pfindpts -m ../../data/square-mixed.mesh -o 2 -mo 3 -ft 2
45 // mpirun -np 2 pfindpts -m ../../data/fichera-mixed.mesh -o 3 -mo 2
46 // mpirun -np 2 pfindpts -m ../../data/inline-pyramid.mesh -o 1 -mo 1
47 // mpirun -np 2 pfindpts -m ../../data/tinyzoo-3d.mesh -o 1 -mo 1
48 
49 #include "mfem.hpp"
50 
51 using namespace mfem;
52 using namespace std;
53 
54 double func_order;
55 
56 // Scalar function to project
57 double field_func(const Vector &x)
58 {
59  const int dim = x.Size();
60  double res = 0.0;
61  for (int d = 0; d < dim; d++) { res += std::pow(x(d), func_order); }
62  return res;
63 }
64 
65 void F_exact(const Vector &p, Vector &F)
66 {
67  F(0) = field_func(p);
68  for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*F(0); }
69 }
70 
71 int main (int argc, char *argv[])
72 {
73  // Initialize MPI and HYPRE.
74  Mpi::Init(argc, argv);
75  int num_procs = Mpi::WorldSize();
76  int myid = Mpi::WorldRank();
77  Hypre::Init();
78 
79  // Set the method's default parameters.
80  const char *mesh_file = "../../data/rt-2d-q3.mesh";
81  int order = 3;
82  int mesh_poly_deg = 3;
83  int rs_levels = 0;
84  int rp_levels = 0;
85  bool visualization = true;
86  int fieldtype = 0;
87  int ncomp = 1;
88  bool search_on_rank_0 = false;
89  bool hrefinement = false;
90  int point_ordering = 0;
91  int gf_ordering = 0;
92 
93  // Parse command-line options.
94  OptionsParser args(argc, argv);
95  args.AddOption(&mesh_file, "-m", "--mesh",
96  "Mesh file to use.");
97  args.AddOption(&order, "-o", "--order",
98  "Finite element order (polynomial degree).");
99  args.AddOption(&mesh_poly_deg, "-mo", "--mesh-order",
100  "Polynomial degree of mesh finite element space.");
101  args.AddOption(&rs_levels, "-rs", "--refine-serial",
102  "Number of times to refine the mesh uniformly in serial.");
103  args.AddOption(&rp_levels, "-rp", "--refine-parallel",
104  "Number of times to refine the mesh uniformly in parallel.");
105  args.AddOption(&fieldtype, "-ft", "--field-type",
106  "Field type: 0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
107  args.AddOption(&ncomp, "-nc", "--ncomp",
108  "Number of components for H1 or L2 GridFunctions");
109  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
110  "--no-visualization",
111  "Enable or disable GLVis visualization.");
112  args.AddOption(&search_on_rank_0, "-sr0", "--search-on-r0", "-no-sr0",
113  "--no-search-on-r0",
114  "Enable search only on rank 0 (disable to search points on all tasks).");
115  args.AddOption(&hrefinement, "-hr", "--h-refinement", "-no-hr",
116  "--no-h-refinement",
117  "Do random h refinements to mesh (does not work for pyramids).");
118  args.AddOption(&point_ordering, "-po", "--point-ordering",
119  "Ordering of points to be found."
120  "0 (default): byNodes, 1: byVDIM");
121  args.AddOption(&gf_ordering, "-gfo", "--gridfunc-ordering",
122  "Ordering of fespace that will be used for gridfunction to be interpolated."
123  "0 (default): byNodes, 1: byVDIM");
124 
125  args.Parse();
126  if (!args.Good())
127  {
128  args.PrintUsage(cout);
129  return 1;
130  }
131  if (myid == 0) { args.PrintOptions(cout); }
132 
133  func_order = std::min(order, 2);
134 
135  // Initialize and refine the starting mesh.
136  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
137  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
138  const int dim = mesh->Dimension();
139  if (myid == 0)
140  {
141  cout << "Mesh curvature of the original mesh: ";
142  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
143  else { cout << "(NONE)"; }
144  cout << endl;
145  }
146 
147  // Mesh bounding box (for the full serial mesh).
148  Vector pos_min, pos_max;
149  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
150  mesh->GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
151  if (myid == 0)
152  {
153  cout << "--- Generating equidistant point for:\n"
154  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
155  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
156  if (dim == 3)
157  {
158  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
159  }
160  }
161 
162  // Distribute the mesh.
163  if (hrefinement) { mesh->EnsureNCMesh(); }
164  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
165  delete mesh;
166  for (int lev = 0; lev < rp_levels; lev++) { pmesh.UniformRefinement(); }
167 
168  // Random h-refinements to mesh
169  if (hrefinement) { pmesh.RandomRefinement(0.5); }
170 
171  // Curve the mesh based on the chosen polynomial degree.
172  H1_FECollection fecm(mesh_poly_deg, dim);
173  ParFiniteElementSpace pfespace(&pmesh, &fecm, dim);
174  pmesh.SetNodalFESpace(&pfespace);
175  if (myid == 0)
176  {
177  cout << "Mesh curvature of the curved mesh: " << fecm.Name() << endl;
178  }
179 
180  MFEM_VERIFY(ncomp > 0, "Invalid number of components.");
181  int vec_dim = ncomp;
182  FiniteElementCollection *fec = NULL;
183  if (fieldtype == 0)
184  {
185  fec = new H1_FECollection(order, dim);
186  if (myid == 0) { cout << "H1-GridFunction\n"; }
187  }
188  else if (fieldtype == 1)
189  {
190  fec = new L2_FECollection(order, dim);
191  if (myid == 0) { cout << "L2-GridFunction\n"; }
192  }
193  else if (fieldtype == 2)
194  {
195  fec = new RT_FECollection(order, dim);
196  ncomp = 1;
197  vec_dim = dim;
198  if (myid == 0) { cout << "H(div)-GridFunction\n"; }
199  }
200  else if (fieldtype == 3)
201  {
202  fec = new ND_FECollection(order, dim);
203  ncomp = 1;
204  vec_dim = dim;
205  if (myid == 0) { cout << "H(curl)-GridFunction\n"; }
206  }
207  else
208  {
209  if (myid == 0) { MFEM_ABORT("Invalid FECollection type."); }
210  }
211  ParFiniteElementSpace sc_fes(&pmesh, fec, ncomp, gf_ordering);
212  ParGridFunction field_vals(&sc_fes);
213 
214  // Project the GridFunction using VectorFunctionCoefficient.
216  field_vals.ProjectCoefficient(F);
217 
218  // Display the mesh and the field through glvis.
219  if (visualization)
220  {
221  char vishost[] = "localhost";
222  int visport = 19916;
223  socketstream sout;
224  sout.open(vishost, visport);
225  if (!sout)
226  {
227  if (myid == 0)
228  {
229  cout << "Unable to connect to GLVis server at "
230  << vishost << ':' << visport << endl;
231  }
232  }
233  else
234  {
235  sout << "parallel " << num_procs << " " << myid << "\n";
236  sout.precision(8);
237  sout << "solution\n" << pmesh << field_vals;
238  if (dim == 2) { sout << "keys RmjA*****\n"; }
239  if (dim == 3) { sout << "keys mA\n"; }
240  sout << flush;
241  }
242  }
243 
244  // Generate equidistant points in physical coordinates over the whole mesh.
245  // Note that some points might be outside, if the mesh is not a box. Note
246  // also that all tasks search the same points (not mandatory).
247  const int pts_cnt_1D = 25;
248  int pts_cnt = pow(pts_cnt_1D, dim);
249  Vector vxyz(pts_cnt * dim);
250  if (dim == 2)
251  {
253  const IntegrationRule &ir = el.GetNodes();
254  for (int i = 0; i < ir.GetNPoints(); i++)
255  {
256  const IntegrationPoint &ip = ir.IntPoint(i);
257  if (point_ordering == Ordering::byNODES)
258  {
259  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
260  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
261  }
262  else
263  {
264  vxyz(i*dim + 0) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
265  vxyz(i*dim + 1) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
266  }
267  }
268  }
269  else
270  {
272  const IntegrationRule &ir = el.GetNodes();
273  for (int i = 0; i < ir.GetNPoints(); i++)
274  {
275  const IntegrationPoint &ip = ir.IntPoint(i);
276  if (point_ordering == Ordering::byNODES)
277  {
278  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
279  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
280  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
281  }
282  else
283  {
284  vxyz(i*dim + 0) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
285  vxyz(i*dim + 1) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
286  vxyz(i*dim + 2) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
287  }
288  }
289  }
290 
291  if ( (myid != 0) && (search_on_rank_0) )
292  {
293  pts_cnt = 0;
294  vxyz.Destroy();
295  }
296 
297  // Find and Interpolate FE function values on the desired points.
298  Vector interp_vals(pts_cnt*vec_dim);
299  FindPointsGSLIB finder(MPI_COMM_WORLD);
300  finder.Setup(pmesh);
301  finder.Interpolate(vxyz, field_vals, interp_vals, point_ordering);
302  Array<unsigned int> code_out = finder.GetCode();
303  Array<unsigned int> task_id_out = finder.GetProc();
304  Vector dist_p_out = finder.GetDist();
305 
306  // Print the results for task 0 since either 1) all tasks have the
307  // same set of points or 2) only task 0 has any points.
308  if (myid == 0 )
309  {
310  int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
311  double err = 0.0, max_err = 0.0, max_dist = 0.0;
312  Vector pos(dim);
313  int npt = 0;
314  for (int j = 0; j < vec_dim; j++)
315  {
316  for (int i = 0; i < pts_cnt; i++)
317  {
318  if (j == 0)
319  {
320  (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
321  }
322 
323  if (code_out[i] < 2)
324  {
325  for (int d = 0; d < dim; d++)
326  {
327  pos(d) = point_ordering == Ordering::byNODES ?
328  vxyz(d*pts_cnt + i) :
329  vxyz(i*dim + d);
330  }
331  Vector exact_val(vec_dim);
332  F_exact(pos, exact_val);
333  err = gf_ordering == Ordering::byNODES ?
334  fabs(exact_val(j) - interp_vals[i + j*pts_cnt]) :
335  fabs(exact_val(j) - interp_vals[i*vec_dim + j]);
336  max_err = std::max(max_err, err);
337  max_dist = std::max(max_dist, dist_p_out(i));
338  if (code_out[i] == 1 && j == 0) { face_pts++; }
339  }
340  else { if (j == 0) { not_found++; } }
341  npt++;
342  }
343  }
344 
345  cout << setprecision(16)
346  << "Searched unique points: " << pts_cnt
347  << "\nFound on local mesh: " << found_loc
348  << "\nFound on other tasks: " << found_away
349  << "\nMax interp error: " << max_err
350  << "\nMax dist (of found): " << max_dist
351  << "\nPoints not found: " << not_found
352  << "\nPoints on faces: " << face_pts << endl;
353  }
354 
355  // Free the internal gslib data.
356  finder.FreeData();
357 
358  delete fec;
359 
360  return 0;
361 }
const char vishost[]
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:803
virtual const Array< unsigned int > & GetProc() const
Return MPI rank on which each point was found by FindPoints.
Definition: gslib.hpp:210
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
double field_func(const Vector &x)
Definition: pfindpts.cpp:57
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Arbitrary order L2 elements in 3D on a cube.
Definition: fe_l2.hpp:78
int Dimension() const
Definition: mesh.hpp:1047
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
void F_exact(const Vector &p, Vector &F)
Definition: pfindpts.cpp:65
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
STL namespace.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
int main(int argc, char *argv[])
Definition: pfindpts.cpp:71
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
virtual void FreeData()
Definition: gslib.cpp:267
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9535
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
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
void RandomRefinement(double prob, bool aniso=false, int nonconforming=-1, int nc_limit=0)
Refine each element with given probability. Uses GeneralRefinement.
Definition: mesh.cpp:9559
static void Init()
Singleton creation with Mpi::Init();.
const int visport
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:373
A general vector function coefficient.
virtual const char * Name() const
Definition: fe_coll.hpp:273
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
double func_order
Definition: pfindpts.cpp:54
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
virtual const Vector & GetDist() const
Definition: gslib.hpp:215
Arbitrary order L2 elements in 2D on a square.
Definition: fe_l2.hpp:44
void Destroy()
Destroy a vector.
Definition: vector.hpp:589
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:447
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe_base.hpp:34
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:52
Vector data type.
Definition: vector.hpp:60
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2076
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual const Array< unsigned int > & GetCode() const
Definition: gslib.hpp:206
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388
Class for parallel meshes.
Definition: pmesh.hpp:32
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320