Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in / Register
Toggle navigation
Menu
Open sidebar
3rd
verdict
Commits
9d11ce5b
Commit
9d11ce5b
authored
3 years ago
by
Spiros Tsalikis
Committed by
Mark Richardson
3 years ago
Browse files
Options
Download
Email Patches
Plain Diff
Use c++ version of math functions
parent
4a176af1
master
julin/master
1.4.2
1.4.1
No related merge requests found
Changes
11
Hide whitespace changes
Inline
Side-by-side
Showing
11 changed files
V_EdgeMetric.cpp
+2
-1
V_EdgeMetric.cpp
V_HexMetric.cpp
+63
-63
V_HexMetric.cpp
V_PyramidMetric.cpp
+3
-4
V_PyramidMetric.cpp
V_QuadMetric.cpp
+23
-23
V_QuadMetric.cpp
V_TetMetric.cpp
+52
-51
V_TetMetric.cpp
V_TriMetric.cpp
+6
-6
V_TriMetric.cpp
V_WedgeMetric.cpp
+21
-21
V_WedgeMetric.cpp
VerdictVector.cpp
+3
-4
VerdictVector.cpp
VerdictVector.hpp
+2
-2
VerdictVector.hpp
unittests/verdict.test.cpp
+24
-26
unittests/verdict.test.cpp
verdict_defines.hpp
+10
-10
verdict_defines.hpp
with
209 additions
and
211 deletions
+209
-211
V_EdgeMetric.cpp
View file @
9d11ce5b
...
...
@@ -19,6 +19,7 @@
*/
#include "verdict.h"
#include <cmath>
namespace
VERDICT_NAMESPACE
...
...
@@ -32,6 +33,6 @@ double edge_length(int /*num_nodes*/, const double coordinates[][3])
double
x
=
coordinates
[
1
][
0
]
-
coordinates
[
0
][
0
];
double
y
=
coordinates
[
1
][
1
]
-
coordinates
[
0
][
1
];
double
z
=
coordinates
[
1
][
2
]
-
coordinates
[
0
][
2
];
return
(
double
)(
sqrt
(
x
*
x
+
y
*
y
+
z
*
z
));
return
(
double
)(
std
::
sqrt
(
x
*
x
+
y
*
y
+
z
*
z
));
}
}
// namespace verdict
This diff is collapsed.
Click to expand it.
V_HexMetric.cpp
View file @
9d11ce5b
...
...
@@ -32,7 +32,7 @@ extern void quad_minimum_maximum_angle(double min_max_angles[2], const double co
static
const
double
one_third
=
1.0
/
3.0
;
static
const
double
two_thirds
=
2.0
/
3.0
;
static
const
double
sqrt3
=
sqrt
(
3.0
);
static
const
double
sqrt3
=
std
::
sqrt
(
3.0
);
//! weights based on the average size of a hex
static
int
hex_get_weight
(
...
...
@@ -47,7 +47,7 @@ static int hex_get_weight(
v2
.
set
(
0
,
1
,
0
);
v3
.
set
(
0
,
0
,
1
);
double
scale
=
pow
(
average_size
/
(
VerdictVector
::
Dot
(
v1
,
(
v2
*
v3
))),
0.33333333333
);
double
scale
=
std
::
pow
(
average_size
/
(
VerdictVector
::
Dot
(
v1
,
(
v2
*
v3
))),
0.33333333333
);
v1
*=
scale
;
v2
*=
scale
;
v3
*=
scale
;
...
...
@@ -320,8 +320,8 @@ static void localize_hex_coordinates(const double coordinates[][3], VerdictVecto
double DY = point_1256.y() - point_0374.y();
double DZ = point_1256.z() - point_0374.z();
double AMAGX = sqrt(DX*DX + DZ*DZ);
double AMAGY = sqrt(DX*DX + DY*DY + DZ*DZ);
double AMAGX =
std::
sqrt(DX*DX + DZ*DZ);
double AMAGY =
std::
sqrt(DX*DX + DY*DY + DZ*DZ);
double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
...
...
@@ -355,7 +355,7 @@ static void localize_hex_coordinates(const double coordinates[][3], VerdictVecto
DY = delta.y();
DZ = delta.z();
AMAGY = sqrt(DY*DY + DZ*DZ);
AMAGY =
std::
sqrt(DY*DY + DZ*DZ);
FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
double CX = DY / (AMAGY + FMAGY) + FMAGY;
...
...
@@ -378,13 +378,13 @@ static double safe_ratio3( const double numerator,
const double filter_n = max_ratio * 1.0e-16;
const double filter_d = 1.0e-16;
if (
f
abs( numerator ) <= filter_n &&
f
abs( denominator ) >= filter_d )
if (
std::
abs( numerator ) <= filter_n &&
std::
abs( denominator ) >= filter_d )
{
return_value = numerator / denominator;
}
else
{
return_value =
f
abs(numerator) / max_ratio >=
f
abs(denominator) ?
return_value =
std::
abs(numerator) / max_ratio >=
std::
abs(denominator) ?
( (numerator >= 0.0 && denominator >= 0.0) ||
(numerator < 0.0 && denominator < 0.0) ?
max_ratio : -max_ratio )
...
...
@@ -400,7 +400,7 @@ static double safe_ratio(const double numerator, const double denominator)
double
return_value
;
const
double
filter_n
=
VERDICT_DBL_MAX
;
const
double
filter_d
=
VERDICT_DBL_MIN
;
if
(
f
abs
(
numerator
)
<=
filter_n
&&
f
abs
(
denominator
)
>=
filter_d
)
if
(
std
::
abs
(
numerator
)
<=
filter_n
&&
std
::
abs
(
denominator
)
>=
filter_d
)
{
return_value
=
numerator
/
denominator
;
}
...
...
@@ -427,7 +427,7 @@ static double condition_comp(
double
term2
=
VerdictVector
::
Dot
((
xxi
*
xet
),
(
xxi
*
xet
))
+
VerdictVector
::
Dot
((
xet
*
xze
),
(
xet
*
xze
))
+
VerdictVector
::
Dot
((
xze
*
xxi
),
(
xze
*
xxi
));
return
sqrt
(
term1
*
term2
)
/
det
;
return
std
::
sqrt
(
term1
*
term2
)
/
det
;
}
static
double
oddy_comp
(
...
...
@@ -451,8 +451,8 @@ static double oddy_comp(
double
norm_J_squared
=
g11
+
g22
+
g33
;
oddy_metric
=
(
norm_G_squared
-
one_third
*
norm_J_squared
*
norm_J_squared
)
/
pow
(
rt_g
,
4.
*
one_third
);
oddy_metric
=
(
norm_G_squared
-
one_third
*
norm_J_squared
*
norm_J_squared
)
/
std
::
pow
(
rt_g
,
4.
*
one_third
);
}
else
{
...
...
@@ -473,84 +473,84 @@ static double hex_edge_length(int max_min, const double coordinates[][3])
temp
[
i
]
=
coordinates
[
1
][
i
]
-
coordinates
[
0
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
0
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
0
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
2
][
i
]
-
coordinates
[
1
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
1
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
1
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
3
][
i
]
-
coordinates
[
2
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
2
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
2
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
0
][
i
]
-
coordinates
[
3
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
3
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
3
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
5
][
i
]
-
coordinates
[
4
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
4
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
4
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
6
][
i
]
-
coordinates
[
5
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
5
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
5
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
7
][
i
]
-
coordinates
[
6
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
6
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
6
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
4
][
i
]
-
coordinates
[
7
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
7
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
7
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
4
][
i
]
-
coordinates
[
0
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
8
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
8
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
5
][
i
]
-
coordinates
[
1
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
9
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
9
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
6
][
i
]
-
coordinates
[
2
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
10
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
10
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
7
][
i
]
-
coordinates
[
3
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
edge
[
11
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
edge
[
11
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
double
_edge
=
edge
[
0
];
...
...
@@ -583,28 +583,28 @@ static double diag_length(int max_min, const double coordinates[][3])
temp
[
i
]
=
coordinates
[
6
][
i
]
-
coordinates
[
0
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
diag
[
0
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
diag
[
0
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
4
][
i
]
-
coordinates
[
2
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
diag
[
1
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
diag
[
1
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
7
][
i
]
-
coordinates
[
1
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
diag
[
2
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
diag
[
2
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
temp
[
i
]
=
coordinates
[
5
][
i
]
-
coordinates
[
3
][
i
];
temp
[
i
]
=
temp
[
i
]
*
temp
[
i
];
}
diag
[
3
]
=
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
diag
[
3
]
=
std
::
sqrt
(
temp
[
0
]
+
temp
[
1
]
+
temp
[
2
]);
double
diagonal
=
diag
[
0
];
if
(
max_min
==
0
)
// Return min diagonal
...
...
@@ -828,7 +828,7 @@ double hex_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
M2
=
Mab
>
Mcd
?
Mab
:
Mcd
;
M2
=
M2
>
Mef
?
M2
:
Mef
;
double
edge_ratio
=
sqrt
(
M2
/
m2
);
double
edge_ratio
=
std
::
sqrt
(
M2
/
m2
);
if
(
edge_ratio
>
0
)
{
...
...
@@ -1055,9 +1055,9 @@ double hex_skew(int /*num_nodes*/, const double coordinates[][3])
return
VERDICT_DBL_MAX
;
}
skew_1
=
f
abs
(
VerdictVector
::
Dot
(
efg1
,
efg2
));
skew_2
=
f
abs
(
VerdictVector
::
Dot
(
efg1
,
efg3
));
skew_3
=
f
abs
(
VerdictVector
::
Dot
(
efg2
,
efg3
));
skew_1
=
std
::
abs
(
VerdictVector
::
Dot
(
efg1
,
efg2
));
skew_2
=
std
::
abs
(
VerdictVector
::
Dot
(
efg1
,
efg3
));
skew_3
=
std
::
abs
(
VerdictVector
::
Dot
(
efg2
,
efg3
));
double
skew
=
std
::
max
({
skew_1
,
skew_2
,
skew_3
});
...
...
@@ -1086,9 +1086,9 @@ double hex_taper(int /*num_nodes*/, const double coordinates[][3])
VerdictVector
efg13
=
calc_hex_efg
(
13
,
node_pos
);
VerdictVector
efg23
=
calc_hex_efg
(
23
,
node_pos
);
double
taper_1
=
f
abs
(
safe_ratio
(
efg12
.
length
(),
std
::
min
(
efg1
.
length
(),
efg2
.
length
())));
double
taper_2
=
f
abs
(
safe_ratio
(
efg13
.
length
(),
std
::
min
(
efg1
.
length
(),
efg3
.
length
())));
double
taper_3
=
f
abs
(
safe_ratio
(
efg23
.
length
(),
std
::
min
(
efg2
.
length
(),
efg3
.
length
())));
double
taper_1
=
std
::
abs
(
safe_ratio
(
efg12
.
length
(),
std
::
min
(
efg1
.
length
(),
efg2
.
length
())));
double
taper_2
=
std
::
abs
(
safe_ratio
(
efg13
.
length
(),
std
::
min
(
efg1
.
length
(),
efg3
.
length
())));
double
taper_3
=
std
::
abs
(
safe_ratio
(
efg23
.
length
(),
std
::
min
(
efg2
.
length
(),
efg3
.
length
())));
double
taper
=
std
::
max
({
taper_1
,
taper_2
,
taper_3
});
...
...
@@ -1469,7 +1469,7 @@ double hex_dimension(int /*num_nodes*/, const double coordinates[][3])
SQR
(
gradop
[
1
][
3
])
+
SQR
(
gradop
[
2
][
3
])
+
SQR
(
gradop
[
3
][
3
])
+
SQR
(
gradop
[
4
][
3
])
+
SQR
(
gradop
[
5
][
3
])
+
SQR
(
gradop
[
6
][
3
])
+
SQR
(
gradop
[
7
][
3
])
+
SQR
(
gradop
[
8
][
3
]));
return
(
double
)
sqrt
(
aspect
);
return
(
double
)
std
::
sqrt
(
aspect
);
}
/*!
...
...
@@ -2016,7 +2016,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
VerdictVector xet(jacobian[1]);
VerdictVector xze(jacobian[2]);
jacobi = xxi % ( xet * xze );
double scale = sqrt(xxi.length_squared() * xet.length_squared() * xze.length_squared());
double scale =
std::
sqrt(xxi.length_squared() * xet.length_squared() * xze.length_squared());
jacobi /= scale;
if(jacobi < min_norm_jac)
{
...
...
@@ -2054,7 +2054,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
...
...
@@ -2086,7 +2086,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2117,7 +2117,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2148,7 +2148,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2179,7 +2179,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2210,7 +2210,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2241,7 +2241,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2272,7 +2272,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2303,7 +2303,7 @@ double hex_scaled_jacobian(int num_nodes, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
temp_norm_jac
=
jacobi
/
lengths
;
if
(
temp_norm_jac
<
min_norm_jac
)
{
...
...
@@ -2406,7 +2406,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2430,7 +2430,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2454,7 +2454,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2478,7 +2478,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2502,7 +2502,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2526,7 +2526,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2550,7 +2550,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2574,7 +2574,7 @@ double hex_shear(int /*num_nodes*/, const double coordinates[][3])
return
0
;
}
lengths
=
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
lengths
=
std
::
sqrt
(
len1_sq
*
len2_sq
*
len3_sq
);
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
<
VERDICT_DBL_MIN
)
{
...
...
@@ -2619,7 +2619,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2640,7 +2640,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2661,7 +2661,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2682,7 +2682,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2703,7 +2703,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2724,7 +2724,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2745,7 +2745,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -2766,7 +2766,7 @@ double hex_shape(int /*num_nodes*/, const double coordinates[][3])
det
=
VerdictVector
::
Dot
(
xxi
,
(
xet
*
xze
));
if
(
det
>
VERDICT_DBL_MIN
)
{
shape
=
3
*
pow
(
det
,
two_thirds
)
/
shape
=
3
*
std
::
pow
(
det
,
two_thirds
)
/
(
VerdictVector
::
Dot
(
xxi
,
xxi
)
+
VerdictVector
::
Dot
(
xet
,
xet
)
+
VerdictVector
::
Dot
(
xze
,
xze
));
}
else
...
...
@@ -3036,7 +3036,7 @@ double hex_distortion(int num_nodes, const double coordinates[][3])
minimum_jacobian
=
jacobian
;
}
}
if
(
f
abs
(
element_volume
)
>
0.0
)
if
(
std
::
abs
(
element_volume
)
>
0.0
)
{
distortion
=
minimum_jacobian
/
element_volume
*
8.
;
}
...
...
@@ -3061,7 +3061,7 @@ double hex_timestep(int num_nodes, const double coordinates[][3], double density
double
char_length
=
hex_dimension
(
num_nodes
,
coordinates
);
double
M
=
youngs_modulus
*
(
1
-
poissons_ratio
)
/
((
1
-
2
*
poissons_ratio
)
*
(
1
+
poissons_ratio
));
double
denominator
=
sqrt
(
M
/
density
);
double
denominator
=
std
::
sqrt
(
M
/
density
);
return
char_length
/
denominator
;
}
...
...
This diff is collapsed.
Click to expand it.
V_PyramidMetric.cpp
View file @
9d11ce5b
...
...
@@ -24,13 +24,12 @@
#include <algorithm>
#include <array>
namespace
VERDICT_NAMESPACE
{
extern
double
quad_equiangle_skew
(
int
num_nodes
,
const
double
coordinates
[][
3
]);
extern
double
tri_equiangle_skew
(
int
num_nodes
,
const
double
coordinates
[][
3
]);
static
const
double
sqrt2_2
=
sqrt
(
2.0
)
/
2.0
;
static
const
double
sqrt2_2
=
std
::
sqrt
(
2.0
)
/
2.0
;
// local methods
void
make_pyramid_tets
(
const
double
coordinates
[][
3
],
double
tet1_coords
[][
3
],
...
...
@@ -212,7 +211,7 @@ double pyramid_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
return
0.0
;
}
double
min_jac
=
(
*
iter
)
*
2.0
/
sqrt
(
2.0
);
double
min_jac
=
(
*
iter
)
*
2.0
/
std
::
sqrt
(
2.0
);
return
min_jac
<
1.0
?
min_jac
:
1.0
-
(
min_jac
-
1.0
);
}
...
...
@@ -444,7 +443,7 @@ double largest_pyramid_edge(const double coordinates[][3])
max
=
std
::
max
(
max
,
l6
);
max
=
std
::
max
(
max
,
l7
);
return
sqrt
(
max
);
return
std
::
sqrt
(
max
);
}
double
distance_point_to_pyramid_base
(
...
...
This diff is collapsed.
Click to expand it.
V_QuadMetric.cpp
View file @
9d11ce5b
...
...
@@ -27,7 +27,7 @@
namespace
VERDICT_NAMESPACE
{
static
const
double
sqrt2
=
sqrt
(
2.0
);
static
const
double
sqrt2
=
std
::
sqrt
(
2.0
);
static
const
double
radius_ratio_normal_coeff
=
1.
/
(
2.
*
sqrt2
);
/*!
...
...
@@ -42,7 +42,7 @@ static int quad_get_weight(
m12
=
0
;
m22
=
1
;
double
scale
=
sqrt
(
average_quad_size
/
(
m11
*
m22
-
m21
*
m12
));
double
scale
=
std
::
sqrt
(
average_quad_size
/
(
m11
*
m22
-
m21
*
m12
));
m11
*=
scale
;
m21
*=
scale
;
...
...
@@ -305,19 +305,19 @@ void quad_minimum_maximum_angle(double min_max_angles[2], const double coordinat
return
;
}
angle
=
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
angle
=
std
::
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
angle
=
std
::
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
angle
=
std
::
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
angle
=
std
::
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
max_angle
=
std
::
max
(
angle
,
max_angle
);
...
...
@@ -393,7 +393,7 @@ double quad_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
}
else
{
double
edge_ratio
=
sqrt
(
M2
/
m2
);
double
edge_ratio
=
std
::
sqrt
(
M2
/
m2
);
if
(
edge_ratio
>
0
)
{
...
...
@@ -529,7 +529,7 @@ double quad_radius_ratio(int /*num_nodes*/, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MAX
;
}
double
radius_ratio
=
radius_ratio_normal_coeff
*
sqrt
((
a2
+
b2
+
c2
+
d2
)
*
h2
)
/
t0
;
double
radius_ratio
=
radius_ratio_normal_coeff
*
std
::
sqrt
((
a2
+
b2
+
c2
+
d2
)
*
h2
)
/
t0
;
if
(
radius_ratio
>
0
)
{
...
...
@@ -666,7 +666,7 @@ double quad_skew(int /*num_nodes*/, const double coordinates[][3])
return
0.0
;
}
double
skew
=
f
abs
(
principle_axes
[
0
]
%
principle_axes
[
1
]);
double
skew
=
std
::
abs
(
principle_axes
[
0
]
%
principle_axes
[
1
]);
return
(
double
)
std
::
min
(
skew
,
VERDICT_DBL_MAX
);
}
...
...
@@ -729,8 +729,8 @@ double quad_warpage(int /*num_nodes*/, const double coordinates[][3])
return
(
double
)
VERDICT_DBL_MIN
;
}
double
warpage
=
pow
(
std
::
min
(
corner_normals
[
0
]
%
corner_normals
[
2
],
corner_normals
[
1
]
%
corner_normals
[
3
]),
3
);
double
warpage
=
std
::
pow
(
std
::
min
(
corner_normals
[
0
]
%
corner_normals
[
2
],
corner_normals
[
1
]
%
corner_normals
[
3
]),
3
);
if
(
warpage
>
0
)
{
...
...
@@ -792,8 +792,8 @@ double quad_stretch(int /*num_nodes*/, const double coordinates[][3])
else
{
double
stretch
=
(
double
)(
sqrt2
*
sqrt
(
std
::
min
(
std
::
min
(
lengths_squared
[
0
],
lengths_squared
[
1
]),
std
::
min
(
lengths_squared
[
2
],
lengths_squared
[
3
]))
/
std
::
sqrt
(
std
::
min
(
std
::
min
(
lengths_squared
[
0
],
lengths_squared
[
1
]),
std
::
min
(
lengths_squared
[
2
],
lengths_squared
[
3
]))
/
diag02
));
return
(
double
)
std
::
min
(
stretch
,
VERDICT_DBL_MAX
);
...
...
@@ -841,16 +841,16 @@ double quad_maximum_angle(int /*num_nodes*/, const double coordinates[][3])
return
0.0
;
}
angle
=
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
angle
=
std
::
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
angle
=
std
::
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
angle
=
std
::
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
max_angle
=
std
::
max
(
angle
,
max_angle
);
angle
=
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
angle
=
std
::
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
max_angle
=
std
::
max
(
angle
,
max_angle
);
max_angle
=
max_angle
*
180.0
/
VERDICT_PI
;
...
...
@@ -912,16 +912,16 @@ double quad_minimum_angle(int /*num_nodes*/, const double coordinates[][3])
return
360.0
;
}
angle
=
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
angle
=
std
::
acos
(
-
(
edges
[
0
]
%
edges
[
1
])
/
(
length
[
0
]
*
length
[
1
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
angle
=
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
angle
=
std
::
acos
(
-
(
edges
[
1
]
%
edges
[
2
])
/
(
length
[
1
]
*
length
[
2
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
angle
=
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
angle
=
std
::
acos
(
-
(
edges
[
2
]
%
edges
[
3
])
/
(
length
[
2
]
*
length
[
3
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
angle
=
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
angle
=
std
::
acos
(
-
(
edges
[
3
]
%
edges
[
0
])
/
(
length
[
3
]
*
length
[
0
]));
min_angle
=
std
::
min
(
angle
,
min_angle
);
min_angle
=
min_angle
*
180.0
/
VERDICT_PI
;
...
...
@@ -1382,7 +1382,7 @@ double quad_distortion(int num_nodes, const double coordinates[][3])
for
(
ja
=
0
;
ja
<
num_nodes
;
ja
++
)
{
dot_product
=
normal_at_nodes
[
0
]
%
normal_at_nodes
[
ja
];
if
(
f
abs
(
dot_product
)
<
0.99
)
if
(
std
::
abs
(
dot_product
)
<
0.99
)
{
flat_element
=
false
;
break
;
...
...
@@ -1392,7 +1392,7 @@ double quad_distortion(int num_nodes, const double coordinates[][3])
// take into consideration the thickness of the element
double
thickness
;
// get_quad_thickness(face, element_area, thickness );
thickness
=
0.001
*
sqrt
(
element_area
);
thickness
=
0.001
*
std
::
sqrt
(
element_area
);
// set thickness gauss point location
double
zl
=
0.5773502691896
;
...
...
This diff is collapsed.
Click to expand it.
V_TetMetric.cpp
View file @
9d11ce5b
...
...
@@ -32,9 +32,9 @@ static const double one_third = 1.0 / 3.0;
static
const
double
two_thirds
=
2.0
/
3.0
;
static
const
double
one_fourth
=
1.0
/
4.0
;
static
const
double
four_ninths
=
4.0
/
9.0
;
static
const
double
sqrt2
=
sqrt
(
2.0
);
static
const
double
sqrt3
=
sqrt
(
3.0
);
static
const
double
sqrt6
=
sqrt
(
6.0
);
static
const
double
sqrt2
=
std
::
sqrt
(
2.0
);
static
const
double
sqrt3
=
std
::
sqrt
(
3.0
);
static
const
double
sqrt6
=
std
::
sqrt
(
6.0
);
static
const
double
three_times_1plussqrt3
=
3.0
*
(
1
+
sqrt3
);
static
const
double
normal_coeff
=
180.
*
.3183098861837906715377675267450287
;
static
const
double
aspect_ratio_normal_coeff
=
sqrt6
/
12.
;
...
...
@@ -98,12 +98,12 @@ double tet_equiangle_skew(int /*num_nodes*/, const double coordinates[][3])
VerdictVector
bcd
=
bc
*
cd
;
bcd
.
normalize
();
double
alpha
=
acos
(
-
(
abc
%
abd
));
double
beta
=
acos
(
-
(
abc
%
acd
));
double
gamma
=
acos
(
-
(
abc
%
bcd
));
double
delta
=
acos
(
-
(
abd
%
acd
));
double
epsilon
=
acos
(
-
(
abd
%
bcd
));
double
zeta
=
acos
(
-
(
acd
%
bcd
));
double
alpha
=
std
::
acos
(
-
(
abc
%
abd
));
double
beta
=
std
::
acos
(
-
(
abc
%
acd
));
double
gamma
=
std
::
acos
(
-
(
abc
%
bcd
));
double
delta
=
std
::
acos
(
-
(
abd
%
acd
));
double
epsilon
=
std
::
acos
(
-
(
abd
%
bcd
));
double
zeta
=
std
::
acos
(
-
(
acd
%
bcd
));
double
min_angle
=
alpha
;
min_angle
=
min_angle
<
beta
?
min_angle
:
beta
;
...
...
@@ -121,7 +121,7 @@ double tet_equiangle_skew(int /*num_nodes*/, const double coordinates[][3])
max_angle
=
max_angle
>
zeta
?
max_angle
:
zeta
;
max_angle
*=
normal_coeff
;
double
theta
=
acos
(
1
/
3.0
)
*
normal_coeff
;
// 70.528779365509308630754000660038;
double
theta
=
std
::
acos
(
1
/
3.0
)
*
normal_coeff
;
// 70.528779365509308630754000660038;
double
dihedral_skew_max
=
(
max_angle
-
theta
)
/
(
180
-
theta
);
double
dihedral_skew_min
=
(
theta
-
min_angle
)
/
theta
;
...
...
@@ -130,18 +130,18 @@ double tet_equiangle_skew(int /*num_nodes*/, const double coordinates[][3])
max_angle
=
0.0
;
double
angles
[
12
];
angles
[
0
]
=
acos
(
-
(
ab
%
bc
));
angles
[
1
]
=
acos
((
bc
%
ac
));
angles
[
2
]
=
acos
((
ac
%
ab
));
angles
[
3
]
=
acos
(
-
(
ab
%
bd
));
angles
[
4
]
=
acos
((
bd
%
ad
));
angles
[
5
]
=
acos
((
ad
%
ab
));
angles
[
6
]
=
acos
(
-
(
bc
%
cd
));
angles
[
7
]
=
acos
((
cd
%
bd
));
angles
[
8
]
=
acos
((
bd
%
bc
));
angles
[
9
]
=
acos
((
ad
%
cd
));
angles
[
10
]
=
acos
(
-
(
cd
%
ac
));
angles
[
11
]
=
acos
((
ac
%
ad
));
angles
[
0
]
=
std
::
acos
(
-
(
ab
%
bc
));
angles
[
1
]
=
std
::
acos
((
bc
%
ac
));
angles
[
2
]
=
std
::
acos
((
ac
%
ab
));
angles
[
3
]
=
std
::
acos
(
-
(
ab
%
bd
));
angles
[
4
]
=
std
::
acos
((
bd
%
ad
));
angles
[
5
]
=
std
::
acos
((
ad
%
ab
));
angles
[
6
]
=
std
::
acos
(
-
(
bc
%
cd
));
angles
[
7
]
=
std
::
acos
((
cd
%
bd
));
angles
[
8
]
=
std
::
acos
((
bd
%
bc
));
angles
[
9
]
=
std
::
acos
((
ad
%
cd
));
angles
[
10
]
=
std
::
acos
(
-
(
cd
%
ac
));
angles
[
11
]
=
std
::
acos
((
ac
%
ad
));
for
(
int
a
=
0
;
a
<
12
;
a
++
)
{
...
...
@@ -178,7 +178,7 @@ static int tet_get_weight(
w2
.
set
(
0.5
,
0.5
*
sqrt3
,
0
);
w3
.
set
(
0.5
,
sqrt3
/
6.0
,
sqrt2
/
sqrt3
);
double
scale
=
pow
(
6.
*
average_tet_volume
/
determinant
(
w1
,
w2
,
w3
),
0.3333333333333
);
double
scale
=
std
::
pow
(
6.
*
average_tet_volume
/
determinant
(
w1
,
w2
,
w3
),
0.3333333333333
);
w1
*=
scale
;
w2
*=
scale
;
...
...
@@ -267,7 +267,7 @@ double tet_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
M2
=
Mab
>
Mcd
?
Mab
:
Mcd
;
M2
=
M2
>
Mef
?
M2
:
Mef
;
const
double
edge_ratio
=
sqrt
(
M2
/
m2
);
const
double
edge_ratio
=
std
::
sqrt
(
M2
/
m2
);
return
fix_range
(
edge_ratio
);
}
...
...
@@ -322,10 +322,10 @@ double tet_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
which_node
=
3
;
}
double
length_product
=
sqrt
(
length_squared
[
which_node
]);
if
(
length_product
<
f
abs
(
jacobi
))
double
length_product
=
std
::
sqrt
(
length_squared
[
which_node
]);
if
(
length_product
<
std
::
abs
(
jacobi
))
{
length_product
=
f
abs
(
jacobi
);
length_product
=
std
::
abs
(
jacobi
);
}
if
(
length_product
<
VERDICT_DBL_MIN
)
...
...
@@ -379,7 +379,7 @@ double tet_radius_ratio(int /*num_nodes*/, const double coordinates[][3])
double
volume
=
tet_volume
(
4
,
coordinates
);
if
(
f
abs
(
volume
)
<
VERDICT_DBL_MIN
)
if
(
std
::
abs
(
volume
)
<
VERDICT_DBL_MIN
)
{
return
(
double
)
VERDICT_DBL_MAX
;
}
...
...
@@ -417,7 +417,7 @@ double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
double
detTet
=
ab
%
(
ac
*
ad
);
if
(
f
abs
(
detTet
)
<
VERDICT_DBL_MIN
)
if
(
std
::
abs
(
detTet
)
<
VERDICT_DBL_MIN
)
{
return
(
double
)
VERDICT_DBL_MAX
;
}
...
...
@@ -442,7 +442,7 @@ double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
double
B
=
ac2
>
ad2
?
ac2
:
ad2
;
double
C
=
bd2
>
cd2
?
bd2
:
cd2
;
double
D
=
A
>
B
?
A
:
B
;
double
hm
=
D
>
C
?
sqrt
(
D
)
:
sqrt
(
C
);
double
hm
=
D
>
C
?
std
::
sqrt
(
D
)
:
std
::
sqrt
(
C
);
bd
=
ab
*
bc
;
A
=
bd
.
length
();
...
...
@@ -453,7 +453,7 @@ double tet_aspect_ratio(int /*num_nodes*/, const double coordinates[][3])
bd
=
bc
*
cd
;
D
=
bd
.
length
();
const
double
aspect_ratio
=
aspect_ratio_normal_coeff
*
hm
*
(
A
+
B
+
C
+
D
)
/
f
abs
(
detTet
);
const
double
aspect_ratio
=
aspect_ratio_normal_coeff
*
hm
*
(
A
+
B
+
C
+
D
)
/
std
::
abs
(
detTet
);
return
fix_range
(
aspect_ratio
);
}
...
...
@@ -486,7 +486,7 @@ double tet_aspect_gamma(int /*num_nodes*/, const double coordinates[][3])
side5
.
set
(
coordinates
[
3
][
0
]
-
coordinates
[
2
][
0
],
coordinates
[
3
][
1
]
-
coordinates
[
2
][
1
],
coordinates
[
3
][
2
]
-
coordinates
[
2
][
2
]);
double
volume
=
f
abs
(
tet_volume
(
4
,
coordinates
));
double
volume
=
std
::
abs
(
tet_volume
(
4
,
coordinates
));
if
(
volume
<
VERDICT_DBL_MIN
)
{
...
...
@@ -494,11 +494,12 @@ double tet_aspect_gamma(int /*num_nodes*/, const double coordinates[][3])
}
else
{
double
srms
=
sqrt
((
side0
.
length_squared
()
+
side1
.
length_squared
()
+
side2
.
length_squared
()
+
side3
.
length_squared
()
+
side4
.
length_squared
()
+
side5
.
length_squared
())
/
6.0
);
double
srms
=
std
::
sqrt
((
side0
.
length_squared
()
+
side1
.
length_squared
()
+
side2
.
length_squared
()
+
side3
.
length_squared
()
+
side4
.
length_squared
()
+
side5
.
length_squared
())
/
6.0
);
double
aspect_ratio_gamma
=
pow
(
srms
,
3
)
/
(
8.48528137423857
*
volume
);
double
aspect_ratio_gamma
=
std
::
pow
(
srms
,
3
)
/
(
8.48528137423857
*
volume
);
return
(
double
)
aspect_ratio_gamma
;
}
}
...
...
@@ -525,7 +526,7 @@ double tet_aspect_frobenius(int /*num_nodes*/, const double coordinates[][3])
double
denominator
=
ab
%
(
ac
*
ad
);
denominator
*=
denominator
;
denominator
*=
2.
;
denominator
=
3.
*
pow
(
denominator
,
one_third
);
denominator
=
3.
*
std
::
pow
(
denominator
,
one_third
);
if
(
denominator
<
VERDICT_DBL_MIN
)
{
...
...
@@ -584,12 +585,12 @@ double tet_minimum_angle(int /*num_nodes*/, const double coordinates[][3])
VerdictVector
bcd
=
bc
*
cd
;
double
nbcd
=
bcd
.
length
();
double
alpha
=
acos
((
abc
%
abd
)
/
(
nabc
*
nabd
));
double
beta
=
acos
((
abc
%
acd
)
/
(
nabc
*
nacd
));
double
gamma
=
acos
((
abc
%
bcd
)
/
(
nabc
*
nbcd
));
double
delta
=
acos
((
abd
%
acd
)
/
(
nabd
*
nacd
));
double
epsilon
=
acos
((
abd
%
bcd
)
/
(
nabd
*
nbcd
));
double
zeta
=
acos
((
acd
%
bcd
)
/
(
nacd
*
nbcd
));
double
alpha
=
std
::
acos
((
abc
%
abd
)
/
(
nabc
*
nabd
));
double
beta
=
std
::
acos
((
abc
%
acd
)
/
(
nabc
*
nacd
));
double
gamma
=
std
::
acos
((
abc
%
bcd
)
/
(
nabc
*
nbcd
));
double
delta
=
std
::
acos
((
abd
%
acd
)
/
(
nabd
*
nacd
));
double
epsilon
=
std
::
acos
((
abd
%
bcd
)
/
(
nabd
*
nbcd
));
double
zeta
=
std
::
acos
((
acd
%
bcd
)
/
(
nacd
*
nbcd
));
alpha
=
alpha
<
beta
?
alpha
:
beta
;
alpha
=
alpha
<
gamma
?
alpha
:
gamma
;
...
...
@@ -717,8 +718,8 @@ double tet_equivolume_skew(int num_nodes, const double coordinates[][3])
double
circumradius
=
num
.
length
()
/
den
;
double
volume
=
tet_volume
(
num_nodes
,
coordinates
);
double
optimal_length
=
circumradius
/
sqrt
(
double
(
3.0
)
/
8.0
);
double
optimal_volume
=
(
1.0
/
12.0
)
*
sqrt
(
double
(
2.0
))
*
pow
(
optimal_length
,
3
);
double
optimal_length
=
circumradius
/
std
::
sqrt
(
double
(
3.0
)
/
8.0
);
double
optimal_volume
=
(
1.0
/
12.0
)
*
std
::
sqrt
(
double
(
2.0
))
*
std
::
pow
(
optimal_length
,
3
);
const
double
eq_v_skew
=
(
optimal_volume
-
volume
)
/
optimal_volume
;
return
fix_range
(
eq_v_skew
);
...
...
@@ -1167,13 +1168,13 @@ double tet_condition(int /*num_nodes*/, const double coordinates[][3])
term2
=
(
c_1
*
c_2
)
%
(
c_1
*
c_2
)
+
(
c_2
*
c_3
)
%
(
c_2
*
c_3
)
+
(
c_1
*
c_3
)
%
(
c_1
*
c_3
);
det
=
c_1
%
(
c_2
*
c_3
);
if
(
f
abs
(
det
)
<=
VERDICT_DBL_MIN
)
if
(
std
::
abs
(
det
)
<=
VERDICT_DBL_MIN
)
{
return
VERDICT_DBL_MAX
;
}
else
{
condition
=
sqrt
(
term1
*
term2
)
/
(
3.0
*
det
);
condition
=
std
::
sqrt
(
term1
*
term2
)
/
(
3.0
*
det
);
}
return
(
double
)
condition
;
...
...
@@ -1256,7 +1257,7 @@ double tet_shape(int /*num_nodes*/, const double coordinates[][3])
{
return
(
double
)
0.0
;
}
double
num
=
3
*
pow
(
sqrt2
*
jacobian
,
two_thirds
);
double
num
=
3
*
std
::
pow
(
sqrt2
*
jacobian
,
two_thirds
);
double
den
=
1.5
*
(
edge0
%
edge0
+
edge2
%
edge2
+
edge3
%
edge3
)
-
(
edge0
%
-
edge2
+
-
edge2
%
edge3
+
edge3
%
edge0
);
...
...
@@ -1423,7 +1424,7 @@ double tet_distortion(int num_nodes, const double coordinates[][3])
minimum_jacobian
=
jacobian
;
}
}
if
(
f
abs
(
element_volume
)
>
0.0
)
if
(
std
::
abs
(
element_volume
)
>
0.0
)
{
distortion
=
minimum_jacobian
/
element_volume
;
}
...
...
@@ -1511,7 +1512,7 @@ double tet_timestep(int num_nodes, const double coordinates[][3], double density
double
M
=
youngs_modulus
*
(
1
-
poissons_ratio
)
/
((
1
-
2
*
poissons_ratio
)
*
(
1
+
poissons_ratio
));
double
denominator
=
sqrt
(
M
/
density
);
double
denominator
=
std
::
sqrt
(
M
/
density
);
return
char_length
/
denominator
;
}
...
...
@@ -1596,7 +1597,7 @@ double calculate_tet4_outer_radius(const double coordinates[][3])
double
BC
=
(
nE
[
3
]
-
nE
[
1
]).
length
();
double
CC
=
(
nE
[
2
]
-
nE
[
1
]).
length
();
double
VP
=
((
nE
[
1
]
-
nE
[
0
])
*
(
nE
[
2
]
-
nE
[
0
]))
%
(
nE
[
3
]
-
nE
[
0
])
/
6
;
double
outer_radius
=
sqrt
((
aC
*
AC
+
bC
*
BC
+
cC
*
CC
)
*
(
aC
*
AC
+
bC
*
BC
-
cC
*
CC
)
*
double
outer_radius
=
std
::
sqrt
((
aC
*
AC
+
bC
*
BC
+
cC
*
CC
)
*
(
aC
*
AC
+
bC
*
BC
-
cC
*
CC
)
*
(
aC
*
AC
-
bC
*
BC
+
cC
*
CC
)
*
(
-
aC
*
AC
+
bC
*
BC
+
cC
*
CC
))
/
24
/
VP
;
...
...
This diff is collapsed.
Click to expand it.
V_TriMetric.cpp
View file @
9d11ce5b
...
...
@@ -27,7 +27,7 @@
namespace
VERDICT_NAMESPACE
{
static
const
double
sqrt3
=
sqrt
(
3.0
);
static
const
double
sqrt3
=
std
::
sqrt
(
3.0
);
static
const
double
aspect_ratio_normal_coeff
=
sqrt3
/
6.
;
static
const
double
two_times_sqrt3
=
2
*
sqrt3
;
static
const
double
two_over_sqrt3
=
2.
/
sqrt3
;
...
...
@@ -43,7 +43,7 @@ static int tri_get_weight(
m21
=
0
;
m12
=
0.5
;
m22
=
0.5
*
sqrt3
;
double
scale
=
sqrt
(
2.0
*
average_tri_area
/
(
m11
*
m22
-
m21
*
m12
));
double
scale
=
std
::
sqrt
(
2.0
*
average_tri_area
/
(
m11
*
m22
-
m21
*
m12
));
m11
*=
scale
;
m21
*=
scale
;
...
...
@@ -128,7 +128,7 @@ double tri_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
else
{
double
edge_ratio
;
edge_ratio
=
sqrt
(
M2
/
m2
);
edge_ratio
=
std
::
sqrt
(
M2
/
m2
);
if
(
edge_ratio
>
0
)
{
...
...
@@ -589,7 +589,7 @@ double tri_relative_size_squared(
return
0.0
;
}
double
size
=
pow
(
deta
/
detw
,
2
);
double
size
=
std
::
pow
(
deta
/
detw
,
2
);
double
rel_size
=
std
::
min
(
size
,
1.0
/
size
);
...
...
@@ -718,7 +718,7 @@ double tri_distortion(int num_nodes, const double coordinates[][3])
for
(
ja
=
0
;
ja
<
num_nodes
;
ja
++
)
{
dot_product
=
normal_at_nodes
[
0
]
%
normal_at_nodes
[
ja
];
if
(
f
abs
(
dot_product
)
<
0.99
)
if
(
std
::
abs
(
dot_product
)
<
0.99
)
{
flat_element
=
false
;
break
;
...
...
@@ -729,7 +729,7 @@ double tri_distortion(int num_nodes, const double coordinates[][3])
double
thickness
,
thickness_gauss
;
double
distrt
;
// get_tri_thickness(tri, element_area, thickness );
thickness
=
0.001
*
sqrt
(
element_area
);
thickness
=
0.001
*
std
::
sqrt
(
element_area
);
// set thickness gauss point location
double
zl
=
0.5773502691896
;
...
...
This diff is collapsed.
Click to expand it.
V_WedgeMetric.cpp
View file @
9d11ce5b
...
...
@@ -506,7 +506,7 @@ double wedge_edge_ratio(int /*num_nodes*/, const double coordinates[][3])
min
=
i2
;
}
double
edge_ratio
=
sqrt
(
max
/
min
);
double
edge_ratio
=
std
::
sqrt
(
max
/
min
);
if
(
std
::
isnan
(
edge_ratio
))
{
...
...
@@ -886,7 +886,7 @@ double wedge_distortion(int num_nodes, const double coordinates[][3])
double
master_volume
=
0.433013
;
double
current_volume
=
wedge_volume
(
num_nodes
,
coordinates
);
double
distortion
=
VERDICT_DBL_MAX
;
if
(
f
abs
(
current_volume
)
>
0.0
)
if
(
std
::
abs
(
current_volume
)
>
0.0
)
distortion
=
jacobian
*
master_volume
/
current_volume
/
0.866025
;
if
(
std
::
isnan
(
distortion
))
...
...
@@ -1022,7 +1022,7 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
2
][
0
]
-
coordinates
[
0
][
0
],
coordinates
[
2
][
1
]
-
coordinates
[
0
][
1
],
coordinates
[
2
][
2
]
-
coordinates
[
0
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
(
vec2
%
(
vec1
*
vec3
));
min_jacobian
=
current_jacobian
/
lengths
;
...
...
@@ -1037,7 +1037,7 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
0
][
0
]
-
coordinates
[
1
][
0
],
coordinates
[
0
][
1
]
-
coordinates
[
1
][
1
],
coordinates
[
0
][
2
]
-
coordinates
[
1
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
min_jacobian
=
std
::
min
(
current_jacobian
/
lengths
,
min_jacobian
);
...
...
@@ -1052,7 +1052,7 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
1
][
0
]
-
coordinates
[
2
][
0
],
coordinates
[
1
][
1
]
-
coordinates
[
2
][
1
],
coordinates
[
1
][
2
]
-
coordinates
[
2
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
min_jacobian
=
std
::
min
(
current_jacobian
/
lengths
,
min_jacobian
);
...
...
@@ -1067,7 +1067,7 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
5
][
0
]
-
coordinates
[
3
][
0
],
coordinates
[
5
][
1
]
-
coordinates
[
3
][
1
],
coordinates
[
5
][
2
]
-
coordinates
[
3
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
min_jacobian
=
std
::
min
(
current_jacobian
/
lengths
,
min_jacobian
);
...
...
@@ -1082,7 +1082,7 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
3
][
0
]
-
coordinates
[
4
][
0
],
coordinates
[
3
][
1
]
-
coordinates
[
4
][
1
],
coordinates
[
3
][
2
]
-
coordinates
[
4
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
min_jacobian
=
std
::
min
(
current_jacobian
/
lengths
,
min_jacobian
);
...
...
@@ -1097,12 +1097,12 @@ double wedge_scaled_jacobian(int /*num_nodes*/, const double coordinates[][3])
vec3
.
set
(
coordinates
[
2
][
0
]
-
coordinates
[
5
][
0
],
coordinates
[
2
][
1
]
-
coordinates
[
5
][
1
],
coordinates
[
2
][
2
]
-
coordinates
[
5
][
2
]);
lengths
=
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
lengths
=
std
::
sqrt
(
vec1
.
length_squared
()
*
vec2
.
length_squared
()
*
vec3
.
length_squared
());
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
min_jacobian
=
std
::
min
(
current_jacobian
/
lengths
,
min_jacobian
);
min_jacobian
*=
2
/
sqrt
(
3.0
);
min_jacobian
*=
2
/
std
::
sqrt
(
3.0
);
if
(
min_jacobian
>
0
)
{
...
...
@@ -1145,8 +1145,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
@@ -1168,8 +1168,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
@@ -1191,8 +1191,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
@@ -1214,8 +1214,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
@@ -1237,8 +1237,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
@@ -1260,8 +1260,8 @@ double wedge_shape(int /*num_nodes*/, const double coordinates[][3])
current_jacobian
=
vec2
%
(
vec1
*
vec3
);
if
(
current_jacobian
>
VERDICT_DBL_MIN
)
{
norm_jacobi
=
current_jacobian
*
2.0
/
sqrt
(
3.0
);
current_shape
=
3
*
pow
(
norm_jacobi
,
two_thirds
)
/
norm_jacobi
=
current_jacobian
*
2.0
/
std
::
sqrt
(
3.0
);
current_shape
=
3
*
std
::
pow
(
norm_jacobi
,
two_thirds
)
/
(
vec1
.
length_squared
()
+
vec2
.
length_squared
()
+
vec3
.
length_squared
());
min_shape
=
std
::
min
(
current_shape
,
min_shape
);
}
...
...
This diff is collapsed.
Click to expand it.
VerdictVector.cpp
View file @
9d11ce5b
...
...
@@ -21,7 +21,6 @@
#include "VerdictVector.hpp"
#include "verdict.h"
#include <cfloat>
#include <cmath>
namespace
VERDICT_NAMESPACE
...
...
@@ -54,16 +53,16 @@ double VerdictVector::interior_angle(const VerdictVector& otherVector)
if
((
cosAngle
>
1.0
)
&&
(
cosAngle
<
1.0001
))
{
cosAngle
=
1.0
;
angleRad
=
acos
(
cosAngle
);
angleRad
=
std
::
acos
(
cosAngle
);
}
else
if
(
cosAngle
<
-
1.0
&&
cosAngle
>
-
1.0001
)
{
cosAngle
=
-
1.0
;
angleRad
=
acos
(
cosAngle
);
angleRad
=
std
::
acos
(
cosAngle
);
}
else
if
(
cosAngle
>=
-
1.0
&&
cosAngle
<=
1.0
)
{
angleRad
=
acos
(
cosAngle
);
angleRad
=
std
::
acos
(
cosAngle
);
}
else
{
...
...
This diff is collapsed.
Click to expand it.
VerdictVector.hpp
View file @
9d11ce5b
...
...
@@ -347,7 +347,7 @@ inline VerdictVector& VerdictVector::operator/=(const double scalar)
// Returns the normalized 'this'.
inline
VerdictVector
operator
~
(
const
VerdictVector
&
vec
)
{
double
mag
=
sqrt
(
vec
.
xVal
*
vec
.
xVal
+
vec
.
yVal
*
vec
.
yVal
+
vec
.
zVal
*
vec
.
zVal
);
double
mag
=
std
::
sqrt
(
vec
.
xVal
*
vec
.
xVal
+
vec
.
yVal
*
vec
.
yVal
+
vec
.
zVal
*
vec
.
zVal
);
VerdictVector
temp
=
vec
;
if
(
mag
!=
0.0
)
...
...
@@ -423,7 +423,7 @@ inline double VerdictVector::length_squared() const
inline
double
VerdictVector
::
length
()
const
{
return
(
sqrt
(
xVal
*
xVal
+
yVal
*
yVal
+
zVal
*
zVal
));
return
(
std
::
sqrt
(
xVal
*
xVal
+
yVal
*
yVal
+
zVal
*
zVal
));
}
inline
double
VerdictVector
::
normalize
()
...
...
This diff is collapsed.
Click to expand it.
unittests/verdict.test.cpp
View file @
9d11ce5b
...
...
@@ -19,15 +19,12 @@
*/
#include "verdict.h"
#include "V_HexMetric.hpp"
#include "gtest/gtest.h"
#include <cmath>
#include <fstream>
#include <functional>
#include <iostream>
#include <stdexcept>
#include <vector>
#define MAX_NODES_PER_ELEMENT 27
...
...
@@ -71,7 +68,7 @@ void runtest(test_case& this_testcase)
<<
std
::
endl
;
EXPECT_NEAR
(
calculated_answer
,
expected_answer
,
f
abs
(
expected_answer
)
*
VERDICT_RELATIVE_TOL
+
VERDICT_ABSOLUTE_TOL
);
std
::
abs
(
expected_answer
)
*
VERDICT_RELATIVE_TOL
+
VERDICT_ABSOLUTE_TOL
);
// old test using strings
/*
...
...
@@ -88,7 +85,7 @@ void runtest(test_case& this_testcase)
// for expected values of zero, we only care about absolute tolerance, not relative tolerance,
std::string expected_v_str = expected.str();
std::string answer_v_str = answer.str();
if (
f
abs(answer_from_lib) < VERDICT_ABSOLUTE_TOL )
if (
std::
abs(answer_from_lib) < VERDICT_ABSOLUTE_TOL )
{
std::stringstream expected_abs;
expected_abs.setf(std::ios::scientific, std::ios::floatfield);
...
...
@@ -128,8 +125,9 @@ void runtest(test_case& this_testcase)
TEST
(
verdict
,
tet_incircle_right
)
{
test_case
testcase
=
{
"tet_incircle_right"
,
{
{
verdict
::
tet_inradius
,
0.5
-
1.
/
sqrt
(
12
)
}
},
4
,
{
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
1
}
}
};
test_case
testcase
=
{
"tet_incircle_right"
,
{
{
verdict
::
tet_inradius
,
0.5
-
1.
/
std
::
sqrt
(
12
)
}
},
4
,
{
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
1
}
}
};
runtest
(
testcase
);
}
...
...
@@ -137,7 +135,7 @@ TEST(verdict, tet_incircle_right)
TEST
(
verdict
,
tet_incircle_right2
)
{
test_case
testcase
=
{
"tet_incircle_right2"
,
{
{
verdict
::
tet_inradius
,
2.0
*
(
0.5
-
1.
/
sqrt
(
12
))
}
},
4
,
{
{
verdict
::
tet_inradius
,
2.0
*
(
0.5
-
1.
/
std
::
sqrt
(
12
))
}
},
4
,
{
{
0
,
0
,
0
},
{
2
,
0
,
0
},
{
0
,
2
,
0
},
{
0
,
0
,
2
}
}
};
runtest
(
testcase
);
...
...
@@ -148,29 +146,29 @@ TEST(verdict, tet_incircle_equilateral)
const
double
pi
=
verdict
::
VERDICT_PI
;
// equilateral tet, with side length 1
double
l2
=
sin
(
30
*
pi
/
180
)
/
sin
(
120
*
pi
/
180
);
double
h
=
sqrt
(
1
-
l2
*
l2
);
double
l2
=
std
::
sin
(
30
*
pi
/
180
)
/
std
::
sin
(
120
*
pi
/
180
);
double
h
=
std
::
sqrt
(
1
-
l2
*
l2
);
test_case
testcase
=
{
"tet_incircle_equilateral"
,
{
{
verdict
::
tet_inradius
,
sqrt
(
6.
)
/
12.
}
},
4
,
{
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
cos
(
60
*
pi
/
180
),
sin
(
60
*
pi
/
180
),
0
},
{
l2
*
cos
(
30
*
pi
/
180
),
l2
*
sin
(
30
*
pi
/
180
),
h
}
}
};
test_case
testcase
=
{
"tet_incircle_equilateral"
,
{
{
verdict
::
tet_inradius
,
std
::
sqrt
(
6.
)
/
12.
}
},
4
,
{
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
std
::
cos
(
60
*
pi
/
180
),
std
::
sin
(
60
*
pi
/
180
),
0
},
{
l2
*
std
::
cos
(
30
*
pi
/
180
),
l2
*
std
::
sin
(
30
*
pi
/
180
),
h
}
}
};
runtest
(
testcase
);
}
TEST
(
verdict
,
tet_incircle_equilateral2
)
{
const
double
pi
=
2.
*
asin
(
1.
);
const
double
pi
=
2.
*
std
::
asin
(
1.
);
// equilateral tet, with side length 2
double
l2
=
2.
*
sin
(
30
*
pi
/
180
)
/
sin
(
120
*
pi
/
180
);
double
h
=
sqrt
(
4.
-
l2
*
l2
);
double
l2
=
2.
*
std
::
sin
(
30
*
pi
/
180
)
/
std
::
sin
(
120
*
pi
/
180
);
double
h
=
std
::
sqrt
(
4.
-
l2
*
l2
);
test_case
testcase
=
{
"tet_incircle_equilateral2"
,
{
{
verdict
::
tet_inradius
,
2
*
sqrt
(
6.
)
/
12.
}
},
4
,
{
{
0
,
0
,
0
},
{
2
,
0
,
0
},
{
2
*
cos
(
60
*
pi
/
180
),
2
*
sin
(
60
*
pi
/
180
),
0
},
{
l2
*
cos
(
30
*
pi
/
180
),
l2
*
sin
(
30
*
pi
/
180
),
h
}
}
};
{
{
verdict
::
tet_inradius
,
2
*
std
::
sqrt
(
6.
)
/
12.
}
},
4
,
{
{
0
,
0
,
0
},
{
2
,
0
,
0
},
{
2
*
std
::
cos
(
60
*
pi
/
180
),
2
*
std
::
sin
(
60
*
pi
/
180
),
0
},
{
l2
*
std
::
cos
(
30
*
pi
/
180
),
l2
*
std
::
sin
(
30
*
pi
/
180
),
h
}
}
};
runtest
(
testcase
);
}
...
...
@@ -186,7 +184,7 @@ TEST(verdict, tet_incircle_flat)
TEST
(
verdict
,
tet_incircle_inverted
)
{
test_case
testcase
=
{
"tet_incircle_inverted"
,
{
{
verdict
::
tet_inradius
,
-
0.5
+
1.
/
sqrt
(
12
)
}
},
4
,
{
{
verdict
::
tet_inradius
,
-
0.5
+
1.
/
std
::
sqrt
(
12
)
}
},
4
,
{
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
0
,
1
,
0
},
{
0
,
0
,
-
1
}
}
};
runtest
(
testcase
);
...
...
@@ -271,7 +269,7 @@ TEST(verdict, wedge_simple)
{
test_case
testcase
=
{
"wedge_simple"
,
{
{
verdict
::
wedge_volume
,
0.5
},
{
verdict
::
wedge_equiangle_skew
,
0.25
},
{
verdict
::
wedge_edge_ratio
,
sqrt
(
2.
)
},
{
verdict
::
wedge_edge_ratio
,
std
::
sqrt
(
2.
)
},
{
verdict
::
wedge_max_aspect_frobenius
,
1.1357042248
},
{
verdict
::
wedge_mean_aspect_frobenius
,
1.0978474174
},
{
verdict
::
wedge_distortion
,
1
},
{
verdict
::
wedge_max_stretch
,
1
}
},
...
...
@@ -489,7 +487,7 @@ TEST(verdict, tet_test1)
{
{
verdict
::
tet_volume
,
166.66666667
},
{
verdict
::
tet_condition
,
1.2247448714
},
{
verdict
::
tet_jacobian
,
1000
},
{
verdict
::
tet_shape
,
0.83994736660
},
{
verdict
::
tet_distortion
,
1
},
/* 6 */
{
verdict
::
tet_edge_ratio
,
sqrt
(
2.
)
},
/* 6 */
{
verdict
::
tet_edge_ratio
,
std
::
sqrt
(
2.
)
},
/* 7 */
{
verdict
::
tet_radius_ratio
,
1.3660254038
},
/* 8 */
{
verdict
::
tet_aspect_ratio
,
1.3660254038
},
/* 9 */
{
verdict
::
tet_aspect_frobenius
,
1.1905507890
},
...
...
@@ -608,7 +606,7 @@ TEST(verdict, hex_test2_perfect_cube)
test_case
testcase
=
{
"hex_test2_perfect_cube"
,
{
{
verdict
::
hex_skew
,
0.
},
{
verdict
::
hex_taper
,
0.
},
{
verdict
::
hex_volume
,
1.0
},
{
verdict
::
hex_stretch
,
1.0
},
{
verdict
::
hex_diagonal
,
1.0
},
{
verdict
::
hex_dimension
,
1.
/
sqrt
(
3.
)
},
{
verdict
::
hex_condition
,
1.0
},
{
verdict
::
hex_dimension
,
1.
/
std
::
sqrt
(
3.
)
},
{
verdict
::
hex_condition
,
1.0
},
{
verdict
::
hex_med_aspect_frobenius
,
1
},
{
verdict
::
hex_jacobian
,
1.0
},
{
verdict
::
hex_shear
,
1.0
},
{
verdict
::
hex_shape
,
1.0
},
{
verdict
::
hex_distortion
,
1.0
},
{
verdict
::
hex_nodal_jacobian_ratio
,
1.0
},
{
verdict
::
hex_oddy
,
0.
},
...
...
@@ -633,7 +631,7 @@ TEST(verdict, hex_test3_flat)
{
verdict
::
hex_distortion
,
verdict
::
VERDICT_DBL_MAX
},
{
verdict
::
hex_nodal_jacobian_ratio
,
-
verdict
::
VERDICT_DBL_MAX
},
{
verdict
::
hex_oddy
,
verdict
::
VERDICT_DBL_MAX
},
{
verdict
::
hex_edge_ratio
,
1.0
/
sqrt
(
0.02
)
},
{
verdict
::
hex_equiangle_skew
,
0.5
}
},
{
verdict
::
hex_edge_ratio
,
1.0
/
std
::
sqrt
(
0.02
)
},
{
verdict
::
hex_equiangle_skew
,
0.5
}
},
8
,
{
// squashed flat
{
0
,
0
,
0
},
{
1
,
0
,
0
},
{
1
,
1
,
0
},
{
0
,
1
,
0
},
{
0.1
,
0.1
,
0
},
{
0.9
,
0.1
,
0
},
...
...
@@ -651,7 +649,7 @@ TEST(verdict, hex_test4_inside_out)
/* 3 */
{
verdict
::
hex_volume
,
-
1.
},
// !=
/* 4 */
{
verdict
::
hex_stretch
,
1
},
/* 5 */
{
verdict
::
hex_diagonal
,
1.0
},
/* 6 */
{
verdict
::
hex_dimension
,
1.
/
sqrt
(
3.
)
},
/* 6 */
{
verdict
::
hex_dimension
,
1.
/
std
::
sqrt
(
3.
)
},
/* 7 */
{
verdict
::
hex_condition
,
verdict
::
VERDICT_DBL_MAX
},
// !=
/* 8 */
{
verdict
::
hex_med_aspect_frobenius
,
verdict
::
VERDICT_DBL_MAX
},
// !=
/* 9 */
{
verdict
::
hex_jacobian
,
-
1.
},
// !=
...
...
This diff is collapsed.
Click to expand it.
verdict_defines.hpp
View file @
9d11ce5b
...
...
@@ -43,7 +43,7 @@ inline double determinant(VerdictVector v1, VerdictVector v2, VerdictVector v3)
return
VerdictVector
::
Dot
(
v1
,
(
v2
*
v3
));
}
static
const
double
sqrt_2
=
sqrt
(
2.0
);
static
const
double
sqrt_2
=
std
::
sqrt
(
2.0
);
inline
double
normalize_jacobian
(
double
jacobi
,
VerdictVector
&
v1
,
VerdictVector
&
v2
,
VerdictVector
&
v3
,
int
tet_flag
=
0
)
...
...
@@ -60,13 +60,13 @@ inline double normalize_jacobian(
l1
=
v1
.
length_squared
();
l2
=
v2
.
length_squared
();
l3
=
v3
.
length_squared
();
length_product
=
sqrt
(
l1
*
l2
*
l3
);
length_product
=
std
::
sqrt
(
l1
*
l2
*
l3
);
// if some numerical scaling problem, or just plain roundoff,
// then push back into range [-1,1].
if
(
length_product
<
f
abs
(
jacobi
))
if
(
length_product
<
std
::
abs
(
jacobi
))
{
length_product
=
f
abs
(
jacobi
);
length_product
=
std
::
abs
(
jacobi
);
}
if
(
tet_flag
==
1
)
...
...
@@ -89,7 +89,7 @@ inline double norm_squared(double m11, double m21, double m12, double m22)
inline
int
skew_matrix
(
double
gm11
,
double
gm12
,
double
gm22
,
double
det
,
double
&
qm11
,
double
&
qm21
,
double
&
qm12
,
double
&
qm22
)
{
double
tmp
=
sqrt
(
gm11
*
gm22
);
double
tmp
=
std
::
sqrt
(
gm11
*
gm22
);
if
(
tmp
==
0
)
{
return
false
;
...
...
@@ -133,14 +133,14 @@ inline void form_Q(const VerdictVector& v1, const VerdictVector& v2, const Verdi
g23
=
VerdictVector
::
Dot
(
v2
,
v3
);
g33
=
VerdictVector
::
Dot
(
v3
,
v3
);
double
rtg11
=
sqrt
(
g11
);
double
rtg22
=
sqrt
(
g22
);
double
rtg33
=
sqrt
(
g33
);
double
rtg11
=
std
::
sqrt
(
g11
);
double
rtg22
=
std
::
sqrt
(
g22
);
double
rtg33
=
std
::
sqrt
(
g33
);
VerdictVector
temp1
;
temp1
=
v1
*
v2
;
double
cross
=
sqrt
(
VerdictVector
::
Dot
(
temp1
,
temp1
));
double
cross
=
std
::
sqrt
(
VerdictVector
::
Dot
(
temp1
,
temp1
));
double
q11
,
q21
,
q31
;
double
q12
,
q22
,
q32
;
...
...
@@ -196,7 +196,7 @@ inline double skew_x(VerdictVector& q1, VerdictVector& q2, VerdictVector& q3, Ve
inverse
(
x1
,
x2
,
x3
,
u1
,
u2
,
u3
);
normsq1
=
norm_squared
(
x1
,
x2
,
x3
);
normsq2
=
norm_squared
(
u1
,
u2
,
u3
);
kappa
=
sqrt
(
normsq1
*
normsq2
);
kappa
=
std
::
sqrt
(
normsq1
*
normsq2
);
double
skew
=
0
;
if
(
kappa
>
VERDICT_DBL_MIN
)
...
...
This diff is collapsed.
Click to expand it.
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment
Menu
Projects
Groups
Snippets
Help