36int main (
int argc,
char *argv[])
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;
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",
60 "Enable or disable GLVis visualization.");
70 Mesh mesh_1(mesh_file_1, 1, 1,
false);
71 Mesh mesh_2(mesh_file_2, 1, 1,
false);
73 MFEM_VERIFY(
dim > 1,
"GSLIB requires a 2D or a 3D mesh" );
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;
90 MFEM_VERIFY(mesh_poly_deg > 0,
"The order of the mesh must be positive.");
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";
97 cout <<
"z in [" << pos_min(2) <<
", " << pos_max(2) <<
"]\n";
100 ifstream mat_stream_1(sltn_file_1), mat_stream_2(sltn_file_2);
114 cout <<
"Unable to connect to GLVis server at "
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"; }
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"; }
140 const int pts_cnt = pow(pts_cnt_1D,
dim);
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));
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));
167 Vector interp_vals_1(pts_cnt), interp_vals_2(pts_cnt);
170 finder1.
Interpolate(mesh_1, vxyz, func_1, interp_vals_1);
173 finder2.
Interpolate(mesh_2, vxyz, func_2, interp_vals_2);
176 double avg_diff = 0.0, max_diff = 0.0, diff_p;
177 for (
int p = 0;
p < pts_cnt;
p++)
179 diff_p = fabs(interp_vals_1(
p) - interp_vals_2(
p));
181 if (diff_p > max_diff) { max_diff = diff_p; }
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())
191 for (
int i = 0; i < node_cnt; i++)
194 for (
int d = 0; d <
dim; d++)
196 const int j = i + d * node_cnt;
197 diff_i += (nd1[j] - nd2[j]) * (nd1[j] - nd2[j]);
199 avg_dist += sqrt(diff_i);
201 avg_dist /= node_cnt;
203 else { avg_dist = -1.0; }
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;
216 const int nodes_cnt = vxyz.
Size() /
dim;
217 interp_vals_2.
SetSize(nodes_cnt);
220 for (
int n = 0; n < nodes_cnt; n++)
222 diff(n) = fabs(diff(n) - interp_vals_2(n));
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"; }
245 const double vol_diff = diff * lf;
246 std::cout <<
"Vol diff: " << vol_diff << std::endl;
@ ClosedUniform
Nodes: x_i = i/(n-1), i=0,...,n-1.
A coefficient that is constant across space and time.
Class for domain integration .
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points....
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
FiniteElementSpace * FESpace()
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Arbitrary order L2 elements in 3D on a cube.
Arbitrary order L2 elements in 2D on a square.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
int Dimension() const
Dimension of the reference space used within the elements.
void GetNodes(Vector &node_coord) const
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.
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...
bool Good() const
Return true if the command line options were parsed successfully.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
real_t p(const Vector &x, real_t t)