1 // Gmsh - Copyright (C) 1997-2020 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Author: Maxence Reberol
7 
8 #pragma once
9 
10 #include <float.h>
11 #include "qmtMeshUtils.h"
12 
13 class SurfaceProjector;
14 
15 /**
16  * @brief Mesh optimization statistics
17  */
18 struct GeomOptimStats {
19   size_t nFree = 0;
20   size_t nLock = 0;
21   double timeCpu = 0.;
22   size_t outerLoopIter = 0;
23   double sicnMinBefore = 0.;
24   double sicnMinAfter = 0.;
25   double sicnAvgBefore = 0.;
26   double sicnAvgAfter = 0.;
27   void print();
28 };
29 
30 /**
31  * @brief Optimize the vertex position inside a patch by
32  *        smoothing the CAD parametric coordinates with
33  *        a Laplacian solver: Du=0, Dv=0
34  *        The patch boundary is fixed.
35  *        Use direct solver and arithmetric average to ensure
36  *        maximum principle.
37  *        Not always applicable: if no param, if the patch contains a CAD uv singularity
38  *
39  * @param[in,out] patch The mesh patch to smooth. The new positions and uv are directly updated
40  *                      in the MVertex instances.
41  * @param[out] stats Some statistics on the smoothing
42  *
43  * @return 0 if success
44  */
45 int patchOptimizeGeometryGlobal(
46     GFaceMeshPatch& patch,
47     GeomOptimStats& stats);
48 
49 /**
50  * @brief Options for the kernel-based explicit untangling/smoothing.
51  * Kernels:
52  * - DMO (brute-force uv sampling) if low quality
53  * - Winslow if regular vertex
54  * - Angle-based smoothing is irregular vertex
55  */
56 enum SmoothingKernel {
57   Laplacian,
58   WinslowFDM,
59   AngleBased
60 };
61 
62 struct GeomOptimOptions {
63   double useDmoIfSICNbelow = 0.1; /* use the DMO kernel if SICN quality below threshold */
64   size_t outerLoopIterMax = 100; /* maximum number of loops over all vertices */
65   double timeMax = DBL_MAX; /* stop smoothing is timeMax elapsed */
66   bool invertCADNormals = false; /* invert the CAD normals for the quality computation */
67   SurfaceProjector* sp = nullptr; /* if present, surface projection is used instead of CAD */
68   double smartMinThreshold = -DBL_MAX; /* do not displace if inducing min(SICN) < smartMinThreshold */
69   bool qualityRangeTechnique = false;
70   bool localLocking = false; /* Lock if small displacement, unlocked neighbors else */
71   double dxLocalMax = 1.e-5; /* lock a vertex if dx < dxLocalMax * local_size */
72   double dxGlobalMax = 1.e-5; /* stop if sum(dx) < dxGlobalMax * sum(dx_0) */
73   double qualityRangeMin = 0.5;
74   double qualityRangeMax = 0.8;
75   bool withBackup = true; /* save the geometry before, restore if quality decreased */
76   bool project = true; /* project with SurfaceProjector/CAD */
77   bool force3DwithProjection = false;
78   SmoothingKernel kernelRegular = SmoothingKernel::WinslowFDM;
79   SmoothingKernel kernelIrregular = SmoothingKernel::Laplacian;
80 };
81 
82 /**
83  * @brief Optimize the mesh by iteratively moving the vertices (explicit approach).
84  *        The patch boundary is fixed.
85  *        Require a parametrization on the face.
86  *
87  * @param[in,out] patch The mesh patch to smooth. The new positions and uv are directly updated
88  *                      in the MVertex instances.
89  * @param[in] opt The optimization parameters
90  * @param[out] stats Some statistics on the smoothing
91  *
92  * @return true if success
93  */
94 bool patchOptimizeGeometryWithKernel(
95     GFaceMeshPatch& patch,
96     const GeomOptimOptions& opt,
97     GeomOptimStats& stats);
98 
99 /**
100  * @brief Compute minimum and average SICN quality of elements
101  *
102  * @param[in] elements The elements on which to measure the quality
103  * @param[out] sicnMin Minimum element SICN quality
104  * @param[out] sicnAvg Average element SICN quality
105  */
106 void computeSICN(const std::vector<MElement*>& elements, double& sicnMin, double& sicnAvg);
107 
108 
109 /**
110  * @brief Project the patch interior vertices on the surface.
111  *        If no parametrization, only the SurfaceProjector is used.
112  *        If param, the SurfaceProjector is used as a initial guess.
113  *
114  * @param patch The patch containing the vertices
115  * @param sp If not nullptr, use it to project
116  *
117  * @return true if success
118  */
119 bool patchProjectOnSurface(GFaceMeshPatch& patch, SurfaceProjector* sp = nullptr);
120 
121 /**
122  * @brief High-level function which try to make good parameter choices
123  *        automatically.
124  *
125  * @param gf The face containing the quad mesh to smooth
126  * @param sp Surface projector (faster than CAD projection)
127  * @param timeMax Time budget for the smoothing
128  * @param withBackup Rollback if quality decreased
129  *
130  * @return true if success
131  */
132 bool optimizeGeometryQuadMesh(GFace* gf, SurfaceProjector* sp = nullptr, double timeMax = DBL_MAX, bool withBackup = true);
133 
134 /**
135  * @brief High-level function which try to make good parameter choices
136  *        automatically.
137  *
138  * @param gf The face containing the quad-tri mesh to smooth
139  * @param sp Surface projector (faster than CAD projection)
140  * @param timeMax Time budget for the smoothing
141  *
142  * @return true if success
143  */
144 bool optimizeGeometryQuadTriMesh(GFace* gf, SurfaceProjector* sp = nullptr, double timeMax = DBL_MAX);
145 
146 
147 class GeometryOptimizer {
148   public:
149     enum PlanarMethod {
150       MeanPlane,
151       ParamCAD
152     };
153 
GeometryOptimizer()154     GeometryOptimizer() {};
155     bool initialize(GFaceMeshPatch& patch, SurfaceProjector* sp = nullptr);
156 
157     bool smoothWithKernel(
158         SmoothingKernel kernelRegular,
159         SmoothingKernel kernelIrregular,
160         double timeMax = DBL_MAX,
161         int iterMax = 1000,
162         double withBackup = false,
163         bool localLocking = false, /* Lock if small displacement, unlocked neighbors else */
164         double dxLocalMax = 1.e-5, /* lock a vertex if dx < dxLocalMax * local_size */
165         double dxGlobalMax = 1.e-5, /* stop if sum(dx) < dxGlobalMax * sum(dx_0) */
166         bool project = true, /* project with SurfaceProjector */
167         bool finalCADprojection = true, /* project with real CAD at the end */
168         bool smartVariant = false  /* compute SICN before/after each point displacement, reject quality decrease */
169         );
170 
171     bool smoothWithWinslowUntangler(
172         PlanarMethod planar,
173         int iterMax = 1000,
174         double withBackup = false,
175         bool finalCADprojection = true,
176         double nonPlanarRatioMax = 0.1 /* if using mean plane, check the orthogonal deviation is less than the ratio */
177         );
178 
179   public:
180     GFaceMeshPatch* patchPtr = nullptr;
181     SurfaceProjector* sp = nullptr;
182   protected:
183     std::vector<std::array<double,3> > points;
184     std::vector<std::array<double,2> > uvs;
185     std::vector<std::array<uint32_t,4> > quads;
186     std::vector<size_t> one_ring_first;
187     std::vector<uint32_t> one_ring_values;
188     std::vector<MVertex*> new2old;
189     std::unordered_map<MVertex*,uint32_t> old2new;
190     size_t nFree =0;
191 };
192 
193 bool optimizeGeometryQuadqs(GModel* gm);
194