MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pfindpts.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 // ----------------------------------------------------------------------
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-tet.mesh -o 3
34 // mpirun -np 2 pfindpts -m ../../data/inline-hex.mesh -o 3
35 // mpirun -np 2 pfindpts -m ../../data/inline-wedge.mesh -o 3
36 // mpirun -np 2 pfindpts -m ../../data/amr-quad.mesh -o 2
37 // mpirun -np 2 pfindpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
38 
39 
40 #include "mfem.hpp"
41 
42 using namespace mfem;
43 using namespace std;
44 
45 // Scalar function to project
46 double field_func(const Vector &x)
47 {
48  const int dim = x.Size();
49  double res = 0.0;
50  for (int d = 0; d < dim; d++) { res += x(d) * x(d); }
51  return res;
52 }
53 
54 void F_exact(const Vector &p, Vector &F)
55 {
56  F(0) = field_func(p);
57  for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*F(0); }
58 }
59 
60 int main (int argc, char *argv[])
61 {
62  // Initialize MPI.
63  int num_procs, myid;
64  MPI_Init(&argc, &argv);
65  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
66  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
67 
68  // Set the method's default parameters.
69  const char *mesh_file = "../../data/rt-2d-q3.mesh";
70  int order = 3;
71  int mesh_poly_deg = 3;
72  int rs_levels = 0;
73  int rp_levels = 0;
74  bool visualization = true;
75  int fieldtype = 0;
76  int ncomp = 1;
77 
78  // Parse command-line options.
79  OptionsParser args(argc, argv);
80  args.AddOption(&mesh_file, "-m", "--mesh",
81  "Mesh file to use.");
82  args.AddOption(&order, "-o", "--order",
83  "Finite element order (polynomial degree).");
84  args.AddOption(&mesh_poly_deg, "-mo", "--mesh-order",
85  "Polynomial degree of mesh finite element space.");
86  args.AddOption(&rs_levels, "-rs", "--refine-serial",
87  "Number of times to refine the mesh uniformly in serial.");
88  args.AddOption(&rp_levels, "-rp", "--refine-parallel",
89  "Number of times to refine the mesh uniformly in parallel.");
90  args.AddOption(&fieldtype, "-ft", "--field-type",
91  "Field type: 0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
92  args.AddOption(&ncomp, "-nc", "--ncomp",
93  "Number of components for H1 or L2 GridFunctions");
94  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
95  "--no-visualization",
96  "Enable or disable GLVis visualization.");
97  args.Parse();
98  if (!args.Good())
99  {
100  args.PrintUsage(cout);
101  return 1;
102  }
103  if (myid == 0) { args.PrintOptions(cout); }
104 
105  // Initialize and refine the starting mesh.
106  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
107  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
108  const int dim = mesh->Dimension();
109  if (myid == 0)
110  {
111  cout << "Mesh curvature of the original mesh: ";
112  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
113  else { cout << "(NONE)"; }
114  cout << endl;
115  }
116 
117  // Mesh bounding box (for the full serial mesh).
118  Vector pos_min, pos_max;
119  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
120  mesh->GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
121  if (myid == 0)
122  {
123  cout << "--- Generating equidistant point for:\n"
124  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
125  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
126  if (dim == 3)
127  {
128  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
129  }
130  }
131 
132  // Distribute the mesh.
133  ParMesh pmesh(MPI_COMM_WORLD, *mesh);
134  delete mesh;
135  for (int lev = 0; lev < rp_levels; lev++) { pmesh.UniformRefinement(); }
136 
137  // Curve the mesh based on the chosen polynomial degree.
138  H1_FECollection fecm(mesh_poly_deg, dim);
139  ParFiniteElementSpace pfespace(&pmesh, &fecm, dim);
140  pmesh.SetNodalFESpace(&pfespace);
141  if (myid == 0)
142  {
143  cout << "Mesh curvature of the curved mesh: " << fecm.Name() << endl;
144  }
145 
146  MFEM_VERIFY(ncomp > 0, "Invalid number of components.");
147  int vec_dim = ncomp;
148  FiniteElementCollection *fec = NULL;
149  if (fieldtype == 0)
150  {
151  fec = new H1_FECollection(order, dim);
152  if (myid == 0) { cout << "H1-GridFunction\n"; }
153  }
154  else if (fieldtype == 1)
155  {
156  fec = new L2_FECollection(order, dim);
157  if (myid == 0) { cout << "L2-GridFunction\n"; }
158  }
159  else if (fieldtype == 2)
160  {
161  fec = new RT_FECollection(order, dim);
162  ncomp = 1;
163  vec_dim = dim;
164  if (myid == 0) { cout << "H(div)-GridFunction\n"; }
165  }
166  else if (fieldtype == 3)
167  {
168  fec = new ND_FECollection(order, dim);
169  ncomp = 1;
170  vec_dim = dim;
171  if (myid == 0) { cout << "H(curl)-GridFunction\n"; }
172  }
173  else
174  {
175  if (myid == 0) { MFEM_ABORT("Invalid FECollection type."); }
176  }
177  ParFiniteElementSpace sc_fes(&pmesh, fec, ncomp);
178  ParGridFunction field_vals(&sc_fes);
179 
180  // Project the GridFunction using VectorFunctionCoefficient.
182  field_vals.ProjectCoefficient(F);
183 
184  // Display the mesh and the field through glvis.
185  if (visualization)
186  {
187  char vishost[] = "localhost";
188  int visport = 19916;
189  socketstream sout;
190  sout.open(vishost, visport);
191  if (!sout)
192  {
193  if (myid == 0)
194  {
195  cout << "Unable to connect to GLVis server at "
196  << vishost << ':' << visport << endl;
197  }
198  }
199  else
200  {
201  sout << "parallel " << num_procs << " " << myid << "\n";
202  sout.precision(8);
203  sout << "solution\n" << pmesh << field_vals;
204  if (dim == 2) { sout << "keys RmjA*****\n"; }
205  if (dim == 3) { sout << "keys mA\n"; }
206  sout << flush;
207  }
208  }
209 
210  // Generate equidistant points in physical coordinates over the whole mesh.
211  // Note that some points might be outside, if the mesh is not a box. Note
212  // also that all tasks search the same points (not mandatory).
213  const int pts_cnt_1D = 10;
214  const int pts_cnt = pow(pts_cnt_1D, dim);
215  Vector vxyz(pts_cnt * dim);
216  if (dim == 2)
217  {
219  const IntegrationRule &ir = el.GetNodes();
220  for (int i = 0; i < ir.GetNPoints(); i++)
221  {
222  const IntegrationPoint &ip = ir.IntPoint(i);
223  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
224  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
225  }
226  }
227  else
228  {
230  const IntegrationRule &ir = el.GetNodes();
231  for (int i = 0; i < ir.GetNPoints(); i++)
232  {
233  const IntegrationPoint &ip = ir.IntPoint(i);
234  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
235  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
236  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
237  }
238  }
239 
240  // Find and Interpolate FE function values on the desired points.
241  Vector interp_vals(pts_cnt*vec_dim);
242  FindPointsGSLIB finder(MPI_COMM_WORLD);
243  finder.Setup(pmesh);
244  finder.Interpolate(vxyz, field_vals, interp_vals);
245  Array<unsigned int> code_out = finder.GetCode();
246  Array<unsigned int> task_id_out = finder.GetProc();
247  Vector dist_p_out = finder.GetDist();
248 
249  int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
250  double max_err = 0.0, max_dist = 0.0;
251  Vector pos(dim);
252  int npt = 0;
253  for (int j = 0; j < vec_dim; j++)
254  {
255  for (int i = 0; i < pts_cnt; i++)
256  {
257  if (j == 0)
258  {
259  (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
260  }
261 
262  if (code_out[i] < 2)
263  {
264  for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
265  Vector exact_val(vec_dim);
266  F_exact(pos, exact_val);
267  max_err = std::max(max_err, fabs(exact_val(j) - interp_vals(npt)));
268  max_dist = std::max(max_dist, dist_p_out(i));
269  if (code_out[i] == 1 && j == 0) { face_pts++; }
270  }
271  else { if (j == 0) { not_found++; } }
272  npt++;
273  }
274  }
275 
276  // Print the results for task 0 since all tasks have the same set of points.
277  if (myid == 0)
278  {
279  cout << setprecision(16)
280  << "Searched unique points: " << pts_cnt
281  << "\nFound on local mesh: " << found_loc
282  << "\nFound on other tasks: " << found_away
283  << "\nMax interp error: " << max_err
284  << "\nMax dist (of found): " << max_dist
285  << "\nPoints not found: " << not_found
286  << "\nPoints on faces: " << face_pts << endl;
287  }
288 
289  // Free the internal gslib data.
290  finder.FreeData();
291 
292  delete fec;
293 
294  MPI_Finalize();
295  return 0;
296 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:575
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const Vector & GetDist() const
Definition: gslib.hpp:171
Arbitrary order L2 elements in 3D on a cube.
Definition: fe.hpp:2499
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:118
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:474
const Array< unsigned int > & GetCode() const
Definition: gslib.hpp:162
int main(int argc, char *argv[])
Definition: ex1.cpp:66
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe.hpp:37
const Array< unsigned int > & GetProc() const
Return MPI rank on which each point was found by FindPoints.
Definition: gslib.hpp:166
double field_func(const Vector &x)
Definition: findpts.cpp:45
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
void F_exact(const Vector &, Vector &)
Definition: ex4.cpp:266
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
virtual const char * Name() const
Definition: fe_coll.hpp:180
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4126
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:272
A general vector function coefficient.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
Arbitrary order L2 elements in 2D on a square.
Definition: fe.hpp:2460
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:333
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:45
Vector data type.
Definition: vector.hpp:51
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:221
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145