MFEM  v4.5.2
Finite element discretization library
field-diff.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // 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-2.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);
48  args.AddOption(&mesh_file_1, "-m1", "--mesh1",
49  "Mesh file for solution 1.");
50  args.AddOption(&mesh_file_2, "-m2", "--mesh2",
51  "Mesh file for solution 2.");
52  args.AddOption(&sltn_file_1, "-s1", "--solution1",
53  "Grid function for solution 1.");
54  args.AddOption(&sltn_file_2, "-s2", "--solution2",
55  "Grid function for solution 2.");
56  args.AddOption(&pts_cnt_1D, "-p", "--points1D",
57  "Number of comparison points in one direction");
58  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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  MFEM_VERIFY(dim > 1, "GSLIB requires a 2D or a 3D mesh" );
74 
75  // The Nodes GridFunctions for each mesh are required.
76  if (mesh_1.GetNodes() == NULL) { mesh_1.SetCurvature(1, false, dim, 0); }
77  if (mesh_2.GetNodes() == NULL) { mesh_2.SetCurvature(1, false, dim, 0); }
78  const int mesh_poly_deg = mesh_1.GetNodes()->FESpace()->GetElementOrder(0);
79  cout << "Mesh curvature: "
80  << mesh_1.GetNodes()->OwnFEC()->Name() << " " << mesh_poly_deg << endl;
81 
82  // The visualization of the difference assumes byNODES ordering.
83  if (mesh_1.GetNodes()->FESpace()->GetOrdering() == Ordering::byVDIM)
84  { mesh_1.SetCurvature(mesh_poly_deg, false, dim, Ordering::byNODES); }
85  if (mesh_2.GetNodes()->FESpace()->GetOrdering() == Ordering::byVDIM)
86  { mesh_2.SetCurvature(mesh_poly_deg, false, dim, Ordering::byNODES); }
87 
88  // Mesh bounding box.
89  Vector pos_min, pos_max;
90  MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
91  mesh_1.GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
92  cout << "Generating equidistant points for:\n"
93  << " x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
94  << " y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
95  if (dim == 3)
96  {
97  cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
98  }
99 
100  ifstream mat_stream_1(sltn_file_1), mat_stream_2(sltn_file_2);
101  GridFunction func_1(&mesh_1, mat_stream_1);
102  GridFunction func_2(&mesh_2, mat_stream_2);
103 
104  // Display the meshes and the fields through glvis.
105  if (visualization)
106  {
107  char vishost[] = "localhost";
108  int visport = 19916;
109  socketstream sout1, sout2;
110  sout1.open(vishost, visport);
111  sout2.open(vishost, visport);
112  if (!sout1)
113  {
114  cout << "Unable to connect to GLVis server at "
115  << vishost << ':' << visport << endl;
116  }
117  else
118  {
119  sout1.precision(8);
120  sout1 << "solution\n" << mesh_1 << func_1
121  << "window_title 'Solution 1'"
122  << "window_geometry 0 0 600 600";
123  if (dim == 2) { sout1 << "keys RmjAc"; }
124  if (dim == 3) { sout1 << "keys mA\n"; }
125  sout1 << flush;
126 
127  sout2.precision(8);
128  sout2 << "solution\n" << mesh_2 << func_2
129  << "window_title 'Solution 2'"
130  << "window_geometry 600 0 600 600";
131  if (dim == 2) { sout2 << "keys RmjAc"; }
132  if (dim == 3) { sout2 << "keys mA\n"; }
133  sout2 << flush;
134  }
135  }
136 
137  // Generate equidistant points in physical coordinates over the whole mesh.
138  // Note that some points might be outside, if the mesh is not a box. Note
139  // also that all tasks search the same points (not mandatory).
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  FindPointsGSLIB finder1, finder2;
167  Vector interp_vals_1(pts_cnt), interp_vals_2(pts_cnt);
168 
169  // First solution.
170  finder1.Interpolate(mesh_1, vxyz, func_1, interp_vals_1);
171 
172  // Second solution.
173  finder2.Interpolate(mesh_2, vxyz, func_2, interp_vals_2);
174 
175  // Compute differences between the two sets of values.
176  double avg_diff = 0.0, max_diff = 0.0, diff_p;
177  for (int p = 0; p < pts_cnt; p++)
178  {
179  diff_p = fabs(interp_vals_1(p) - interp_vals_2(p));
180  avg_diff += diff_p;
181  if (diff_p > max_diff) { max_diff = diff_p; }
182  }
183  avg_diff /= pts_cnt;
184 
185  GridFunction *n1 = mesh_1.GetNodes(), *n2 = mesh_2.GetNodes();
186  double *nd1 = n1->GetData(), *nd2 = n2->GetData();
187  double avg_dist = 0.0;
188  const int node_cnt = n1->Size() / dim;
189  if (n1->Size() == n2->Size())
190  {
191  for (int i = 0; i < node_cnt; i++)
192  {
193  double diff_i = 0.0;
194  for (int d = 0; d < dim; d++)
195  {
196  const int j = i + d * node_cnt;
197  diff_i += (nd1[j] - nd2[j]) * (nd1[j] - nd2[j]);
198  }
199  avg_dist += sqrt(diff_i);
200  }
201  avg_dist /= node_cnt;
202  }
203  else { avg_dist = -1.0; }
204 
205  std::cout << "Avg position difference: " << avg_dist << std::endl
206  << "Searched " << pts_cnt << " points.\n"
207  << "Max diff: " << max_diff << std::endl
208  << "Avg diff: " << avg_diff << std::endl;
209 
210  // Visualize the difference at the nodes of mesh 1.
211  FiniteElementSpace m1_fes(&mesh_1, mesh_1.GetNodes()->FESpace()->FEColl());
212  GridFunction diff(&m1_fes);
213  GridFunctionCoefficient f1c(&func_1);
214  diff.ProjectDiscCoefficient(f1c, GridFunction::ARITHMETIC);
215  vxyz = *mesh_1.GetNodes();
216  const int nodes_cnt = vxyz.Size() / dim;
217  interp_vals_2.SetSize(nodes_cnt);
218  finder2.Interpolate(vxyz, func_2, interp_vals_2);
219 
220  for (int n = 0; n < nodes_cnt; n++)
221  {
222  diff(n) = fabs(diff(n) - interp_vals_2(n));
223  }
224 
225  if (visualization)
226  {
227  char vishost[] = "localhost";
228  int visport = 19916;
229  socketstream sout3;
230  sout3.open(vishost, visport);
231  sout3.precision(8);
232  sout3 << "solution\n" << mesh_1 << diff
233  << "window_title 'Difference'"
234  << "window_geometry 1200 0 600 600";
235  if (dim == 2) { sout3 << "keys RmjAcpppppppppppppppppppppp"; }
236  if (dim == 3) { sout3 << "keys mA\n"; }
237  sout3 << flush;
238  }
239 
240  ConstantCoefficient coeff1(1.0);
241  DomainLFIntegrator *lf_integ = new DomainLFIntegrator(coeff1);
242  LinearForm lf(func_1.FESpace());
243  lf.AddDomainIntegrator(lf_integ);
244  lf.Assemble();
245  const double vol_diff = diff * lf;
246  std::cout << "Vol diff: " << vol_diff << std::endl;
247 
248  // Free the internal gslib data.
249  finder1.FreeData();
250  finder2.FreeData();
251 
252  return 0;
253 }
const char vishost[]
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:803
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
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:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Arbitrary order L2 elements in 3D on a cube.
Definition: fe_l2.hpp:78
int Dimension() const
Definition: mesh.hpp:1047
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
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:199
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
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:267
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5356
const int visport
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
int main(int argc, char *argv[])
Definition: field-diff.cpp:36
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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_l2.hpp:44
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
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:52
Vector data type.
Definition: vector.hpp:60
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388