1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <vector>
7 #include <algorithm>
8 #include "GmshConfig.h"
9 #include "MFace.h"
10 #include "Numeric.h"
11 #include "ElementType.h"
12 #include "nodalBasis.h"
13 #include "BasisFactory.h"
14 
compare(const MVertex * const v0,const MVertex * const v1)15 bool compare(const MVertex *const v0, const MVertex *const v1)
16 {
17   return v0->getNum() < v1->getNum();
18 }
19 
sortVertices(const std::vector<MVertex * > & v,std::vector<char> & s)20 void sortVertices(const std::vector<MVertex *> &v, std::vector<char> &s)
21 {
22   if(v.size() == 3) {
23     s.resize(3);
24     if(v[0]->getNum() < v[1]->getNum() && v[0]->getNum() < v[2]->getNum()) {
25       s[0] = 0;
26       s[1] = 1;
27       s[2] = 2;
28     }
29     else if(v[1]->getNum() < v[0]->getNum() &&
30             v[1]->getNum() < v[2]->getNum()) {
31       s[0] = 1;
32       s[1] = 0;
33       s[2] = 2;
34     }
35     else {
36       s[0] = 2;
37       s[1] = 0;
38       s[2] = 1;
39     }
40 
41     if(v[s[2]]->getNum() < v[s[1]]->getNum()) {
42       char temp = s[1];
43       s[1] = s[2];
44       s[2] = temp;
45     }
46     return;
47   }
48   else if(v.size() == 4) {
49     // avoid overhead of general case below
50     MVertex * sorted[4] {v[0], v[1], v[2], v[3]};
51     std::sort(&sorted[0], &sorted[4], compare);
52     s.reserve(4);
53     for(int i = 0; i < 4; ++i) {
54       s.push_back(
55         std::distance(v.begin(), std::find(v.begin(), v.end(), sorted[i])));
56     }
57     return;
58   }
59 
60   std::vector<MVertex *> sorted = v;
61   std::sort(sorted.begin(), sorted.end(), compare);
62   s.reserve(sorted.size());
63   for(std::size_t i = 0; i < sorted.size(); i++)
64     s.push_back(
65       std::distance(v.begin(), std::find(v.begin(), v.end(), sorted[i])));
66 }
67 
MFace(MVertex * v0,MVertex * v1,MVertex * v2,MVertex * v3)68 MFace::MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3)
69 {
70   _v.reserve(v3 ? 4 : 3);
71   _v.push_back(v0);
72   _v.push_back(v1);
73   _v.push_back(v2);
74   if(v3) _v.push_back(v3);
75   sortVertices(_v, _si);
76 }
77 
MFace(const std::vector<MVertex * > & v)78 MFace::MFace(const std::vector<MVertex *> &v)
79 {
80   _v.reserve(v.size());
81   for(std::size_t i = 0; i < v.size(); i++) _v.push_back(v[i]);
82   sortVertices(_v, _si);
83 }
getOrientationFlagForFace(std::vector<int> & faceOrientationFlag)84 void MFace::getOrientationFlagForFace(std::vector<int> &faceOrientationFlag)
85 {
86   // Reference : "Higher-Order Finite Element Methods"; Pavel Solin, Karel
87   // Segeth, Ivo Dolezel, Chapman and Hall/CRC (2003)
88   if(_v.size() == 3) { // triangular face
89     if(_v[int(_si[0])]->getNum() == _v[0]->getNum()) {
90       faceOrientationFlag[0] = 0;
91       if(_v[int(_si[1])]->getNum() == _v[1]->getNum()) {
92         faceOrientationFlag[1] = 1;
93       }
94       else {
95         faceOrientationFlag[1] = -1;
96       }
97     }
98     else {
99       if(_v[1]->getNum() == _v[int(_si[0])]->getNum()) {
100         faceOrientationFlag[0] = 1;
101         if(_v[0]->getNum() == _v[int(_si[2])]->getNum()) {
102           faceOrientationFlag[1] = 1;
103         }
104         else {
105           faceOrientationFlag[1] = -1;
106         }
107       }
108       else {
109         faceOrientationFlag[0] = 2;
110         if(_v[1]->getNum() == _v[int(_si[2])]->getNum()) {
111           faceOrientationFlag[1] = 1;
112         }
113         else {
114           faceOrientationFlag[1] = -1;
115         }
116       }
117     }
118   }
119   else { // quadrilateral face
120     int c = 0;
121     for(int i = 0; i < 4; i++) {
122       if(_v[int(_si[0])]->getNum() == unsigned(_v[i]->getNum())) { c = i; }
123     }
124     int indexopposedVertex = 0;
125     switch(c) {
126     case(0): indexopposedVertex = 3; break;
127     case(1): indexopposedVertex = 2; break;
128     case(2): indexopposedVertex = 1; break;
129     case(3): indexopposedVertex = 0; break;
130     }
131     int numVertexOpposed = _v[indexopposedVertex]->getNum();
132 
133     int axis1A = _v[int(_si[0])]->getNum();
134     int axis1B = 0;
135     if(_v[int(_si[1])]->getNum() == unsigned(numVertexOpposed)) {
136       axis1B = _v[int(_si[2])]->getNum();
137     }
138     else {
139       axis1B = _v[int(_si[1])]->getNum();
140     }
141     if(unsigned(axis1A) == _v[0]->getNum()) {
142       if(unsigned(axis1B) == _v[1]->getNum()) {
143         faceOrientationFlag[0] = 1;
144         faceOrientationFlag[1] = 1;
145         faceOrientationFlag[2] = 1;
146       }
147       else {
148         faceOrientationFlag[0] = 1;
149         faceOrientationFlag[1] = 1;
150         faceOrientationFlag[2] = -1;
151       }
152     }
153     else {
154       if(unsigned(axis1A) == _v[1]->getNum()) {
155         if(unsigned(axis1B) == _v[0]->getNum()) {
156           faceOrientationFlag[0] = -1;
157           faceOrientationFlag[1] = 1;
158           faceOrientationFlag[2] = 1;
159         }
160         else {
161           faceOrientationFlag[0] = -1;
162           faceOrientationFlag[1] = 1;
163           faceOrientationFlag[2] = -1;
164         }
165       }
166       else {
167         if(unsigned(axis1A) == _v[2]->getNum()) {
168           if(unsigned(axis1B) == _v[3]->getNum()) {
169             faceOrientationFlag[0] = 1;
170             faceOrientationFlag[1] = -1;
171             faceOrientationFlag[2] = 1;
172           }
173           else {
174             faceOrientationFlag[0] = 1;
175             faceOrientationFlag[1] = -1;
176             faceOrientationFlag[2] = -1;
177           }
178         }
179         else {
180           if(unsigned(axis1B) == _v[2]->getNum()) {
181             faceOrientationFlag[0] = -1;
182             faceOrientationFlag[1] = -1;
183             faceOrientationFlag[2] = 1;
184           }
185           else {
186             faceOrientationFlag[0] = -1;
187             faceOrientationFlag[1] = -1;
188             faceOrientationFlag[2] = -1;
189           }
190         }
191       }
192     }
193   }
194 }
approximateArea() const195 double MFace::approximateArea() const
196 {
197   SPoint3 p0 = _v[0]->point(), p1 = _v[1]->point(), p2 = _v[2]->point();
198   double a = triangle_area(p0, p1, p2);
199   if(_v.size() == 3) return a;
200   a += triangle_area(p0, p2, _v[3]->point());
201   return a;
202 }
203 
normal() const204 SVector3 MFace::normal() const
205 {
206   double n[3];
207   normal3points(_v[0]->x(), _v[0]->y(), _v[0]->z(), _v[1]->x(), _v[1]->y(),
208                 _v[1]->z(), _v[2]->x(), _v[2]->y(), _v[2]->z(), n);
209   return SVector3(n[0], n[1], n[2]);
210 }
211 
computeCorrespondence(const MFace & other,int & rotation,bool & swap) const212 bool MFace::computeCorrespondence(const MFace &other, int &rotation,
213                                   bool &swap) const
214 {
215   rotation = 0;
216   swap = false;
217 
218   if(*this == other) {
219     for(std::size_t i = 0; i < getNumVertices(); i++) {
220       if(_v[0] == other.getVertex(i)) {
221         rotation = i;
222         break;
223       }
224     }
225     if(_v[1] == other.getVertex((rotation + 1) % getNumVertices()))
226       swap = false;
227     else
228       swap = true;
229     return true;
230   }
231   return false;
232 }
233 
MFaceN(int type,int order,const std::vector<MVertex * > & v)234 MFaceN::MFaceN(int type, int order, const std::vector<MVertex *> &v)
235   : _type(type), _order(order)
236 {
237   _v.resize(v.size());
238   for(std::size_t i = 0; i < v.size(); i++) _v[i] = v[i];
239 }
240 
getHighOrderEdge(int num,int sign) const241 MEdgeN MFaceN::getHighOrderEdge(int num, int sign) const
242 {
243   int nCorner = getNumCorners();
244   std::vector<MVertex *> vertices(static_cast<std::size_t>(_order) + 1);
245   if(sign == 1) {
246     vertices[0] = _v[num];
247     vertices[1] = _v[(num + 1) % nCorner];
248   }
249   else {
250     vertices[0] = _v[(num + 1) % nCorner];
251     vertices[1] = _v[num];
252   }
253   int start = nCorner + num * (_order - 1);
254   int end = nCorner + (num + 1) * (_order - 1);
255   int k = 1;
256   if(sign == 1) {
257     for(int i = start; i < end; ++i) vertices[++k] = _v[i];
258   }
259   else {
260     for(int i = end - 1; i >= start; --i) vertices[++k] = _v[i];
261   }
262   return MEdgeN(vertices);
263 }
264 
getFace() const265 MFace MFaceN::getFace() const
266 {
267   if(_type == TYPE_TRI)
268     return MFace(_v[0], _v[1], _v[2]);
269   else
270     return MFace(_v[0], _v[1], _v[2], _v[3]);
271 }
272 
pnt(double u,double v) const273 SPoint3 MFaceN::pnt(double u, double v) const
274 {
275   int tag = ElementType::getType(_type, _order);
276   const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
277 
278   double f[100];
279   fs->f(u, v, 0, f);
280 
281   double x = 0, y = 0, z = 0;
282   for(int j = 0; j < fs->getNumShapeFunctions(); j++) {
283     x += f[j] * _v[j]->x();
284     y += f[j] * _v[j]->y();
285     z += f[j] * _v[j]->z();
286   }
287   return SPoint3(x, y, z);
288 }
289 
tangent(double u,double v,int num) const290 SVector3 MFaceN::tangent(double u, double v, int num) const
291 {
292   if(num != 0 && num != 1) num = 0;
293 
294   int tag = ElementType::getType(_type, _order);
295   const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
296 
297   double sf[100][3];
298   fs->df(u, v, 0, sf);
299 
300   double dx = 0, dy = 0, dz = 0;
301   for(int j = 0; j < fs->getNumShapeFunctions(); j++) {
302     dx += sf[j][num] * _v[j]->x();
303     dy += sf[j][num] * _v[j]->y();
304     dz += sf[j][num] * _v[j]->z();
305   }
306   return SVector3(dx, dy, dz).unit();
307 }
308 
normal(double u,double v) const309 SVector3 MFaceN::normal(double u, double v) const
310 {
311   int tag = ElementType::getType(_type, _order);
312   const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
313 
314   double sf[100][3];
315   fs->df(u, v, 0, sf);
316 
317   double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
318   for(int j = 0; j < fs->getNumShapeFunctions(); j++) {
319     for(int k = 0; k < 1; ++k) {
320       dx[k] += sf[j][k] * _v[j]->x();
321       dy[k] += sf[j][k] * _v[j]->y();
322       dz[k] += sf[j][k] * _v[j]->z();
323     }
324   }
325 
326   SVector3 t0 = SVector3(dx[0], dy[0], dz[0]);
327   SVector3 t1 = SVector3(dx[1], dy[1], dz[1]);
328 
329   return crossprod(t0, t1).unit();
330 }
331 
frame(double u,double v,SVector3 & t0,SVector3 & t1,SVector3 & n) const332 void MFaceN::frame(double u, double v, SVector3 &t0, SVector3 &t1,
333                    SVector3 &n) const
334 {
335   int tag = ElementType::getType(_type, _order);
336   const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
337 
338   double sf[100][3];
339   fs->df(u, v, 0, sf);
340 
341   double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
342   for(int j = 0; j < fs->getNumShapeFunctions(); j++) {
343     for(int k = 0; k < 2; ++k) {
344       dx[k] += sf[j][k] * _v[j]->x();
345       dy[k] += sf[j][k] * _v[j]->y();
346       dz[k] += sf[j][k] * _v[j]->z();
347     }
348   }
349 
350   t0 = SVector3(dx[0], dy[0], dz[0]).unit();
351   t1 = SVector3(dx[1], dy[1], dz[1]).unit();
352   n = crossprod(t0, t1);
353 }
354 
frame(double u,double v,SPoint3 & p,SVector3 & t0,SVector3 & t1,SVector3 & n) const355 void MFaceN::frame(double u, double v, SPoint3 &p, SVector3 &t0, SVector3 &t1,
356                    SVector3 &n) const
357 {
358   int tag = ElementType::getType(_type, _order);
359   const nodalBasis *fs = BasisFactory::getNodalBasis(tag);
360 
361   double f[100];
362   double sf[100][3];
363   fs->f(u, v, 0, f);
364   fs->df(u, v, 0, sf);
365 
366   double x = 0, y = 0, z = 0;
367   double dx[2] = {0, 0}, dy[2] = {0, 0}, dz[2] = {0, 0};
368   for(int j = 0; j < fs->getNumShapeFunctions(); j++) {
369     x += f[j] * _v[j]->x();
370     y += f[j] * _v[j]->y();
371     z += f[j] * _v[j]->z();
372     for(int k = 0; k < 2; ++k) {
373       dx[k] += sf[j][k] * _v[j]->x();
374       dy[k] += sf[j][k] * _v[j]->y();
375       dz[k] += sf[j][k] * _v[j]->z();
376     }
377   }
378 
379   p = SPoint3(x, y, z);
380   t0 = SVector3(dx[0], dy[0], dz[0]).unit();
381   t1 = SVector3(dx[1], dy[1], dz[1]).unit();
382   n = crossprod(t0, t1);
383 }
384 
repositionInnerVertices(const fullMatrix<double> * placement) const385 void MFaceN::repositionInnerVertices(const fullMatrix<double> *placement) const
386 {
387   int nCorner = getNumCorners();
388   int start = nCorner + (_order - 1) * nCorner;
389   for(int i = start; i < (int)_v.size(); ++i) {
390     MVertex *v = _v[i];
391     v->x() = 0;
392     v->y() = 0;
393     v->z() = 0;
394     for(int j = 0; j < placement->size2(); ++j) {
395       const double coeff = (*placement)(i - start, j);
396       v->x() += coeff * _v[j]->x();
397       v->y() += coeff * _v[j]->y();
398       v->z() += coeff * _v[j]->z();
399     }
400   }
401 }
402