MFEM  v4.1.0 Finite element discretization library
pfindpts.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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 // Find Points Miniapp: Evaluate grid function in physical space - Parallel Ver.
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-q3.mesh -o 3
31 // mpirun -np 2 pfindpts -m ../../data/fichera.mesh -o 3
32
33 #include "mfem.hpp"
34
35 using namespace mfem;
36 using namespace std;
37
38 // Scalar function to project
39 double field_func(const Vector &x)
40 {
41  const int dim = x.Size();
42  double res = 0.0;
43  for (int d = 0; d < dim; d++) { res += x(d) * x(d); }
44  return res;
45 }
46
47 int main (int argc, char *argv[])
48 {
49  // Initialize MPI.
50  int num_procs, myid;
51  MPI_Init(&argc, &argv);
52  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
53  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
54
55  // Set the method's default parameters.
56  const char *mesh_file = "../../data/rt-2d-q3.mesh";
57  int mesh_poly_deg = 3;
58  int rs_levels = 0;
59  int rp_levels = 0;
60  bool visualization = true;
61
62  // Parse command-line options.
63  OptionsParser args(argc, argv);
65  "Mesh file to use.");
67  "Polynomial degree of mesh finite element space.");
69  "Number of times to refine the mesh uniformly in serial.");
71  "Number of times to refine the mesh uniformly in parallel.");
73  "--no-visualization",
74  "Enable or disable GLVis visualization.");
75  args.Parse();
76  if (!args.Good())
77  {
78  args.PrintUsage(cout);
79  return 1;
80  }
81  if (myid == 0) { args.PrintOptions(cout); }
82
83  // Initialize and refine the starting mesh.
84  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
85  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
86  const int dim = mesh->Dimension();
87  if (myid == 0)
88  {
89  cout << "Mesh curvature of the original mesh: ";
90  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
91  else { cout << "(NONE)"; }
92  cout << endl;
93  }
94
95  // Mesh bounding box (for the full serial mesh).
96  Vector pos_min, pos_max;
97  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
98  mesh->GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
99  if (myid == 0)
100  {
101  cout << "--- Generating equidistant point for:\n"
102  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
103  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
104  if (dim == 3)
105  {
106  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
107  }
108  }
109
110  // Distribute the mesh.
111  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
112  delete mesh;
113  for (int lev = 0; lev < rp_levels; lev++) { pmesh.UniformRefinement(); }
114
115  // Curve the mesh based on the chosen polynomial degree.
116  H1_FECollection fec(mesh_poly_deg, dim);
117  ParFiniteElementSpace pfespace(&pmesh, &fec, dim);
118  pmesh.SetNodalFESpace(&pfespace);
119  if (myid == 0)
120  {
121  cout << "Mesh curvature of the curved mesh: " << fec.Name() << endl;
122  }
123
124  // Define a scalar function on the mesh.
125  ParFiniteElementSpace sc_fes(&pmesh, &fec, 1);
126  GridFunction field_vals(&sc_fes);
128  field_vals.ProjectCoefficient(fc);
129
130  // Display the mesh and the field through glvis.
131  if (visualization)
132  {
133  char vishost[] = "localhost";
134  int visport = 19916;
135  socketstream sout;
136  sout.open(vishost, visport);
137  if (!sout)
138  {
139  if (myid == 0)
140  {
141  cout << "Unable to connect to GLVis server at "
142  << vishost << ':' << visport << endl;
143  }
144  }
145  else
146  {
147  sout << "parallel " << num_procs << " " << myid << "\n";
148  sout.precision(8);
149  sout << "solution\n" << pmesh << field_vals;
150  if (dim == 2) { sout << "keys RmjA*****\n"; }
151  if (dim == 3) { sout << "keys mA\n"; }
152  sout << flush;
153  }
154  }
155
156  // Setup the gslib mesh.
157  FindPointsGSLIB finder(MPI_COMM_WORLD);
158  const double rel_bbox_el = 0.05;
159  const double newton_tol = 1.0e-12;
160  const int npts_at_once = 256;
161  finder.Setup(pmesh, rel_bbox_el, newton_tol, npts_at_once);
162
163  // Generate equidistant points in physical coordinates over the whole mesh.
164  // Note that some points might be outside, if the mesh is not a box. Note
165  // also that all tasks search the same points (not mandatory).
166  const int pts_cnt_1D = 5;
167  const int pts_cnt = pow(pts_cnt_1D, dim);
168  Vector vxyz(pts_cnt * dim);
169  if (dim == 2)
170  {
172  const IntegrationRule &ir = el.GetNodes();
173  for (int i = 0; i < ir.GetNPoints(); i++)
174  {
175  const IntegrationPoint &ip = ir.IntPoint(i);
176  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
177  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
178  }
179  }
180  else
181  {
183  const IntegrationRule &ir = el.GetNodes();
184  for (int i = 0; i < ir.GetNPoints(); i++)
185  {
186  const IntegrationPoint &ip = ir.IntPoint(i);
187  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
188  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
189  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
190  }
191  }
192
193  Array<unsigned int> el_id_out(pts_cnt), code_out(pts_cnt),
195  Vector pos_r_out(pts_cnt * dim), dist_p_out(pts_cnt);
196
197  // Finds points stored in vxyz.
199  el_id_out, pos_r_out, dist_p_out);
200
201  // Interpolate FE function values on the found points.
202  Vector interp_vals(pts_cnt);
204  pos_r_out, field_vals, interp_vals);
205
206  // Free the internal gslib data.
207  finder.FreeData();
208
209  int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
210  double max_err = 0.0, max_dist = 0.0;
211  Vector pos(dim);
212  for (int i = 0; i < pts_cnt; i++)
213  {
214  (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
215
216  if (code_out[i] < 2)
217  {
218  for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
219  const double exact_val = field_func(pos);
220
221  max_err = std::max(max_err, fabs(exact_val - interp_vals[i]));
222  max_dist = std::max(max_dist, dist_p_out(i));
223  if (code_out[i] == 1) { face_pts++; }
224  }
225  else { not_found++; }
226  }
227
228  // We print only the task 0 result (other tasks should be identical except
229  // the number of points found locally).
230  if (myid == 0)
231  {
232  cout << setprecision(16) << "--- Task " << myid << ": "
233  << "\nSearched points: " << pts_cnt
234  << "\nFound on local mesh: " << found_loc
235  << "\nFound on other tasks: " << found_away
236  << "\nMax interp error: " << max_err
237  << "\nMax dist (of found): " << max_dist
239  << "\nPoints on faces: " << face_pts << endl;
240  }
241
242  MPI_Finalize();
243  return 0;
244 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
void FindPoints(Vector &point_pos, Array< unsigned int > &codes, Array< unsigned int > &proc_ids, Array< unsigned int > &elem_ids, Vector &ref_pos, Vector &dist)
Definition: gslib.cpp:115
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:27
void Setup(Mesh &m, double bb_t, double newt_tol, int npt_max)
Definition: gslib.cpp:58
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:120
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
Definition: ex1.cpp:62
double field_func(const Vector &x)
Definition: findpts.cpp:39
void Interpolate(Array< unsigned int > &codes, Array< unsigned int > &proc_ids, Array< unsigned int > &elem_ids, Vector &ref_pos, const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:155
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
const IntegrationRule & GetNodes() const
Definition: fe.hpp:367
virtual const char * Name() const
Definition: fe_coll.hpp:104
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3885
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe.hpp:37
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
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)
Definition: optparser.hpp:76
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1685
int open(const char hostname[], int port)
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Class for parallel meshes.
Definition: pmesh.hpp:32
bool Good() const
Definition: optparser.hpp:125