MFEM  v4.1.0 Finite element discretization library
findpts.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
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 findpts
28 //
29 // Sample runs:
30 // findpts -m ../../data/rt-2d-q3.mesh -o 3
31 // findpts -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  // Set the method's default parameters.
50  const char *mesh_file = "../../data/rt-2d-q3.mesh";
51  int mesh_poly_deg = 3;
52  int rs_levels = 0;
53  bool visualization = true;
54
55  // Parse command-line options.
56  OptionsParser args(argc, argv);
58  "Mesh file to use.");
60  "Polynomial degree of mesh finite element space.");
62  "Number of times to refine the mesh uniformly in serial.");
64  "--no-visualization",
65  "Enable or disable GLVis visualization.");
66  args.Parse();
67  if (!args.Good())
68  {
69  args.PrintUsage(cout);
70  return 1;
71  }
72  args.PrintOptions(cout);
73
74  // Initialize and refine the starting mesh.
75  Mesh mesh(mesh_file, 1, 1, false);
76  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
77  const int dim = mesh.Dimension();
78  cout << "Mesh curvature of the original mesh: ";
79  if (mesh.GetNodes()) { cout << mesh.GetNodes()->OwnFEC()->Name(); }
80  else { cout << "(NONE)"; }
81  cout << endl;
82
83  // Mesh bounding box.
84  Vector pos_min, pos_max;
85  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
86  mesh.GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
87  cout << "--- Generating equidistant point for:\n"
88  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
89  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
90  if (dim == 3)
91  {
92  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
93  }
94
95  // Curve the mesh based on the chosen polynomial degree.
96  H1_FECollection fec(mesh_poly_deg, dim);
97  FiniteElementSpace fespace(&mesh, &fec, dim);
98  mesh.SetNodalFESpace(&fespace);
99  cout << "Mesh curvature of the curved mesh: " << fec.Name() << endl;
100
101  // Define a scalar function on the mesh.
102  FiniteElementSpace sc_fes(&mesh, &fec, 1);
103  GridFunction field_vals(&sc_fes);
105  field_vals.ProjectCoefficient(fc);
106
107  // Display the mesh and the field through glvis.
108  if (visualization)
109  {
110  char vishost[] = "localhost";
111  int visport = 19916;
112  socketstream sout;
113  sout.open(vishost, visport);
114  if (!sout)
115  {
116  cout << "Unable to connect to GLVis server at "
117  << vishost << ':' << visport << endl;
118  }
119  else
120  {
121  sout.precision(8);
122  sout << "solution\n" << mesh << field_vals;
123  if (dim == 2) { sout << "keys RmjA*****\n"; }
124  if (dim == 3) { sout << "keys mA\n"; }
125  sout << flush;
126  }
127  }
128
129  // Setup the gslib mesh.
130  FindPointsGSLIB finder;
131  const double rel_bbox_el = 0.05;
132  const double newton_tol = 1.0e-12;
133  const int npts_at_once = 256;
134  finder.Setup(mesh, rel_bbox_el, newton_tol, npts_at_once);
135
136  // Generate equidistant points in physical coordinates over the whole mesh.
137  // Note that some points might be outside, if the mesh is not a box. Note
138  // also that all tasks search the same points (not mandatory).
139  const int pts_cnt_1D = 5;
140  const int pts_cnt = pow(pts_cnt_1D, dim);
141  Vector vxyz(pts_cnt * dim);
142  if (dim == 2)
143  {
145  const IntegrationRule &ir = el.GetNodes();
146  for (int i = 0; i < ir.GetNPoints(); i++)
147  {
148  const IntegrationPoint &ip = ir.IntPoint(i);
149  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
150  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
151  }
152  }
153  else
154  {
156  const IntegrationRule &ir = el.GetNodes();
157  for (int i = 0; i < ir.GetNPoints(); i++)
158  {
159  const IntegrationPoint &ip = ir.IntPoint(i);
160  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
161  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
162  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
163  }
164  }
165
166  Array<unsigned int> el_id_out(pts_cnt), code_out(pts_cnt),
168  Vector pos_r_out(pts_cnt * dim), dist_p_out(pts_cnt);
169
170  // Finds points stored in vxyz.
172  el_id_out, pos_r_out, dist_p_out);
173
174  // Interpolate FE function values on the found points.
175  Vector interp_vals(pts_cnt);
177  pos_r_out, field_vals, interp_vals);
178
179  // Free the internal gslib data.
180  finder.FreeData();
181
182  int face_pts = 0, not_found = 0, found = 0;
183  double max_err = 0.0, max_dist = 0.0;
184  Vector pos(dim);
185  for (int i = 0; i < pts_cnt; i++)
186  {
187  if (code_out[i] < 2)
188  {
189  found++;
190  for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
191  const double exact_val = field_func(pos);
192
193  max_err = std::max(max_err, fabs(exact_val - interp_vals[i]));
194  max_dist = std::max(max_dist, dist_p_out(i));
195  if (code_out[i] == 1) { face_pts++; }
196  }
197  else { not_found++; }
198  }
199
200  cout << setprecision(16)
201  << "Searched points: " << pts_cnt
202  << "\nFound points: " << found
203  << "\nMax interp error: " << max_err
204  << "\nMax dist (of found): " << max_dist
206  << "\nPoints on faces: " << face_pts << endl;
207
208  return 0;
209 }
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
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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
bool Good() const
Definition: optparser.hpp:125