MFEM  v4.2.0 Finite element discretization library
field-diff.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 // Field Diff Miniapp: Compare grid functions on different meshes
14 // --------------------------------------------------------------
15 //
16 // This miniapp compares two different high-order grid functions, defined on two
17 // different high-order meshes, based on the GSLIB-FindPoints general off-grid
18 // interpolation utility. Using a set of points defined within the bounding box
19 // of the domain, FindPoints is used to interpolate the grid functions from the
20 // two different meshes and output the difference between the interpolated
21 // values. The miniapp also uses FindPoints to interpolate the solution from one
22 // mesh onto another, and visualize the difference using GLVis.
23 //
24 // Compile with: make field-diff
25 //
26 // Sample runs:
27 // field-diff
28 // field-diff -m1 triple-pt-1.mesh -s1 triple-pt-1.gf -m2 triple-pt-2.mesh -s2 triple-pt-1.gf -p 200
29
30 #include "mfem.hpp"
31 #include <fstream>
32
33 using namespace mfem;
34 using namespace std;
35
36 int main (int argc, char *argv[])
37 {
38  // Set the method's default parameters.
39  const char *mesh_file_1 = "triple-pt-1.mesh";
40  const char *mesh_file_2 = "triple-pt-2.mesh";
41  const char *sltn_file_1 = "triple-pt-1.gf";
42  const char *sltn_file_2 = "triple-pt-2.gf";
43  bool visualization = true;
44  int pts_cnt_1D = 100;
45
46  // Parse command-line options.
47  OptionsParser args(argc, argv);
49  "Mesh file for solution 1.");
51  "Mesh file for solution 2.");
53  "Grid function for solution 1.");
55  "Grid function for solution 2.");
57  "Number of comparison points in one direction");
59  "--no-visualization",
60  "Enable or disable GLVis visualization.");
61  args.Parse();
62  if (!args.Good())
63  {
64  args.PrintUsage(cout);
65  return 1;
66  }
67  args.PrintOptions(cout);
68
69  // Input meshes.
70  Mesh mesh_1(mesh_file_1, 1, 1, false);
71  Mesh mesh_2(mesh_file_2, 1, 1, false);
72  const int dim = mesh_1.Dimension();
73
74  MFEM_VERIFY(dim > 1, "GSLIB requires a 2D or a 3D mesh" );
75  if (mesh_1.GetNodes() == NULL) { mesh_1.SetCurvature(1); }
76  if (mesh_2.GetNodes() == NULL) { mesh_2.SetCurvature(1); }
77  const int mesh_poly_deg = mesh_1.GetNodes()->FESpace()->GetOrder(0);
78  cout << "Mesh curvature: "
79  << mesh_1.GetNodes()->OwnFEC()->Name() << " " << mesh_poly_deg << endl;
80
81  // Mesh bounding box.
82  Vector pos_min, pos_max;
83  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
84  mesh_1.GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
85  cout << "Generating equidistant points for:\n"
86  << " x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
87  << " y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
88  if (dim == 3)
89  {
90  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
91  }
92
93  ifstream mat_stream_1(sltn_file_1), mat_stream_2(sltn_file_2);
94  GridFunction func_1(&mesh_1, mat_stream_1);
95  GridFunction func_2(&mesh_2, mat_stream_2);
96
97  // Display the meshes and the fields through glvis.
98  if (visualization)
99  {
100  char vishost[] = "localhost";
101  int visport = 19916;
102  socketstream sout1, sout2;
103  sout1.open(vishost, visport);
104  sout2.open(vishost, visport);
105  if (!sout1)
106  {
107  cout << "Unable to connect to GLVis server at "
108  << vishost << ':' << visport << endl;
109  }
110  else
111  {
112  sout1.precision(8);
113  sout1 << "solution\n" << mesh_1 << func_1
114  << "window_title 'Solution 1'"
115  << "window_geometry 0 0 600 600";
116  if (dim == 2) { sout1 << "keys RmjAc"; }
117  if (dim == 3) { sout1 << "keys mA\n"; }
118  sout1 << flush;
119
120  sout2.precision(8);
121  sout2 << "solution\n" << mesh_2 << func_2
122  << "window_title 'Solution 2'"
123  << "window_geometry 600 0 600 600";
124  if (dim == 2) { sout2 << "keys RmjAc"; }
125  if (dim == 3) { sout2 << "keys mA\n"; }
126  sout2 << flush;
127  }
128  }
129
130  // Generate equidistant points in physical coordinates over the whole mesh.
131  // Note that some points might be outside, if the mesh is not a box. Note
132  // also that all tasks search the same points (not mandatory).
133  const int pts_cnt = pow(pts_cnt_1D, dim);
134  Vector vxyz(pts_cnt * dim);
135  if (dim == 2)
136  {
138  const IntegrationRule &ir = el.GetNodes();
139  for (int i = 0; i < ir.GetNPoints(); i++)
140  {
141  const IntegrationPoint &ip = ir.IntPoint(i);
142  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
143  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
144  }
145  }
146  else
147  {
149  const IntegrationRule &ir = el.GetNodes();
150  for (int i = 0; i < ir.GetNPoints(); i++)
151  {
152  const IntegrationPoint &ip = ir.IntPoint(i);
153  vxyz(i) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
154  vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
155  vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
156  }
157  }
158
159  FindPointsGSLIB finder1, finder2;
160  Vector interp_vals_1(pts_cnt), interp_vals_2(pts_cnt);
161
162  // First solution.
163  finder1.Interpolate(mesh_1, vxyz, func_1, interp_vals_1);
164
165  // Second solution.
166  finder2.Interpolate(mesh_2, vxyz, func_2, interp_vals_2);
167
168  // Compute differences between the two sets of values.
169  double avg_diff = 0.0, max_diff = 0.0, diff_p;
170  for (int p = 0; p < pts_cnt; p++)
171  {
172  diff_p = fabs(interp_vals_1(p) - interp_vals_2(p));
173  avg_diff += diff_p;
174  if (diff_p > max_diff) { max_diff = diff_p; }
175  }
176  avg_diff /= pts_cnt;
177
178  GridFunction *n1 = mesh_1.GetNodes(), *n2 = mesh_2.GetNodes();
179  double *nd1 = n1->GetData(), *nd2 = n2->GetData();
180  double avg_dist = 0.0;
181  const int node_cnt = n1->Size() / dim;
182  if (n1->Size() == n2->Size())
183  {
184  for (int i = 0; i < node_cnt; i++)
185  {
186  double diff_i = 0.0;
187  for (int d = 0; d < dim; d++)
188  {
189  const int j = i + d * node_cnt;
190  diff_i += (nd1[j] - nd2[j]) * (nd1[j] - nd2[j]);
191  }
192  avg_dist += sqrt(diff_i);
193  }
194  avg_dist /= node_cnt;
195  }
196  else { avg_dist = -1.0; }
197
198  std::cout << "Avg position difference: " << avg_dist << std::endl
199  << "Searched " << pts_cnt << " points.\n"
200  << "Max diff: " << max_diff << std::endl
201  << "Avg diff: " << avg_diff << std::endl;
202
203  // This is used only for visualization and approximating the volume of the
204  // differences.
205  GridFunction diff(func_1.FESpace());
206  vxyz = *mesh_1.GetNodes();
207  const int nodes_cnt = vxyz.Size() / dim;
208
209  // Difference at the nodes of mesh 1.
210  interp_vals_2.SetSize(nodes_cnt);
211  finder2.Interpolate(vxyz, func_2, interp_vals_2);
212  for (int n = 0; n < nodes_cnt; n++)
213  {
214  diff(n) = fabs(func_1(n) - interp_vals_2(n));
215  }
216
217  if (visualization)
218  {
219  char vishost[] = "localhost";
220  int visport = 19916;
221  socketstream sout3;
222  sout3.open(vishost, visport);
223  sout3.precision(8);
224  sout3 << "solution\n" << mesh_1 << diff
225  << "window_title 'Difference'"
226  << "window_geometry 1200 0 600 600";
227  if (dim == 2) { sout3 << "keys RmjAcpppppppppppppppppppppp"; }
228  if (dim == 3) { sout3 << "keys mA\n"; }
229  sout3 << flush;
230  }
231
232  ConstantCoefficient coeff1(1.0);
233  DomainLFIntegrator *lf_integ = new DomainLFIntegrator(coeff1);
234  LinearForm lf(func_1.FESpace());
236  lf.Assemble();
237  const double vol_diff = diff * lf;
238  std::cout << "Vol diff: " << vol_diff << std::endl;
239
240  // Free the internal gslib data.
241  finder1.FreeData();
242  finder2.FreeData();
243
244  return 0;
245 }
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 domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
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
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
Arbitrary order L2 elements in 3D on a cube.
Definition: fe.hpp:2499
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
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
int main(int argc, char *argv[])
Definition: ex1.cpp:66
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe.hpp:37
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[]
constexpr int visport
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4165
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.hpp:382
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:595
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
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;.
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
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145