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