MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
findpts.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 // 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-tet.mesh -o 3
34 // findpts -m ../../data/inline-hex.mesh -o 3
35 // findpts -m ../../data/inline-wedge.mesh -o 3
36 // findpts -m ../../data/amr-quad.mesh -o 2
37 // findpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
38 
39 #include "mfem.hpp"
40 
41 using namespace mfem;
42 using namespace std;
43 
44 // Scalar function to project
45 double field_func(const Vector &x)
46 {
47  const int dim = x.Size();
48  double res = 0.0;
49  for (int d = 0; d < dim; d++) { res += x(d) * x(d); }
50  return res;
51 }
52 
53 void F_exact(const Vector &p, Vector &F)
54 {
55  F(0) = field_func(p);
56  for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*F(0); }
57 }
58 
59 int main (int argc, char *argv[])
60 {
61  // Set the method's default parameters.
62  const char *mesh_file = "../../data/rt-2d-q3.mesh";
63  int order = 3;
64  int mesh_poly_deg = 3;
65  int rs_levels = 0;
66  bool visualization = true;
67  int fieldtype = 0;
68  int ncomp = 1;
69 
70  // Parse command-line options.
71  OptionsParser args(argc, argv);
72  args.AddOption(&mesh_file, "-m", "--mesh",
73  "Mesh file to use.");
74  args.AddOption(&order, "-o", "--order",
75  "Finite element order (polynomial degree).");
76  args.AddOption(&mesh_poly_deg, "-mo", "--mesh-order",
77  "Polynomial degree of mesh finite element space.");
78  args.AddOption(&rs_levels, "-rs", "--refine-serial",
79  "Number of times to refine the mesh uniformly in serial.");
80  args.AddOption(&fieldtype, "-ft", "--field-type",
81  "Field type: 0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
82  args.AddOption(&ncomp, "-nc", "--ncomp",
83  "Number of components for H1 or L2 GridFunctions");
84  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
85  "--no-visualization",
86  "Enable or disable GLVis visualization.");
87  args.Parse();
88  if (!args.Good())
89  {
90  args.PrintUsage(cout);
91  return 1;
92  }
93  args.PrintOptions(cout);
94 
95  // Initialize and refine the starting mesh.
96  Mesh mesh(mesh_file, 1, 1, false);
97  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
98  const int dim = mesh.Dimension();
99  cout << "Mesh curvature of the original mesh: ";
100  if (mesh.GetNodes()) { cout << mesh.GetNodes()->OwnFEC()->Name(); }
101  else { cout << "(NONE)"; }
102  cout << endl;
103 
104  // Mesh bounding box.
105  Vector pos_min, pos_max;
106  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
107  mesh.GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
108  cout << "--- Generating equidistant point for:\n"
109  << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
110  << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
111  if (dim == 3)
112  {
113  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
114  }
115 
116  // Curve the mesh based on the chosen polynomial degree.
117  H1_FECollection fecm(mesh_poly_deg, dim);
118  FiniteElementSpace fespace(&mesh, &fecm, dim);
119  mesh.SetNodalFESpace(&fespace);
120  cout << "Mesh curvature of the curved mesh: " << fecm.Name() << endl;
121 
122  MFEM_VERIFY(ncomp > 0, "Invalid number of components.");
123  int vec_dim = ncomp;
124  FiniteElementCollection *fec = NULL;
125  if (fieldtype == 0)
126  {
127  fec = new H1_FECollection(order, dim);
128  cout << "H1-GridFunction\n";
129  }
130  else if (fieldtype == 1)
131  {
132  fec = new L2_FECollection(order, dim);
133  cout << "L2-GridFunction\n";
134  }
135  else if (fieldtype == 2)
136  {
137  fec = new RT_FECollection(order, dim);
138  ncomp = 1;
139  vec_dim = dim;
140  cout << "H(div)-GridFunction\n";
141  }
142  else if (fieldtype == 3)
143  {
144  fec = new ND_FECollection(order, dim);
145  ncomp = 1;
146  vec_dim = dim;
147  cout << "H(curl)-GridFunction\n";
148  }
149  else
150  {
151  MFEM_ABORT("Invalid field type.");
152  }
153  FiniteElementSpace sc_fes(&mesh, fec, ncomp);
154  GridFunction field_vals(&sc_fes);
155 
156  // Project the GridFunction using VectorFunctionCoefficient.
158  field_vals.ProjectCoefficient(F);
159 
160  // Display the mesh and the field through glvis.
161  if (visualization)
162  {
163  char vishost[] = "localhost";
164  int visport = 19916;
165  socketstream sout;
166  sout.open(vishost, visport);
167  if (!sout)
168  {
169  cout << "Unable to connect to GLVis server at "
170  << vishost << ':' << visport << endl;
171  }
172  else
173  {
174  sout.precision(8);
175  sout << "solution\n" << mesh << field_vals;
176  if (dim == 2) { sout << "keys RmjA*****\n"; }
177  if (dim == 3) { sout << "keys mA\n"; }
178  sout << flush;
179  }
180  }
181 
182  // Generate equidistant points in physical coordinates over the whole mesh.
183  // Note that some points might be outside, if the mesh is not a box. Note
184  // also that all tasks search the same points (not mandatory).
185  const int pts_cnt_1D = 25;
186  int pts_cnt = pow(pts_cnt_1D, dim);
187  Vector vxyz(pts_cnt * dim);
188  if (dim == 2)
189  {
191  const IntegrationRule &ir = el.GetNodes();
192  for (int i = 0; i < ir.GetNPoints(); i++)
193  {
194  const IntegrationPoint &ip = ir.IntPoint(i);
195  vxyz(i) = 100*pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
196  vxyz(pts_cnt + i) = 100*pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
197  }
198  }
199  else
200  {
202  const IntegrationRule &ir = el.GetNodes();
203  for (int i = 0; i < ir.GetNPoints(); i++)
204  {
205  const IntegrationPoint &ip = ir.IntPoint(i);
206  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
207  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
208  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
209  }
210  }
211 
212  // Find and Interpolate FE function values on the desired points.
213  Vector interp_vals(pts_cnt*vec_dim);
214  FindPointsGSLIB finder;
215  finder.Setup(mesh);
217  finder.Interpolate(vxyz, field_vals, interp_vals);
218  Array<unsigned int> code_out = finder.GetCode();
219  Vector dist_p_out = finder.GetDist();
220 
221  int face_pts = 0, not_found = 0, found = 0;
222  double max_err = 0.0, max_dist = 0.0;
223  Vector pos(dim);
224  int npt = 0;
225  for (int j = 0; j < vec_dim; j++)
226  {
227  for (int i = 0; i < pts_cnt; i++)
228  {
229  if (code_out[i] < 2)
230  {
231  if (j == 0) { found++; }
232  for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
233  Vector exact_val(vec_dim);
234  F_exact(pos, exact_val);
235  max_err = std::max(max_err, fabs(exact_val(j) - interp_vals[npt]));
236  max_dist = std::max(max_dist, dist_p_out(i));
237  if (code_out[i] == 1 && j == 0) { face_pts++; }
238  }
239  else { if (j == 0) { not_found++; } }
240  npt++;
241  }
242  }
243 
244  cout << setprecision(16)
245  << "Searched points: " << pts_cnt
246  << "\nFound points: " << found
247  << "\nMax interp error: " << max_err
248  << "\nMax dist (of found): " << max_dist
249  << "\nPoints not found: " << not_found
250  << "\nPoints on faces: " << face_pts << endl;
251 
252  // Free the internal gslib data.
253  finder.FreeData();
254 
255  delete fec;
256 
257  return 0;
258 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:582
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 const Vector & GetDist() const
Definition: gslib.hpp:178
Arbitrary order L2 elements in 3D on a cube.
Definition: fe.hpp:2643
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe.hpp:37
double field_func(const Vector &x)
Definition: findpts.cpp:45
virtual void SetL2AvgType(AvgType avgtype_)
Definition: gslib.hpp:153
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:150
constexpr char vishost[]
virtual void FreeData()
Definition: gslib.cpp:212
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
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:75
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.hpp:390
virtual const char * Name() const
Definition: fe_coll.hpp:238
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4843
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
A general vector function coefficient.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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
virtual const Array< unsigned int > & GetCode() const
Definition: gslib.hpp:169
Arbitrary order L2 elements in 2D on a square.
Definition: fe.hpp:2604
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2278
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:416
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:47
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
int main()
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150