MFEM  v4.5.1 Finite element discretization library
findpts.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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-p4-tri.mesh -o 4
31 // findpts -m ../../data/inline-tri.mesh -o 3
32 // findpts -m ../../data/inline-quad.mesh -o 3
33 // findpts -m ../../data/inline-quad.mesh -o 3 -po 1
34 // findpts -m ../../data/inline-quad.mesh -o 3 -po 1 -fo 1 -nc 2
35 // findpts -m ../../data/inline-quad.mesh -o 3 -hr -pr
36 // findpts -m ../../data/inline-tet.mesh -o 3
37 // findpts -m ../../data/inline-hex.mesh -o 3
38 // findpts -m ../../data/inline-wedge.mesh -o 3
39 // findpts -m ../../data/amr-quad.mesh -o 2
40 // findpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
41 // findpts -m ../../data/square-mixed.mesh -o 2 -mo 2
42 // findpts -m ../../data/square-mixed.mesh -o 2 -mo 2 -hr -pr
43 // findpts -m ../../data/square-mixed.mesh -o 2 -mo 3 -ft 2
44 // findpts -m ../../data/fichera-mixed.mesh -o 3 -mo 2
45 // findpts -m ../../data/inline-pyramid.mesh -o 1 -mo 1
46 // findpts -m ../../data/tinyzoo-3d.mesh -o 1 -mo 1
47
48 #include "mfem.hpp"
49
50 using namespace mfem;
51 using namespace std;
52
53 // Experimental - required for visualizing functions on p-refined spaces.
54 GridFunction* ProlongToMaxOrder(const GridFunction *x, const int fieldtype)
55 {
56  const FiniteElementSpace *fespace = x->FESpace();
57  Mesh *mesh = fespace->GetMesh();
58  const FiniteElementCollection *fec = fespace->FEColl();
59
60  // find the max order in the space
61  int max_order = 1;
62  for (int i = 0; i < mesh->GetNE(); i++)
63  {
64  max_order = std::max(fespace->GetElementOrder(i), max_order);
65  }
66
67  // create a visualization space of max order for all elements
68  FiniteElementCollection *fecInt = NULL;
69  if (fieldtype == 0)
70  {
71  fecInt = new H1_FECollection(max_order, mesh->Dimension());
72  }
73  else if (fieldtype == 1)
74  {
75  fecInt = new L2_FECollection(max_order, mesh->Dimension());
76  }
77  FiniteElementSpace *spaceInt = new FiniteElementSpace(mesh, fecInt);
78
80  DenseMatrix I;
81
82  GridFunction *xInt = new GridFunction(spaceInt);
83
84  // interpolate solution vector in the larger space
85  for (int i = 0; i < mesh->GetNE(); i++)
86  {
87  Geometry::Type geom = mesh->GetElementGeometry(i);
89
90  Array<int> dofs;
91  fespace->GetElementDofs(i, dofs);
92  Vector elemvect, vectInt;
93  x->GetSubVector(dofs, elemvect);
94
95  const auto *fe = fec->GetFE(geom, fespace->GetElementOrder(i));
96  const auto *feInt = fecInt->GetFE(geom, max_order);
97
98  feInt->GetTransferMatrix(*fe, T, I);
99  spaceInt->GetElementDofs(i, dofs);
100  vectInt.SetSize(dofs.Size());
101
102  I.Mult(elemvect, vectInt);
103  xInt->SetSubVector(dofs, vectInt);
104  }
105
106  xInt->MakeOwner(fecInt);
107  return xInt;
108 }
109
110 double func_order;
111
112 // Scalar function to project
113 double field_func(const Vector &x)
114 {
115  const int dim = x.Size();
116  double res = 0.0;
117  for (int d = 0; d < dim; d++) { res += std::pow(x(d), func_order); }
118  return res;
119 }
120
121 void F_exact(const Vector &p, Vector &F)
122 {
123  F(0) = field_func(p);
124  for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*F(0); }
125 }
126
127 int main (int argc, char *argv[])
128 {
129  // Set the method's default parameters.
130  const char *mesh_file = "../../data/rt-2d-q3.mesh";
131  int order = 3;
132  int mesh_poly_deg = 3;
133  int rs_levels = 0;
134  bool visualization = true;
135  int fieldtype = 0;
136  int ncomp = 1;
137  bool hrefinement = false;
138  bool prefinement = false;
139  int point_ordering = 0;
140  int gf_ordering = 0;
141
142  // Parse command-line options.
143  OptionsParser args(argc, argv);
145  "Mesh file to use.");
147  "Finite element order (polynomial degree).");
149  "Polynomial degree of mesh finite element space.");
151  "Number of times to refine the mesh uniformly in serial.");
153  "Field type: 0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
155  "Number of components for H1 or L2 GridFunctions");
157  "--no-visualization",
158  "Enable or disable GLVis visualization.");
160  "--no-h-refinement",
161  "Do random h refinements to mesh (does not work for pyramids).");
163  "--no-p-refinement",
164  "Do random p refinements to solution field (does not work for pyramids).");
166  "Ordering of points to be found."
167  "0 (default): byNodes, 1: byVDIM");
169  "Ordering of fespace that will be used for gridfunction to be interpolated."
170  "0 (default): byNodes, 1: byVDIM");
171
172  args.Parse();
173  if (!args.Good())
174  {
175  args.PrintUsage(cout);
176  return 1;
177  }
178  args.PrintOptions(cout);
179
180  func_order = std::min(order, 2);
181
182  // Initialize and refine the starting mesh.
183  Mesh mesh(mesh_file, 1, 1, false);
184  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
185  const int dim = mesh.Dimension();
186  cout << "Mesh curvature of the original mesh: ";
187  if (mesh.GetNodes()) { cout << mesh.GetNodes()->OwnFEC()->Name(); }
188  else { cout << "(NONE)"; }
189  cout << endl;
190
191  // Mesh bounding box.
192  Vector pos_min, pos_max;
193  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
194  mesh.GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
195  if (hrefinement || prefinement) { mesh.EnsureNCMesh(); }
196  cout << "--- Generating equidistant point for:\n"
197  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
198  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
199  if (dim == 3)
200  {
201  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
202  }
203
204  // Random h-refinements to mesh
205  if (hrefinement) { mesh.RandomRefinement(0.5); }
206
207  // Curve the mesh based on the chosen polynomial degree.
208  H1_FECollection fecm(mesh_poly_deg, dim);
209  FiniteElementSpace fespace(&mesh, &fecm, dim);
210  mesh.SetNodalFESpace(&fespace);
211  cout << "Mesh curvature of the curved mesh: " << fecm.Name() << endl;
212
213  MFEM_VERIFY(ncomp > 0, "Invalid number of components.");
214  int vec_dim = ncomp;
215  FiniteElementCollection *fec = NULL;
216  if (fieldtype == 0)
217  {
218  fec = new H1_FECollection(order, dim);
219  cout << "H1-GridFunction\n";
220  }
221  else if (fieldtype == 1)
222  {
223  fec = new L2_FECollection(order, dim);
224  cout << "L2-GridFunction\n";
225  }
226  else if (fieldtype == 2)
227  {
228  fec = new RT_FECollection(order, dim);
229  ncomp = 1;
230  vec_dim = dim;
231  cout << "H(div)-GridFunction\n";
232  }
233  else if (fieldtype == 3)
234  {
235  fec = new ND_FECollection(order, dim);
236  ncomp = 1;
237  vec_dim = dim;
238  cout << "H(curl)-GridFunction\n";
239  }
240  else
241  {
242  MFEM_ABORT("Invalid field type.");
243  }
244  FiniteElementSpace sc_fes(&mesh, fec, ncomp, gf_ordering);
245  GridFunction field_vals(&sc_fes);
246
247  // Random p-refinements to the solution field
248  if (prefinement)
249  {
250  for (int e = 0; e < mesh.GetNE(); e++)
251  {
252  if (rand() % 2 == 0)
253  {
254  int element_order = sc_fes.GetElementOrder(e);
255  sc_fes.SetElementOrder(e, element_order + 1);
256  }
257  }
258  sc_fes.Update(false);
259  field_vals.Update();
260  }
261
262  // Project the GridFunction using VectorFunctionCoefficient.
264  field_vals.ProjectCoefficient(F);
265
266  GridFunction *field_vals_pref = prefinement ?
267  ProlongToMaxOrder(&field_vals, fieldtype) :
268  &field_vals;
269
270  // Display the mesh and the field through glvis.
271  if (visualization)
272  {
273  char vishost[] = "localhost";
274  int visport = 19916;
275  socketstream sout;
276  sout.open(vishost, visport);
277  if (!sout)
278  {
279  cout << "Unable to connect to GLVis server at "
280  << vishost << ':' << visport << endl;
281  }
282  else
283  {
284  sout.precision(8);
285  sout << "solution\n" << mesh << *field_vals_pref;
286  if (dim == 2) { sout << "keys RmjA*****\n"; }
287  if (dim == 3) { sout << "keys mA\n"; }
288  sout << flush;
289  }
290  }
291
292  // Generate equidistant points in physical coordinates over the whole mesh.
293  // Note that some points might be outside, if the mesh is not a box. Note
294  // also that all tasks search the same points (not mandatory).
295  const int pts_cnt_1D = 25;
296  int pts_cnt = pow(pts_cnt_1D, dim);
297  Vector vxyz(pts_cnt * dim);
298  if (dim == 2)
299  {
301  const IntegrationRule &ir = el.GetNodes();
302  for (int i = 0; i < ir.GetNPoints(); i++)
303  {
304  const IntegrationPoint &ip = ir.IntPoint(i);
305  if (point_ordering == Ordering::byNODES)
306  {
307  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
308  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
309  }
310  else
311  {
312  vxyz(i*dim + 0) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
313  vxyz(i*dim + 1) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
314  }
315  }
316  }
317  else
318  {
320  const IntegrationRule &ir = el.GetNodes();
321  for (int i = 0; i < ir.GetNPoints(); i++)
322  {
323  const IntegrationPoint &ip = ir.IntPoint(i);
324  if (point_ordering == Ordering::byNODES)
325  {
326  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
327  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
328  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
329  }
330  else
331  {
332  vxyz(i*dim + 0) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
333  vxyz(i*dim + 1) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
334  vxyz(i*dim + 2) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
335  }
336  }
337  }
338
339  // Find and Interpolate FE function values on the desired points.
340  Vector interp_vals(pts_cnt*vec_dim);
341  FindPointsGSLIB finder;
342  finder.Setup(mesh);
344  finder.Interpolate(vxyz, field_vals, interp_vals, point_ordering);
345  Array<unsigned int> code_out = finder.GetCode();
346  Vector dist_p_out = finder.GetDist();
347
348  int face_pts = 0, not_found = 0, found = 0;
349  double err = 0.0, max_err = 0.0, max_dist = 0.0;
350  Vector pos(dim);
351  int npt = 0;
352  for (int j = 0; j < vec_dim; j++)
353  {
354  for (int i = 0; i < pts_cnt; i++)
355  {
356  if (code_out[i] < 2)
357  {
358  if (j == 0) { found++; }
359  for (int d = 0; d < dim; d++)
360  {
361  pos(d) = point_ordering == Ordering::byNODES ?
362  vxyz(d*pts_cnt + i) :
363  vxyz(i*dim + d);
364  }
365  Vector exact_val(vec_dim);
366  F_exact(pos, exact_val);
367  err = gf_ordering == Ordering::byNODES ?
368  fabs(exact_val(j) - interp_vals[i + j*pts_cnt]) :
369  fabs(exact_val(j) - interp_vals[i*vec_dim + j]);
370  max_err = std::max(max_err, err);
371  max_dist = std::max(max_dist, dist_p_out(i));
372  if (code_out[i] == 1 && j == 0) { face_pts++; }
373  }
374  else { if (j == 0) { not_found++; } }
375  npt++;
376  }
377  }
378
379  cout << setprecision(16)
380  << "Searched points: " << pts_cnt
381  << "\nFound points: " << found
382  << "\nMax interp error: " << max_err
383  << "\nMax dist (of found): " << max_dist
385  << "\nPoints on faces: " << face_pts << endl;
386
387  // Free the internal gslib data.
388  finder.FreeData();
389
390  if (prefinement) { delete field_vals_pref; }
391  delete fec;
392
393  return 0;
394 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
const char vishost[]
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:801
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
const FiniteElement * GetFE(Geometry::Type geom, int p) const
Variable order version of FiniteElementForGeometry().
Definition: fe_coll.hpp:165
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:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3334
virtual const Vector & GetDist() const
Definition: gslib.hpp:200
Arbitrary order L2 elements in 3D on a cube.
Definition: fe_l2.hpp:68
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void SetElementOrder(int i, int p)
Sets the order of the i&#39;th finite element.
Definition: fespace.cpp:151
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 GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int GetElementOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:178
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:124
double field_func(const Vector &x)
Definition: findpts.cpp:113
virtual void SetL2AvgType(AvgType avgtype_)
Definition: gslib.hpp:175
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:265
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9473
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
GridFunction * ProlongToMaxOrder(const GridFunction *x, const int fieldtype)
Definition: findpts.cpp:54
void F_exact(const Vector &, Vector &)
Definition: ex4.cpp:268
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
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
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:9497
const int visport
virtual const char * Name() const
Definition: fe_coll.hpp:241
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2678
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
A general vector function coefficient.
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
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
Geometry::Type GetElementGeometry(int i) const
Definition: mesh.hpp:1055
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
Class for integration point with weight.
Definition: intrules.hpp:25
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
int dim
Definition: ex24.cpp:53
virtual const Array< unsigned int > & GetCode() const
Definition: gslib.hpp:191
Arbitrary order L2 elements in 2D on a square.
Definition: fe_l2.hpp:39
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2399
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:415
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:47
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.cpp:118
int main()
double func_order
Definition: findpts.cpp:110
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150