MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh-quality.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// Mesh Quality Miniapp: Visualize and Check Mesh Quality
14// ------------------------------------------------------
15//
16// This miniapp extracts geometric parameters from the Jacobian transformation
17// at each degree-of-freedom of the mesh, and visualizes it. The geometric
18// parameters size, skewness, and aspect-ratio can help assess the quality of a
19// mesh.
20//
21// Compile with: make mesh-quality
22//
23// Sample runs:
24// mesh-quality -m blade.mesh -size -aspr -skew -vis -visit
25// mesh-quality -m ../../data/square-disc.mesh -o 2 -size -aspr -skew -vis
26
27#include "mfem.hpp"
28#include <iostream>
29#include "../common/mfem-common.hpp"
30
31using namespace std;
32using namespace mfem;
33
34int main(int argc, char *argv[])
35{
36 // 1. Parse command-line options.
37 const char *mesh_file = "../../data/inline-quad.mesh";
38 int order = -1;
39 int ref_levels = 0;
40 bool visualization = true;
41 bool visit = false;
42 bool size = true;
43 bool aspect_ratio = true;
44 bool skewness = true;
45
46 OptionsParser args(argc, argv);
47 args.AddOption(&mesh_file, "-m", "--mesh",
48 "Mesh file to use.");
49 args.AddOption(&order, "-o", "--order",
50 "Finite element order for visualization"
51 "(defaults to order of the mesh).");
52 args.AddOption(&ref_levels, "-r", "--ref-levels",
53 "Number of initial uniform refinement levels.");
54 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
55 "--no-visualization",
56 "Enable or disable GLVis visualization.");
57 args.AddOption(&visit, "-visit", "--visit", "-no-visit",
58 "--no-visit",
59 "Enable or disable VisIt.");
60 args.AddOption(&size, "-size", "--size", "-no-size",
61 "--no-size",
62 "Visualize size parameter.");
63 args.AddOption(&aspect_ratio, "-aspr", "--aspect-ratio", "-no-aspr",
64 "--no-aspect-ratio",
65 "Visualize aspect-ratio parameter.");
66 args.AddOption(&skewness, "-skew", "--skew", "-no-skew",
67 "--no-skew",
68 "Visualize skewness parameter.");
69 args.Parse();
70 if (!args.Good())
71 {
72 args.PrintUsage(cout);
73 return 1;
74 }
75 args.PrintOptions(cout);
76
77 Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
78 for (int lev = 0; lev < ref_levels; lev++)
79 {
80 mesh->UniformRefinement();
81 }
82 const int dim = mesh->Dimension();
83
84 int nSize = 1, nAspr = 1, nSkew = 1;
85 if (dim == 3)
86 {
87 nAspr = 2;
88 nSkew = 3;
89 }
90
91 // Total number of geometric parameters; for now we skip orientation.
92 const int nTotalParams = nSize + nAspr + nSkew;
93 if (order < 0)
94 {
95 order = mesh->GetNodes() == NULL ? 1 :
96 mesh->GetNodes()->FESpace()->GetFE(0)->GetOrder();
97 }
98
99 // Define a GridFunction for all geometric parameters associated with the
100 // mesh.
101 L2_FECollection l2fec(order, mesh->Dimension());
102 FiniteElementSpace fespace(mesh, &l2fec, nTotalParams); // must order byNodes
103 GridFunction quality(&fespace);
104
105 DenseMatrix jacobian(dim);
106 Vector geomParams(nTotalParams);
107 Array<int> vdofs;
108 Vector allVals;
109 // Compute the geometric parameter at the dofs of each element.
110 for (int e = 0; e < mesh->GetNE(); e++)
111 {
112 const FiniteElement *fe = fespace.GetFE(e);
113 const IntegrationRule &ir = fe->GetNodes();
114 fespace.GetElementVDofs(e, vdofs);
115 allVals.SetSize(vdofs.Size());
116 for (int q = 0; q < ir.GetNPoints(); q++)
117 {
118 const IntegrationPoint &ip = ir.IntPoint(q);
119 mesh->GetElementJacobian(e, jacobian, &ip);
120 real_t sizeVal;
121 Vector asprVals, skewVals, oriVals;
122 mesh->GetGeometricParametersFromJacobian(jacobian, sizeVal,
123 asprVals, skewVals, oriVals);
124 allVals(q + 0) = sizeVal;
125 for (int n = 0; n < nAspr; n++)
126 {
127 allVals(q + (n+1)*ir.GetNPoints()) = asprVals(n);
128 }
129 for (int n = 0; n < nSkew; n++)
130 {
131 allVals(q + (n+1+nAspr)*ir.GetNPoints()) = skewVals(n);
132 }
133 }
134 quality.SetSubVector(vdofs, allVals);
135 }
136
137 VisItDataCollection visit_dc("quality", mesh);
138
139 // Visualize different parameters
140 int visw = 400;
141 int vish = 400;
142 int cx = 0;
143 int cy = 0;
144 int gap = 10;
145 FiniteElementSpace scfespace(mesh, &l2fec);
146 int ndofs = scfespace.GetNDofs();
147 if (dim == 2)
148 {
149 int idx = 0;
150 GridFunction size_gf, aspr_gf, skew_gf;
151 socketstream vis1, vis2, vis3;
152 if (size)
153 {
154 size_gf.MakeRef(&scfespace, quality.GetData());
155 if (visit) { visit_dc.RegisterField("Size", &size_gf); }
156 if (visualization)
157 {
158 common::VisualizeField(vis1, "localhost", 19916, size_gf,
159 "Size", cx, cy, visw, vish, "Rjmc");
160 }
161 real_t min_size = size_gf.Min(),
162 max_size = size_gf.Max();
163 cout << "Min size: " << min_size << endl;
164 cout << "Max size: " << max_size << endl;
165 }
166 idx++;
167 if (aspect_ratio)
168 {
169 aspr_gf.MakeRef(&scfespace, quality.GetData() + idx*ndofs);
170 if (visit) { visit_dc.RegisterField("Aspect-Ratio", &aspr_gf); }
171 if (visualization)
172 {
173 cx += gap+visw;
174 common::VisualizeField(vis2, "localhost", 19916, aspr_gf,
175 "Aspect-Ratio", cx, cy, visw, vish, "Rjmc");
176 }
177 real_t min_aspr = aspr_gf.Min(),
178 max_aspr = aspr_gf.Max();
179 max_aspr = std::max((real_t) 1.0/min_aspr, max_aspr);
180 cout << "Worst aspect-ratio: " << max_aspr << endl;
181 cout << "(in any direction)" << endl;
182 }
183 idx++;
184 if (skewness)
185 {
186 skew_gf.MakeRef(&scfespace, quality.GetData() + idx*ndofs);
187 if (visit) { visit_dc.RegisterField("Skewness", &skew_gf); }
188 if (visualization)
189 {
190 cx += gap+visw;
191 common::VisualizeField(vis3, "localhost", 19916, skew_gf,
192 "Skewness (radians)", cx, cy, visw, vish,
193 "Rjmc");
194 }
195 real_t min_skew = skew_gf.Min(),
196 max_skew = skew_gf.Max();
197 cout << "Min skew (in deg): " << min_skew*180/M_PI << endl;
198 cout << "Max skew (in deg): " << max_skew*180/M_PI << endl;
199 }
200 if (visit)
201 {
203 visit_dc.Save();
204 }
205 }
206 else if (dim == 3)
207 {
208 int idx = 0;
209 GridFunction size_gf, aspr_gf1, aspr_gf2,
210 skew_gf1, skew_gf2, skew_gf3;
211 socketstream vis1, vis2, vis3, vis4, vis5, vis6;
212
213 if (size)
214 {
215 size_gf.MakeRef(&scfespace, quality.GetData());
216 if (visit) { visit_dc.RegisterField("Size", &size_gf); }
217 if (visualization)
218 {
219 common::VisualizeField(vis1, "localhost", 19916, size_gf,
220 "Size", cx, cy, visw, vish, "Rjmc");
221 }
222 real_t min_size = size_gf.Min(),
223 max_size = size_gf.Max();
224 cout << "Min size: " << min_size << endl;
225 cout << "Max size: " << max_size << endl;
226 }
227 idx++;
228
229 if (aspect_ratio)
230 {
231 aspr_gf1.MakeRef(&scfespace, quality.GetData() + (idx++)*ndofs);
232 aspr_gf2.MakeRef(&scfespace, quality.GetData() + (idx++)*ndofs);
233 if (visit)
234 {
235 visit_dc.RegisterField("Aspect-Ratio", &aspr_gf1);
236 visit_dc.RegisterField("Aspect-Ratio2", &aspr_gf2);
237 }
238 if (visualization)
239 {
240 cx += gap+visw;
241 common::VisualizeField(vis2, "localhost", 19916, aspr_gf1,
242 "Aspect-Ratio", cx, cy, visw, vish, "Rjmc");
243 cx += gap+visw;
244 common::VisualizeField(vis3, "localhost", 19916, aspr_gf2,
245 "Aspect-Ratio2", cx, cy, visw, vish, "Rjmc");
246 }
247 real_t min_aspr1 = aspr_gf1.Min(),
248 max_aspr1 = aspr_gf1.Max();
249 max_aspr1 = std::max((real_t) 1.0/min_aspr1, max_aspr1);
250
251 real_t min_aspr2 = aspr_gf2.Min(),
252 max_aspr2 = aspr_gf2.Max();
253 max_aspr2 = std::max((real_t) 1.0/min_aspr2, max_aspr2);
254 real_t max_aspr = max(max_aspr1, max_aspr2);
255
256 Vector aspr_gf3(aspr_gf1.Size());
257 for (int i = 0; i < aspr_gf1.Size(); i++)
258 {
259 aspr_gf3(i) = 1.0/(aspr_gf1(i)*aspr_gf2(i));
260 }
261 real_t min_aspr3 = aspr_gf3.Min(),
262 max_aspr3 = aspr_gf3.Max();
263 max_aspr3 = std::max((real_t) 1.0/min_aspr3, max_aspr3);
264 max_aspr = std::max(max_aspr, max_aspr3);
265
266 cout << "Worst aspect-ratio: " << max_aspr << endl;
267 cout << "(in any direction)" << endl;
268 }
269 else { idx += 2; }
270
271 if (skewness)
272 {
273 skew_gf1.MakeRef(&scfespace, quality.GetData() + (idx++)*ndofs);
274 skew_gf2.MakeRef(&scfespace, quality.GetData() + (idx++)*ndofs);
275 skew_gf3.MakeRef(&scfespace, quality.GetData() + (idx++)*ndofs);
276 if (visit)
277 {
278 visit_dc.RegisterField("Skewness (radians)", &skew_gf1);
279 visit_dc.RegisterField("Skewness2 (radians)", &skew_gf2);
280 visit_dc.RegisterField("Dihedral (radians)", &skew_gf3);
281 }
282 if (visualization)
283 {
284 cx = 0;
285 cy += 10*gap+vish;
286 common::VisualizeField(vis4, "localhost", 19916, skew_gf1,
287 "Skewness", cx, cy, visw, vish, "Rjmc");
288 cx += gap+visw;
289 common::VisualizeField(vis5, "localhost", 19916, skew_gf2,
290 "Skewness2", cx, cy, visw, vish, "Rjmc");
291 cx += gap+visw;
292 common::VisualizeField(vis6, "localhost", 19916, skew_gf3,
293 "Dihedral", cx, cy, visw, vish, "Rjmc");
294 }
295 real_t min_skew1 = skew_gf1.Min(),
296 max_skew1 = skew_gf1.Max();
297 real_t min_skew2 = skew_gf2.Min(),
298 max_skew2 = skew_gf2.Max();
299 real_t min_skew3 = skew_gf3.Min(),
300 max_skew3 = skew_gf3.Max();
301 cout << "Min skew 1 (in deg): " << min_skew1*180/M_PI << endl;
302 cout << "Max skew 1 (in deg): " << max_skew1*180/M_PI << endl;
303
304 cout << "Min skew 2 (in deg): " << min_skew2*180/M_PI << endl;
305 cout << "Max skew 2 (in deg): " << max_skew2*180/M_PI << endl;
306
307 cout << "Min skew 3 (in deg): " << min_skew3*180/M_PI << endl;
308 cout << "Max skew 3 (in deg): " << max_skew3*180/M_PI << endl;
309 }
310 else { idx += 3; }
311
312 if (visit)
313 {
315 visit_dc.Save();
316 }
317 }
318
319 delete mesh;
320 return 0;
321}
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
virtual void SetFormat(int fmt)
Set the desired output mesh and data format.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:710
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition: fespace.cpp:287
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition: fespace.cpp:3168
Abstract class for all finite elements.
Definition: fe_base.hpp:239
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:395
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:203
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-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Mesh data type.
Definition: mesh.hpp:56
void GetElementJacobian(int i, DenseMatrix &J, const IntegrationPoint *ip=NULL)
Definition: mesh.cpp:62
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
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
void GetGeometricParametersFromJacobian(const DenseMatrix &J, real_t &volume, Vector &aspr, Vector &skew, Vector &ori) const
Computes geometric parameters associated with a Jacobian matrix in 2D/3D. These parameters are (1) Ar...
Definition: mesh.cpp:13195
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10970
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
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.
Definition: optparser.hpp:159
Vector data type.
Definition: vector.hpp:80
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:604
real_t Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:922
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
real_t Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1212
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition: ex24.cpp:53
int main()
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:92
float real_t
Definition: config.hpp:43