FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
mesh_quality_heuristic.hpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5#pragma once
6
7#include <kernel/base_header.hpp>
8#include <kernel/util/math.hpp>
9#include <kernel/util/tiny_algebra.hpp>
10
11namespace FEAT
12{
13 namespace Geometry
14 {
15
27 template<typename Shape_>
28 struct MeshQualityHeuristic
29 {
35 static String description()
36 {
37 return String("UNDEFINED");
38 }
39
58 template<typename IdxType_, typename VtxType_>
59 static typename VtxType_::CoordType compute(const IdxType_& idx, const VtxType_& vtx);
60 };
61
63 template<>
64 struct MeshQualityHeuristic<Shape::Simplex<2>>
65 {
66
67 static String description()
68 {
69 return String("max( h(T) / vol(T)^1/2):");
70 }
71
72 template<typename IdxType_, typename VtxType_>
73 static typename VtxType_::CoordType angle(const IdxType_& idx, const VtxType_& vtx,
74 typename VtxType_::CoordType* cell_worst_angle = nullptr)
75 {
76
77 typedef typename VtxType_::CoordType CoordType;
78
79 CoordType min_angle(Math::huge<CoordType>());
80 CoordType max_angle(0);
81 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
82 {
83 Index i0(idx(cell,Index(0)));
84 Index i1(idx(cell,Index(1)));
85 Index i2(idx(cell,Index(2)));
86
87 auto v1 = vtx[i1] - vtx[i0];
88 auto v2 = vtx[i2] - vtx[i0];
89 auto v3 = vtx[i2] - vtx[i1];
90
91 CoordType n1(v1.norm_euclid());
92 CoordType n2(v2.norm_euclid());
93 CoordType n3(v3.norm_euclid());
94
95 CoordType a0 = Math::acos(Tiny::dot(v1,v2)/(n1*n2));
96 CoordType a1 = Math::acos( - Tiny::dot(v3,v1)/(n3*n1));
97 CoordType a2 = Math::acos( Tiny::dot(v3,v2)/(n3*n2));
98
99 ASSERT(Math::abs(a0+a1+a2-Math::pi<CoordType>()) < Math::sqrt(Math::eps<CoordType>()));
100
101 CoordType my_min_angle(Math::min(a2,Math::min(a0,a1)));
102 CoordType my_max_angle(Math::max(a2,Math::max(a0,a1)));
103
104 if(cell_worst_angle != nullptr)
105 {
106 cell_worst_angle[cell] = CoordType(360)/(CoordType(2)*Math::pi<CoordType>())
107 *Math::min(my_min_angle, Math::abs(Math::pi<CoordType>() - my_max_angle));
108 }
109
110 min_angle = Math::min(my_min_angle, min_angle);
111 max_angle = Math::max(my_max_angle, max_angle);
112 }
113
114 CoordType worst_angle(Math::min(min_angle, Math::abs(Math::pi<CoordType>() - max_angle)));
115
116 return worst_angle*(CoordType(360))/(CoordType(2)*Math::pi<CoordType>());
117 }
118
119 template<typename IdxType_, typename VtxType_>
120 static void compute(typename VtxType_::CoordType& qual_min, typename VtxType_::CoordType& qual_sum,
121 const IdxType_& idx, const VtxType_& vtx, typename VtxType_::CoordType* cell_qual = nullptr)
122 {
123 typedef typename VtxType_::CoordType CoordType;
124
125 CoordType rho_min(Math::huge<CoordType>());
126 CoordType rho_sum(0);
127 CoordType diam(0);
128 CoordType vol(0);
129
130 Tiny::Matrix<CoordType, VtxType_::num_coords, VtxType_::num_coords> A(CoordType(0));
131
132 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
133 {
134 diam = CoordType(0);
135
136 Index i0(idx(cell,Index(0)));
137 Index i1(idx(cell,Index(1)));
138 Index i2(idx(cell,Index(2)));
139
140 A[0] = vtx[i1] - vtx[i0];
141 A[1] = vtx[i2] - vtx[i0];
142
143 vol = A.vol()/CoordType(2);
144
145 CoordType h0((vtx[i2] - vtx[i0]).norm_euclid());
146 diam = Math::max(diam, h0);
147
148 CoordType h1(A[1].norm_euclid());
149 diam = Math::max(diam, h1);
150
151 CoordType h2(A[0].norm_euclid());
152 diam = Math::max(diam, h2);
153
154 CoordType my_rho(Math::sqrt(vol)/diam);
155
156 if(cell_qual != nullptr)
157 {
158 cell_qual[cell] = Math::sqrt(CoordType(4)/Math::sqrt(CoordType(3)))*my_rho;
159 }
160
161 rho_sum += my_rho;
162 rho_min = Math::min(rho_min, my_rho);
163 }
164 // We want the quality to go to zero if rho_max goes to infinity, so we take 1/rho_max. The absolute minimum
165 // of rho is for the triangle with all angles = 60 degrees, so we scale with its rho for normalization.
166 qual_min = Math::sqrt(CoordType(4)/Math::sqrt(CoordType(3)))*rho_min;
167 qual_sum = Math::sqrt(CoordType(4)/Math::sqrt(CoordType(3)))*rho_sum;
168 }
169
170 };
171
172 template<>
173 struct MeshQualityHeuristic<Shape::Simplex<3>>
174 {
175
176 static String description()
177 {
178 return String("max( h(T) / vol(T)^1/2):");
179 }
180
181 template<typename IdxType_, typename VtxType_>
182 static typename VtxType_::CoordType angle(const IdxType_& idx, const VtxType_& vtx,
183 typename VtxType_::CoordType* cell_worst_angle = nullptr)
184 {
185
186 typedef typename VtxType_::CoordType CoordType;
187
188 Tiny::Matrix<CoordType, 6, VtxType_::num_coords> E(CoordType(0));
189 CoordType edgelengths[6];
190 CoordType angle[6];
191
192 CoordType min_angle(Math::huge<CoordType>());
193 CoordType max_angle(0);
194 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
195 {
196
197 CoordType my_min_angle(Math::huge<CoordType>());
198 CoordType my_max_angle(0);
199
200 Index i0(idx(cell,Index(0)));
201 Index i1(idx(cell,Index(1)));
202 Index i2(idx(cell,Index(2)));
203 Index i3(idx(cell,Index(3)));
204
205 E[0] = vtx[i1] - vtx[i0];
206 E[1] = vtx[i2] - vtx[i0];
207 E[2] = vtx[i3] - vtx[i0];
208 E[3] = vtx[i2] - vtx[i1];
209 E[4] = vtx[i3] - vtx[i1];
210 E[5] = vtx[i3] - vtx[i2];
211
212 for(int i(0); i < 6; ++i)
213 {
214 edgelengths[i] = E[i].norm_euclid();
215 }
216
217 angle[0] = Math::acos(Tiny::dot(E[0],E[1])/(edgelengths[0]*edgelengths[1]));
218 angle[1] = Math::acos(Tiny::dot(E[0],E[2])/(edgelengths[0]*edgelengths[2]));
219 angle[2] = Math::acos(Tiny::dot(E[1],E[2])/(edgelengths[1]*edgelengths[2]));
220 angle[3] = Math::acos(Tiny::dot(E[3],E[4])/(edgelengths[3]*edgelengths[4]));
221 angle[4] = Math::acos(Tiny::dot(E[4],E[5])/(edgelengths[4]*edgelengths[5]));
222 angle[5] = Math::acos(Tiny::dot(E[0],E[4])/(edgelengths[0]*edgelengths[4]));
223
224 for(int i(0); i < 6; ++i)
225 {
226 my_min_angle = Math::min(my_min_angle, angle[i]);
227 my_max_angle = Math::max(my_max_angle, angle[i]);
228 }
229
230 min_angle = Math::min(my_min_angle, min_angle);
231 max_angle = Math::max(my_max_angle, max_angle);
232
233 if(cell_worst_angle != nullptr)
234 {
235 cell_worst_angle[cell] = CoordType(360)/(CoordType(2)*Math::pi<CoordType>())
236 *Math::min(my_min_angle, Math::abs(Math::pi<CoordType>() - my_max_angle));
237 }
238
239 }
240
241 CoordType worst_angle(Math::min(min_angle, Math::abs(Math::pi<CoordType>() - max_angle)));
242
243 return worst_angle*(CoordType(360))/(CoordType(2)*Math::pi<CoordType>());
244 }
245
246 template<typename IdxType_, typename VtxType_>
247 static void compute(typename VtxType_::CoordType& qual_min, typename VtxType_::CoordType& qual_sum,
248 const IdxType_& idx, const VtxType_& vtx, typename VtxType_::CoordType* cell_qual = nullptr)
249 {
250 typedef typename VtxType_::CoordType CoordType;
251 const CoordType vol_ref(Math::sqrt(CoordType(3))/CoordType(4)*Math::sqrt(CoordType(6))/CoordType(3*3));
252 const CoordType fac(Math::pow(vol_ref, -CoordType(1)/CoordType(3)));
253
254 CoordType rho_min(Math::huge<CoordType>());
255 CoordType rho_sum(0);
256 CoordType diam(0);
257 CoordType vol(0);
258
259 Tiny::Matrix<CoordType, VtxType_::num_coords, VtxType_::num_coords> A(CoordType(0));
260 Tiny::Matrix<CoordType, 6, VtxType_::num_coords> E(CoordType(0));
261 CoordType edgelengths[6];
262
263 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
264 {
265 diam = CoordType(0);
266
267 Index i0(idx(cell,Index(0)));
268 Index i1(idx(cell,Index(1)));
269 Index i2(idx(cell,Index(2)));
270 Index i3(idx(cell,Index(3)));
271
272 A[0] = vtx[i1] - vtx[i0];
273 A[1] = vtx[i2] - vtx[i0];
274 A[2] = vtx[i3] - vtx[i0];
275
276 vol = A.det()/CoordType(6);
277
278 E[0] = vtx[i1] - vtx[i0];
279 E[1] = vtx[i2] - vtx[i0];
280 E[2] = vtx[i3] - vtx[i0];
281 E[3] = vtx[i2] - vtx[i1];
282 E[4] = vtx[i3] - vtx[i1];
283 E[5] = vtx[i3] - vtx[i2];
284
285 for(int i(0); i < 6; ++i)
286 {
287 edgelengths[i] = E[i].norm_euclid();
288 diam = Math::max(diam, edgelengths[i]);
289 }
290
291 CoordType my_rho(Math::pow(vol, CoordType(1)/CoordType(3))/diam);
292
293 if(cell_qual != nullptr)
294 {
295 cell_qual[cell] = fac*my_rho;
296 }
297
298 rho_sum += my_rho;
299 rho_min = Math::min(rho_min, my_rho);
300 }
301 // We want the quality to go to zero if rho_max goes to infinity, so we take 1/rho_max. The absolute minimum
302 // of rho is for the triangle with all angles = 60 degrees, so we scale with its rho for normalization.
303 qual_min = fac*rho_min;
304 qual_sum = fac*rho_sum;
305 }
306
307 };
308
309 template<>
310 struct MeshQualityHeuristic<Shape::Hypercube<2>>
311 {
312
313 static String description()
314 {
315 return String("h_min/h_max*sqrt(1 - cos(alpha_min)):");
316 }
317
318 template<typename IdxType_, typename VtxType_>
319 static typename VtxType_::CoordType angle(const IdxType_& idx, const VtxType_& vtx,
320 typename VtxType_::CoordType* cell_worst_angle = nullptr)
321 {
322
323 typedef typename VtxType_::CoordType CoordType;
324
325 CoordType min_angle(Math::huge<CoordType>());
326 CoordType max_angle(0);
327 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
328 {
329
330 CoordType my_min_angle(Math::huge<CoordType>());
331 CoordType my_max_angle(0);
332
333 Index i0(idx(cell,Index(0)));
334 Index i1(idx(cell,Index(1)));
335 Index i2(idx(cell,Index(2)));
336 Index i3(idx(cell,Index(3)));
337
338 auto v0 = vtx[i1] - vtx[i0];
339 auto v1 = vtx[i3] - vtx[i2];
340 auto v2 = vtx[i2] - vtx[i0];
341 auto v3 = vtx[i3] - vtx[i1];
342
343 CoordType h0(v0.norm_euclid());
344 CoordType h1(v1.norm_euclid());
345 CoordType h2(v2.norm_euclid());
346 CoordType h3(v3.norm_euclid());
347
348 CoordType a0 = Math::acos(Tiny::dot(v2,v0)/(h2*h0));
349 my_min_angle = Math::min(my_min_angle, Math::abs(a0));
350 my_max_angle = Math::max(my_max_angle, Math::abs(a0));
351
352 CoordType a1 = Math::acos( - Tiny::dot(v0,v3)/(h0*h3));
353 my_min_angle = Math::min(my_min_angle, Math::abs(a1));
354 my_max_angle = Math::max(my_max_angle, Math::abs(a1));
355
356 CoordType a2 = Math::acos( - Tiny::dot(v1,v2)/(h1*h2));
357 my_min_angle = Math::min(my_min_angle, Math::abs(a2));
358 my_max_angle = Math::max(my_max_angle, Math::abs(a2));
359
360 CoordType a3 = Math::acos( Tiny::dot(v3,v1)/(h3*h1));
361 my_min_angle = Math::min(my_min_angle, Math::abs(a3));
362 my_max_angle = Math::max(my_max_angle, Math::abs(a3));
363
364 min_angle = Math::min(my_min_angle, min_angle);
365 max_angle = Math::max(my_max_angle, max_angle);
366
367 if(cell_worst_angle != nullptr)
368 {
369 cell_worst_angle[cell] = CoordType(360)/(CoordType(2)*Math::pi<CoordType>())
370 *Math::min(my_min_angle, Math::abs(Math::pi<CoordType>() - my_max_angle));
371 }
372
373 ASSERT(Math::abs(a0+a1+a2+a3-CoordType(2)*Math::pi<CoordType>()) < Math::sqrt(Math::eps<CoordType>()));
374 }
375
376 CoordType worst_angle(Math::min(min_angle, Math::abs(Math::pi<CoordType>() - max_angle)));
377
378 return worst_angle*CoordType(360)/(CoordType(2)*Math::pi<CoordType>());
379 }
380
381 template<typename IdxType_, typename VtxType_>
382 static void compute(typename VtxType_::CoordType& qual_min, typename VtxType_::CoordType& qual_sum,
383 const IdxType_& idx, const VtxType_& vtx, typename VtxType_::CoordType* cell_qual = nullptr)
384 {
385
386 typedef typename VtxType_::CoordType CoordType;
387
388 qual_min = Math::huge<CoordType>();
389 qual_sum = CoordType(0);
390
391 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
392 {
393 CoordType gamma(0);
394 CoordType angle(Math::huge<CoordType>());
395
396 Index i0(idx(cell,Index(0)));
397 Index i1(idx(cell,Index(1)));
398 Index i2(idx(cell,Index(2)));
399 Index i3(idx(cell,Index(3)));
400
401 auto v0 = vtx[i1] - vtx[i0];
402 auto v1 = vtx[i3] - vtx[i2];
403 auto v2 = vtx[i2] - vtx[i0];
404 auto v3 = vtx[i3] - vtx[i1];
405
406 CoordType h_min(Math::huge<CoordType>());
407 CoordType h_max(Math::eps<CoordType>());
408
409 CoordType h0(v0.norm_euclid());
410 h_min = Math::min(h_min, h0);
411 h_max = Math::max(h_max, h0);
412
413 CoordType h1(v1.norm_euclid());
414 h_min = Math::min(h_min, h1);
415 h_max = Math::max(h_max, h1);
416
417 CoordType h2(v2.norm_euclid());
418 h_min = Math::min(h_min, h2);
419 h_max = Math::max(h_max, h2);
420
421 CoordType h3(v3.norm_euclid());
422 h_min = Math::min(h_min, h3);
423 h_max = Math::max(h_max, h3);
424
425 CoordType a0 = Math::acos(Tiny::dot(v2,v0)/(h2*h0));
426 angle = Math::min(angle, Math::abs(a0));
427 gamma = Math::max(gamma, Tiny::dot(v2,v0)/(h2*h0));
428
429 CoordType a1 = Math::acos( - Tiny::dot(v0,v3)/(h0*h3));
430 angle = Math::min(angle, Math::abs(a1));
431 gamma = Math::max(gamma, - Tiny::dot(v0,v3)/(h0*h3));
432
433 CoordType a2 = Math::acos( - Tiny::dot(v1,v2)/(h1*h2));
434 angle = Math::min(angle, Math::abs(a2));
435 gamma = Math::max(gamma, - Tiny::dot(v1,v2)/(h1*h2));
436
437 CoordType a3 = Math::acos( Tiny::dot(v3,v1)/(h3*h1));
438 angle = Math::min(angle, Math::abs(a3));
439 gamma = Math::max(gamma, Tiny::dot(v3,v1)/(h3*h1));
440
441 ASSERT(Math::abs(a0+a1+a2+a3-CoordType(2)*Math::pi<CoordType>()) < Math::sqrt(Math::eps<CoordType>()));
442
443 CoordType my_qual(h_min/h_max*Math::sqrt(CoordType(1) - gamma));
444
445 if(cell_qual != nullptr)
446 {
447 cell_qual[cell] = my_qual;
448 }
449
450 qual_min = Math::min(qual_min, my_qual);
451 qual_sum += my_qual;
452
453 }
454 }
455 };
456 template<>
457 struct MeshQualityHeuristic<Shape::Hypercube<3>>
458 {
459
460
461 static String description()
462 {
463 return String("h_min*min(det( nabla Phi))/(h_max*max(det nabla Phi)):");
464 }
465
466 template<typename IdxType_, typename VtxType_>
467 static typename VtxType_::CoordType angle(const IdxType_& idx, const VtxType_& vtx,
468 typename VtxType_::CoordType* cell_worst_angle = nullptr)
469 {
470
471 typedef typename VtxType_::CoordType CoordType;
472
473 Tiny::Matrix<CoordType, 12, VtxType_::num_coords> E(CoordType(0));
474 CoordType edgelengths[12];
475 CoordType angle[24];
476
477 CoordType min_angle(Math::huge<CoordType>());
478 CoordType max_angle(0);
479 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
480 {
481
482 CoordType my_min_angle(Math::huge<CoordType>());
483 CoordType my_max_angle(0);
484
485 E[0] = vtx[idx(cell,Index(1))] - vtx[idx(cell,Index(0))];
486 E[1] = vtx[idx(cell,Index(3))] - vtx[idx(cell,Index(2))];
487 E[2] = vtx[idx(cell,Index(5))] - vtx[idx(cell,Index(4))];
488 E[3] = vtx[idx(cell,Index(7))] - vtx[idx(cell,Index(6))];
489 E[4] = vtx[idx(cell,Index(2))] - vtx[idx(cell,Index(0))];
490 E[5] = vtx[idx(cell,Index(3))] - vtx[idx(cell,Index(1))];
491 E[6] = vtx[idx(cell,Index(6))] - vtx[idx(cell,Index(4))];
492 E[7] = vtx[idx(cell,Index(7))] - vtx[idx(cell,Index(5))];
493 E[8] = vtx[idx(cell,Index(4))] - vtx[idx(cell,Index(0))];
494 E[9] = vtx[idx(cell,Index(5))] - vtx[idx(cell,Index(1))];
495 E[10] = vtx[idx(cell,Index(6))] - vtx[idx(cell,Index(2))];
496 E[11] = vtx[idx(cell,Index(7))] - vtx[idx(cell,Index(3))];
497
498 for(int i(0); i < 12; ++i)
499 {
500 edgelengths[i] = E[i].norm_euclid();
501 }
502
503 // Angles at vertex 0
504 angle[0] = Math::acos(Tiny::dot(E[0],E[4])/(edgelengths[0]*edgelengths[4]));
505 angle[1] = Math::acos(Tiny::dot(E[0],E[8])/(edgelengths[0]*edgelengths[8]));
506 angle[2] = Math::acos(Tiny::dot(E[4],E[8])/(edgelengths[4]*edgelengths[8]));
507
508 // Angles at vertex 1
509 angle[3] = Math::acos(-Tiny::dot(E[0],E[5])/(edgelengths[0]*edgelengths[5]));
510 angle[4] = Math::acos(-Tiny::dot(E[0],E[9])/(edgelengths[9]*edgelengths[9]));
511 angle[5] = Math::acos( Tiny::dot(E[5],E[9])/(edgelengths[5]*edgelengths[9]));
512
513 // Angles at vertex 2
514 angle[6] = Math::acos(-Tiny::dot(E[1],E[4])/(edgelengths[1]*edgelengths[ 4]));
515 angle[7] = Math::acos( Tiny::dot(E[1],E[10])/(edgelengths[1]*edgelengths[10]));
516 angle[8] = Math::acos(-Tiny::dot(E[4],E[10])/(edgelengths[4]*edgelengths[10]));
517
518 // Angles at vertex 3
519 angle[9] = Math::acos(-Tiny::dot(E[1],E[ 5])/(edgelengths[1]*edgelengths[ 5]));
520 angle[10] = Math::acos(-Tiny::dot(E[1],E[11])/(edgelengths[1]*edgelengths[11]));
521 angle[11] = Math::acos(-Tiny::dot(E[5],E[11])/(edgelengths[5]*edgelengths[11]));
522
523 // Angles at vertex 4
524 angle[12] = Math::acos( Tiny::dot(E[2],E[6])/(edgelengths[2]*edgelengths[6]));
525 angle[13] = Math::acos(-Tiny::dot(E[2],E[8])/(edgelengths[2]*edgelengths[8]));
526 angle[14] = Math::acos( Tiny::dot(E[6],E[8])/(edgelengths[6]*edgelengths[8]));
527
528 // Angles at vertex 5
529 angle[15] = Math::acos(-Tiny::dot(E[2],E[7])/(edgelengths[2]*edgelengths[7]));
530 angle[16] = Math::acos(-Tiny::dot(E[2],E[9])/(edgelengths[2]*edgelengths[9]));
531 angle[17] = Math::acos(-Tiny::dot(E[7],E[9])/(edgelengths[7]*edgelengths[9]));
532
533 // Angles at vertex 6
534 angle[18] = Math::acos(-Tiny::dot(E[3],E[ 6])/(edgelengths[3]*edgelengths[ 6]));
535 angle[19] = Math::acos(-Tiny::dot(E[3],E[10])/(edgelengths[3]*edgelengths[10]));
536 angle[20] = Math::acos(-Tiny::dot(E[6],E[10])/(edgelengths[6]*edgelengths[10]));
537
538 // Angles at vertex 7
539 angle[21] = Math::acos(Tiny::dot(E[3],E[ 7])/(edgelengths[3]*edgelengths[ 7]));
540 angle[22] = Math::acos(Tiny::dot(E[3],E[11])/(edgelengths[3]*edgelengths[11]));
541 angle[23] = Math::acos(Tiny::dot(E[7],E[11])/(edgelengths[7]*edgelengths[11]));
542
543 for(int i(0); i < 24; ++i)
544 {
545 my_min_angle = Math::min(my_min_angle,angle[i]);
546 my_max_angle = Math::max(my_max_angle,angle[i]);
547 }
548
549 min_angle = Math::min(my_min_angle, min_angle);
550 max_angle = Math::max(my_max_angle, max_angle);
551
552 if(cell_worst_angle != nullptr)
553 {
554 cell_worst_angle[cell] = CoordType(360)/(CoordType(2)*Math::pi<CoordType>())
555 *Math::min(my_min_angle, Math::abs(Math::pi<CoordType>() - my_max_angle));
556 }
557
558 }
559
560 CoordType worst_angle(Math::min(min_angle, Math::abs(Math::pi<CoordType>() - max_angle)));
561
562 return worst_angle*(CoordType(360))/(CoordType(2)*Math::pi<CoordType>());
563 }
564
565 template<typename IdxType_, typename VtxType_>
566 static void compute(typename VtxType_::CoordType& qual_min, typename VtxType_::CoordType& qual_sum,
567 const IdxType_& idx, const VtxType_& vtx, typename VtxType_::CoordType* cell_qual = nullptr)
568 {
569
570 typedef typename VtxType_::CoordType CoordType;
571 typedef typename VtxType_::VertexType DomainPointType;
572
573 static constexpr int num_coords = VtxType_::num_coords;
574
575 typedef Tiny::Matrix<CoordType, 3, num_coords> JacobiMatrixType;
576
577 CoordType coeff[num_coords][8];
578 DomainPointType xq(0);
579 JacobiMatrixType jac_mat;
580 Tiny::Matrix<CoordType, 12, VtxType_::num_coords> E(CoordType(0));
581 CoordType edgelengths[12];
582
583 qual_min = Math::huge<CoordType>();
584 qual_sum = CoordType(0);
585
586 for(Index cell(0); cell < idx.get_num_entities(); ++cell)
587 {
588 CoordType h_min(Math::huge<CoordType>());
589 CoordType h_max(0);
590
591 CoordType jac_det_min(Math::huge<CoordType>());
592 CoordType jac_det_max(0);
593
594 const DomainPointType& v0 = vtx[idx(cell, 0)];
595 const DomainPointType& v1 = vtx[idx(cell, 1)];
596 const DomainPointType& v2 = vtx[idx(cell, 2)];
597 const DomainPointType& v3 = vtx[idx(cell, 3)];
598 const DomainPointType& v4 = vtx[idx(cell, 4)];
599 const DomainPointType& v5 = vtx[idx(cell, 5)];
600 const DomainPointType& v6 = vtx[idx(cell, 6)];
601 const DomainPointType& v7 = vtx[idx(cell, 7)];
602
603 E[0] = v1 - v0;
604 E[1] = v3 - v2;
605 E[2] = v5 - v4;
606 E[3] = v7 - v6;
607 E[4] = v2 - v0;
608 E[5] = v3 - v1;
609 E[6] = v6 - v4;
610 E[7] = v7 - v5;
611 E[8] = v4 - v0;
612 E[9] = v5 - v1;
613 E[10] = v6 - v2;
614 E[11] = v7 - v3;
615
616 for(int i(0); i < 12; ++i)
617 {
618 edgelengths[i] = E[i].norm_euclid();
619 h_min = Math::min(edgelengths[i], h_min);
620 h_max = Math::max(edgelengths[i], h_max);
621 }
622
623 // calculate transformation coefficients
624 // j = _coeff[i][j] for all j = 0....7
625 // v = 0 + 1*x + 2*y + 3*z + x*y*4 + x*z*5 + y*z*6 + x*y*z*7
626 for(int i(0); i < num_coords; ++i)
627 {
628 coeff[i][0] = CoordType(0.125) * CoordType( + v0[i] + v1[i] + v2[i] + v3[i] + v4[i] + v5[i] + v6[i] + v7[i]);
629 coeff[i][1] = CoordType(0.125) * CoordType( - v0[i] + v1[i] - v2[i] + v3[i] - v4[i] + v5[i] - v6[i] + v7[i]);
630 coeff[i][2] = CoordType(0.125) * CoordType( - v0[i] - v1[i] + v2[i] + v3[i] - v4[i] - v5[i] + v6[i] + v7[i]);
631 coeff[i][3] = CoordType(0.125) * CoordType( - v0[i] - v1[i] - v2[i] - v3[i] + v4[i] + v5[i] + v6[i] + v7[i]);
632 coeff[i][4] = CoordType(0.125) * CoordType( + v0[i] - v1[i] - v2[i] + v3[i] + v4[i] - v5[i] - v6[i] + v7[i]);
633 coeff[i][5] = CoordType(0.125) * CoordType( + v0[i] - v1[i] + v2[i] - v3[i] - v4[i] + v5[i] - v6[i] + v7[i]);
634 coeff[i][6] = CoordType(0.125) * CoordType( + v0[i] + v1[i] - v2[i] - v3[i] - v4[i] - v5[i] + v6[i] + v7[i]);
635 coeff[i][7] = CoordType(0.125) * CoordType( - v0[i] + v1[i] + v2[i] - v3[i] + v4[i] - v5[i] - v6[i] + v7[i]);
636 }
637
638 for(int i(0); i < 8; ++i)
639 {
640 // Set point coords using bloody bitshifts
641 for(int j(0); j < 3; ++j)
642 {
643 xq(j) = (CoordType(((i >> j) & 1) << 1) - CoordType(1)) * CoordType(1);
644 }
645
646 for(int j(0); j < jac_mat.m; ++j)
647 {
648 jac_mat(j,0) = coeff[j][1] + xq[1] * coeff[j][4] + xq[2] * (coeff[j][5] + xq[1] * coeff[j][7]);
649 jac_mat(j,1) = coeff[j][2] + xq[0] * coeff[j][4] + xq[2] * (coeff[j][6] + xq[0] * coeff[j][7]);
650 jac_mat(j,2) = coeff[j][3] + xq[0] * coeff[j][5] + xq[1] * (coeff[j][6] + xq[0] * coeff[j][7]);
651 }
652
653 CoordType jac_det(jac_mat.vol());
654 jac_det_min = Math::min(jac_det, jac_det_min);
655 jac_det_max = Math::max(jac_det, jac_det_max);
656
657 }
658
659 CoordType my_qual(h_min*jac_det_min/(h_max*jac_det_max));
660
661 if(cell_qual != nullptr)
662 {
663 cell_qual[cell] = my_qual;
664 }
665
666 qual_min = Math::min(qual_min, my_qual);
667 qual_sum += my_qual;
668
669 }
670 }
671 };
673 } // namespace Geometry
674} // namespace FEAT
ASSERT
#define ASSERT(expr)
Debug-Assertion macro definition.
Definition: assertion.hpp:229
base_header.hpp
FEAT Kernel base header.
FEAT::String
String class implementation.
Definition: string.hpp:47
FEAT::Math::sqrt
T_ sqrt(T_ x)
Returns the square-root of a value.
Definition: math.hpp:300
FEAT::Math::abs
T_ abs(T_ x)
Returns the absolute value.
Definition: math.hpp:275
FEAT::Math::acos
T_ acos(T_ x)
Returns the arccosine of a value.
Definition: math.hpp:967
FEAT::Math::pow
T_ pow(T_ x, T_ y)
Returns x raised to the power of y.
Definition: math.hpp:643
FEAT::Math::min
T_ min(T_ a, T_ b)
Returns the minimum of two values.
Definition: math.hpp:123
FEAT::Math::max
T_ max(T_ a, T_ b)
Returns the maximum of two values.
Definition: math.hpp:137
FEAT::Tiny::dot
CUDA_HOST_DEVICE T_ dot(const T_ &a, const T_ &b)
Computes the dot-product of two scalars.
Definition: tiny_algebra.hpp:3001
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122
FEAT::TrafoTags::jac_mat
@ jac_mat
specifies whether the trafo should supply jacobian matrices
FEAT::TrafoTags::jac_det
@ jac_det
specifies whether the trafo should supply jacobian determinants
FEAT::Geometry::MeshQualityHeuristic
Helper class for computing heuristic mesh quality.
Definition: mesh_quality_heuristic.hpp:29
FEAT::Geometry::MeshQualityHeuristic::description
static String description()
Returns a descriptive String.
Definition: mesh_quality_heuristic.hpp:35
FEAT::Geometry::MeshQualityHeuristic::compute
static VtxType_::CoordType compute(const IdxType_ &idx, const VtxType_ &vtx)
Computes minimum cell quality.