MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
field-diff.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
33using namespace mfem;
34using namespace std;
35
36int 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);
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}
@ ClosedUniform
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition fe_base.hpp:35
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points....
Definition gslib.hpp:53
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition gslib.cpp:856
virtual void FreeData()
Definition gslib.cpp:283
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Arbitrary order L2 elements in 3D on a cube.
Definition fe_l2.hpp:79
Arbitrary order L2 elements in 2D on a square.
Definition fe_l2.hpp:45
Vector with associated FE space and LinearFormIntegrators.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Mesh data type.
Definition mesh.hpp:56
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition mesh.cpp:138
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
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:6211
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
int main()
const int visport
const char vishost[]
real_t p(const Vector &x, real_t t)