1 // MeshOptimizer - Copyright (C) 2013-2019 UCLouvain-ULiege
2 //
3 // Permission is hereby granted, free of charge, to any person
4 // obtaining a copy of this software and associated documentation
5 // files (the "Software"), to deal in the Software without
6 // restriction, including without limitation the rights to use, copy,
7 // modify, merge, publish, distribute, and/or sell copies of the
8 // Software, and to permit persons to whom the Software is furnished
9 // to do so, provided that the above copyright notice(s) and this
10 // permission notice appear in all copies of the Software and that
11 // both the above copyright notice(s) and this permission notice
12 // appear in supporting documentation.
13 //
14 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 // NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
18 // COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
19 // ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
20 // DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
21 // WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
22 // ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
23 // OF THIS SOFTWARE.
24 //
25 // Contributors: Thomas Toulorge, Jonathan Lambrechts
26
27 #ifndef PATCH_H
28 #define PATCH_H
29
30 #include <vector>
31 #include <map>
32 #include <set>
33 #include "fullMatrix.h"
34 #include "SPoint3.h"
35 #include "VertexCoord.h"
36
37 class GEntity;
38 class MVertex;
39 class MElement;
40
41 class Patch {
42 public:
43 Patch(const std::map<MElement *, GEntity *> &element2entity,
44 const std::map<MElement *, GEntity *> &bndEl2Ent,
45 const std::set<MElement *> &els, std::set<MVertex *> &toFix,
46 const std::set<MElement *> &bndEls, bool fixBndNodes);
47
48 // Mesh entities and variables
dim()49 inline const int &dim() { return _dim; }
nPC()50 inline const int &nPC() { return _nPC; }
nVert()51 inline int nVert() { return _vert.size(); }
nFV()52 inline int nFV() { return _freeVert.size(); }
nEl()53 inline int nEl() { return _el.size(); }
nPCFV(int iFV)54 inline const int &nPCFV(int iFV) { return _nPCFV[iFV]; }
indPCFV(int iFV,int iPC)55 inline int indPCFV(int iFV, int iPC) { return _startPCFV[iFV] + iPC; }
nPCEl(int iEl)56 inline int nPCEl(int iEl) { return _indPCEl[iEl].size(); }
nNodEl(int iEl)57 inline const int &nNodEl(int iEl) { return _nNodEl[iEl]; }
indPCEl(int iEl,int iPC)58 inline const int &indPCEl(int iEl, int iPC) { return _indPCEl[iEl][iPC]; }
el2V(int iEl,int i)59 inline const int &el2V(int iEl, int i) { return _el2V[iEl][i]; }
el2FV(int iEl,int i)60 inline const int &el2FV(int iEl, int i) { return _el2FV[iEl][i]; }
fv2V(int iFV)61 inline const int &fv2V(int iFV) { return _fv2V[iFV]; }
xyz(int iV)62 inline const SPoint3 &xyz(int iV) { return _xyz[iV]; }
ixyz(int iV)63 inline const SPoint3 &ixyz(int iV) { return _ixyz[iV]; }
uvwFV(int iFV)64 inline const SPoint3 &uvwFV(int iFV) { return _uvw[iFV]; }
65 inline void gXyz2gUvw(int iFV, const SPoint3 &gXyz, SPoint3 &gUvw);
66 inline void gXyz2gUvw(int iFV, const std::vector<SPoint3> &gXyz,
67 std::vector<SPoint3> &gUvw);
68 void pcScale(
69 int iFV,
70 std::vector<double> &scale); // Calc. scale of param. coord. for precond.
71 void getUvw(double *it);
72 void updateMesh(const double *it);
73 void updateGEntityPositions();
74 void writeMSH(const char *filename);
75
76 // Node distance measure
77 enum LengthScaling {
78 LS_NONE,
79 LS_MAXNODEDIST,
80 LS_MAXOUTERRADIUS,
81 LS_MINEDGELENGTH
82 };
83 void initScaledNodeDispSq(LengthScaling scaling);
invLengthScaleSq()84 inline double invLengthScaleSq() { return _invLengthScaleSq; }
85 double scaledNodeDispSq(int iFV);
86 void gradScaledNodeDispSq(int iFV, std::vector<double> &gDSq);
87
88 // High-order: scaled Jacobian and metric measures, distance to CAD
nBezEl(int iEl)89 inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
indGSJ(int iEl,int l,int iPC)90 inline int indGSJ(int iEl, int l, int iPC) { return iPC * _nBezEl[iEl] + l; }
nBndEl()91 inline int nBndEl() { return _bndEl.size(); }
nNodBndEl(int iBndEl)92 inline int nNodBndEl(int iBndEl) { return _bndEl2V[iBndEl].size(); }
bndEl2FV(int iBndEl,int i)93 inline const int &bndEl2FV(int iBndEl, int i) { return _bndEl2FV[iBndEl][i]; }
94 void initScaledJac();
95 void scaledJacAndGradients(int iEl, std::vector<double> &sJ,
96 std::vector<double> &gSJ);
97 void initMetricMin();
98 void metricMinAndGradients(int iEl, std::vector<double> &sJ,
99 std::vector<double> &gSJ);
100 bool bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF,
101 double eps);
102 void initScaledCADDistSq(double refCADDist);
103 void scaledCADDistSqAndGradients(int iBndEl, double &scaledDist,
104 std::vector<double> &gradScaledDist);
105
106 // Mesh quality
nIJacEl(int iEl)107 inline const int &nIJacEl(int iEl) { return _nIJacEl[iEl]; }
nICNEl(int iEl)108 inline const int &nICNEl(int iEl) { return _nICNEl[iEl]; }
indGIJac(int iEl,int l,int iPC)109 inline int indGIJac(int iEl, int l, int iPC)
110 {
111 return iPC * _nIJacEl[iEl] + l;
112 }
indGICN(int iEl,int l,int iPC)113 inline int indGICN(int iEl, int l, int iPC) { return iPC * _nICNEl[iEl] + l; }
114 void initIdealJac();
115 void idealJacAndGradients(int iEl, std::vector<double> &NCJ,
116 std::vector<double> &gNCJ);
117 void initInvCondNum();
118 void invCondNumAndGradients(int iEl, std::vector<double> &condNum,
119 std::vector<double> &gCondNum);
120
121 private:
122 // Mesh entities and variables
123 int _dim;
124 int _nPC; // Total nb. of parametric coordinates
125 std::vector<MElement *> _el; // List of elements
126 std::vector<GEntity *>
127 _gEnt; // Geometric entity corresponding to each element
128 std::vector<MVertex *> _vert, _freeVert; // List of vert., free vert.
129 std::vector<int> _fv2V; // Index of free vert. -> index of vert.
130 std::vector<SPoint3> _xyz,
131 _ixyz; // Physical coord. of ALL vertices (current, straight, init.)
132 std::vector<SPoint3> _uvw,
133 _iuvw; // Parametric coord. of FREE vertices (current, straight, init.)
134 std::vector<int>
135 _startPCFV; // Start index of parametric coordinates for a free vertex
136 std::vector<int> _nPCFV; // Number of parametric coordinates for a free vertex
137 std::vector<std::vector<int> > _el2FV,
138 _el2V; // Free vertices, vertices in element
139 std::vector<int> _nNodEl; // Number of Bezier poly. and nodes for an el.
140 std::vector<std::vector<int> >
141 _indPCEl; // Index of parametric coord. for an el.
142 std::vector<VertexCoord *> _coordFV; // Parametrization for a free vertex
143 int addVert(MVertex *vert);
144 int addFreeVert(MVertex *vert, const int iV, const int nPCV,
145 VertexCoord *param, std::set<MVertex *> &toFix);
indJB2DBase(int nNod,int l,int i,int j)146 static inline int indJB2DBase(int nNod, int l, int i, int j)
147 {
148 return (l * nNod + i) * nNod + j;
149 }
indJB2D(int iEl,int l,int i,int j)150 inline int indJB2D(int iEl, int l, int i, int j)
151 {
152 return indJB2DBase(_nNodEl[iEl], l, i, j);
153 }
indJB3DBase(int nNod,int l,int i,int j,int m)154 static inline int indJB3DBase(int nNod, int l, int i, int j, int m)
155 {
156 return ((l * nNod + i) * nNod + j) * nNod + m;
157 }
indJB3D(int iEl,int l,int i,int j,int m)158 inline int indJB3D(int iEl, int l, int i, int j, int m)
159 {
160 return indJB3DBase(_nNodEl[iEl], l, i, j, m);
161 }
162
163 // Node displacement
164 LengthScaling _typeLengthScale;
165 double _invLengthScaleSq; // Square inverse of a length for node displacement
166 // scaling
167
168 // High-order: scaled Jacobian and metric measures
169 enum NormalScaling { NS_UNIT, NS_INVNORM, NS_SQRTNORM };
170 std::vector<int> _nBezEl; // Number of Bezier poly. for an el.
171 std::vector<fullMatrix<double> >
172 _jacNormEl; // Normals to 2D elements for Jacobian regularization and
173 // scaling
174 std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
175 std::vector<MElement *> _bndEl; // Boundary elements
176 std::vector<std::vector<int> > _bndEl2V,
177 _bndEl2FV; // Vertices & corresponding free vertices on the boundary
178 // elements
179 std::vector<GEntity *>
180 _bndEl2Ent; // Geometric entities corresponding to the boundary elements
181 double _invRefCADDistSq;
182 void calcNormalEl2D(int iEl, NormalScaling scaling,
183 fullMatrix<double> &elNorm, bool ideal);
184
185 // Mesh quality
186 std::vector<int> _nIJacEl; // Number of NCJ values for an el
187 std::vector<fullMatrix<double> >
188 _IJacNormEl; // Normals to 2D elements for Jacobian regularization and
189 // scaling
190 std::vector<double> _invIJac; // Initial Jacobians for 3D elements
191 std::vector<int> _nICNEl; // Number of inv. cond. number values for an el.
192 std::vector<fullMatrix<double> >
193 _condNormEl; // Normals to 2D elements for inverse conditioning computation
194 };
195
gXyz2gUvw(int iFV,const SPoint3 & gXyz,SPoint3 & gUvw)196 inline void Patch::gXyz2gUvw(int iFV, const SPoint3 &gXyz, SPoint3 &gUvw)
197 {
198 _coordFV[iFV]->gXyz2gUvw(_uvw[iFV], gXyz, gUvw);
199 }
200
gXyz2gUvw(int iFV,const std::vector<SPoint3> & gXyz,std::vector<SPoint3> & gUvw)201 inline void Patch::gXyz2gUvw(int iFV, const std::vector<SPoint3> &gXyz,
202 std::vector<SPoint3> &gUvw)
203 {
204 _coordFV[iFV]->gXyz2gUvw(_uvw[iFV], gXyz, gUvw);
205 }
206
207 #endif
208