feat: add mitc4 geometry directors
This commit is contained in:
NINI
+220
-8
@@ -198,6 +198,22 @@ inline Real norm(const Vec3& value) {
|
||||
return std::sqrt(dot(value, value));
|
||||
}
|
||||
|
||||
inline bool isFinite(Real value) {
|
||||
return std::isfinite(value);
|
||||
}
|
||||
|
||||
inline bool isFinite(const Vec3& value) {
|
||||
return isFinite(value.x) && isFinite(value.y) && isFinite(value.z);
|
||||
}
|
||||
|
||||
inline std::optional<Vec3> normalizedIfValid(const Vec3& value, Real tolerance = 1.0e-12) {
|
||||
const Real length = norm(value);
|
||||
if (!isFinite(length) || length <= tolerance) {
|
||||
return std::nullopt;
|
||||
}
|
||||
return (1.0 / length) * value;
|
||||
}
|
||||
|
||||
inline Vec3 normalized(const Vec3& value) {
|
||||
const Real length = norm(value);
|
||||
if (length <= std::numeric_limits<Real>::epsilon()) {
|
||||
@@ -1222,15 +1238,211 @@ struct LocalBasis {
|
||||
Vec3 e3;
|
||||
};
|
||||
|
||||
struct MITC4NaturalPoint {
|
||||
Real xi = 0.0;
|
||||
Real eta = 0.0;
|
||||
};
|
||||
|
||||
struct MITC4TyingPoint {
|
||||
std::string label;
|
||||
MITC4NaturalPoint natural;
|
||||
std::array<LocalIndex, 2> edge_node_indices{};
|
||||
};
|
||||
|
||||
inline std::array<MITC4NaturalPoint, 4> mitc4NodeNaturalCoordinates() {
|
||||
return {{{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}}};
|
||||
}
|
||||
|
||||
inline std::array<MITC4TyingPoint, 4> mitc4TyingPoints() {
|
||||
return {{{"A", {0.0, -1.0}, {0, 1}},
|
||||
{"B", {-1.0, 0.0}, {0, 3}},
|
||||
{"C", {0.0, 1.0}, {3, 2}},
|
||||
{"D", {1.0, 0.0}, {1, 2}}}};
|
||||
}
|
||||
|
||||
struct MITC4DirectorFrame {
|
||||
Vec3 v1;
|
||||
Vec3 v2;
|
||||
Vec3 vn;
|
||||
};
|
||||
|
||||
struct MITC4MidsurfaceDerivatives {
|
||||
ShapeData shape;
|
||||
Vec3 g1;
|
||||
Vec3 g2;
|
||||
};
|
||||
|
||||
struct MITC4Geometry {
|
||||
std::array<Vec3, 4> coordinates{};
|
||||
Real thickness = 0.0;
|
||||
ShapeData center_shape;
|
||||
Vec3 g1_center;
|
||||
Vec3 g2_center;
|
||||
Vec3 center_normal;
|
||||
std::array<MITC4DirectorFrame, 4> nodal_frames{};
|
||||
std::vector<Diagnostic> diagnostics;
|
||||
|
||||
bool ok() const {
|
||||
return !hasError(diagnostics);
|
||||
}
|
||||
};
|
||||
|
||||
struct MITC4IntegrationBasis {
|
||||
ShapeData shape;
|
||||
Vec3 g1;
|
||||
Vec3 g2;
|
||||
Vec3 g3;
|
||||
Real jacobian = 0.0;
|
||||
LocalBasis local;
|
||||
std::vector<Diagnostic> diagnostics;
|
||||
|
||||
bool ok() const {
|
||||
return !hasError(diagnostics);
|
||||
}
|
||||
};
|
||||
|
||||
inline Vec3 globalEX() {
|
||||
return {1.0, 0.0, 0.0};
|
||||
}
|
||||
|
||||
inline Vec3 globalEY() {
|
||||
return {0.0, 1.0, 0.0};
|
||||
}
|
||||
|
||||
inline Vec3 globalEZ() {
|
||||
return {0.0, 0.0, 1.0};
|
||||
}
|
||||
|
||||
inline Diagnostic mitc4Diagnostic(std::string code, std::string message) {
|
||||
return makeDiagnostic(Severity::Error, std::move(code), std::move(message), "mitc4", "<element>", 0);
|
||||
}
|
||||
|
||||
inline MITC4MidsurfaceDerivatives mitc4MidsurfaceDerivatives(const std::array<Vec3, 4>& coordinates, Real xi, Real eta) {
|
||||
MITC4MidsurfaceDerivatives result;
|
||||
result.shape = shapeFunctions(xi, eta);
|
||||
for (std::size_t i = 0; i < 4; ++i) {
|
||||
result.g1 = result.g1 + result.shape.dr[i] * coordinates[i];
|
||||
result.g2 = result.g2 + result.shape.ds[i] * coordinates[i];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
inline std::optional<Vec3> firstNormalizedCross(const std::array<Vec3, 3>& axes, const Vec3& vector, Real tolerance) {
|
||||
for (const Vec3& axis : axes) {
|
||||
auto candidate = normalizedIfValid(cross(axis, vector), tolerance);
|
||||
if (candidate) {
|
||||
return candidate;
|
||||
}
|
||||
}
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
inline std::optional<MITC4DirectorFrame> buildMITC4DirectorFrame(const Vec3& normal, Real tolerance) {
|
||||
auto v1 = normalizedIfValid(cross(globalEY(), normal), tolerance);
|
||||
if (!v1) {
|
||||
v1 = firstNormalizedCross({globalEZ(), globalEX(), globalEY()}, normal, tolerance);
|
||||
}
|
||||
if (!v1) {
|
||||
return std::nullopt;
|
||||
}
|
||||
auto v2 = normalizedIfValid(cross(normal, *v1), tolerance);
|
||||
if (!v2) {
|
||||
return std::nullopt;
|
||||
}
|
||||
return MITC4DirectorFrame{*v1, *v2, normal};
|
||||
}
|
||||
|
||||
inline MITC4Geometry buildMITC4Geometry(const std::array<Vec3, 4>& coordinates, Real thickness, Real tolerance = 1.0e-12) {
|
||||
MITC4Geometry geometry;
|
||||
geometry.coordinates = coordinates;
|
||||
geometry.thickness = thickness;
|
||||
if (!isFinite(thickness) || thickness <= tolerance) {
|
||||
geometry.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-THICKNESS", "MITC4 shell thickness must be positive and finite"));
|
||||
}
|
||||
for (const Vec3& coordinate : coordinates) {
|
||||
if (!isFinite(coordinate)) {
|
||||
geometry.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-COORDINATE", "MITC4 element coordinates must be finite"));
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
const auto center = mitc4MidsurfaceDerivatives(coordinates, 0.0, 0.0);
|
||||
geometry.center_shape = center.shape;
|
||||
geometry.g1_center = center.g1;
|
||||
geometry.g2_center = center.g2;
|
||||
const auto normal = normalizedIfValid(cross(center.g1, center.g2), tolerance);
|
||||
if (!normal) {
|
||||
geometry.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-NORMAL", "MITC4 element center normal is near zero"));
|
||||
return geometry;
|
||||
}
|
||||
geometry.center_normal = *normal;
|
||||
|
||||
const auto frame = buildMITC4DirectorFrame(*normal, tolerance);
|
||||
if (!frame) {
|
||||
geometry.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 nodal director basis could not be constructed"));
|
||||
return geometry;
|
||||
}
|
||||
geometry.nodal_frames.fill(*frame);
|
||||
return geometry;
|
||||
}
|
||||
|
||||
inline MITC4IntegrationBasis computeMITC4IntegrationBasis(const MITC4Geometry& geometry, Real xi, Real eta, Real zeta,
|
||||
Real tolerance = 1.0e-12) {
|
||||
MITC4IntegrationBasis result;
|
||||
result.diagnostics = geometry.diagnostics;
|
||||
result.shape = shapeFunctions(xi, eta);
|
||||
for (std::size_t i = 0; i < 4; ++i) {
|
||||
const Vec3& coordinate = geometry.coordinates[i];
|
||||
const Vec3& normal = geometry.nodal_frames[i].vn;
|
||||
result.g1 = result.g1 + result.shape.dr[i] * coordinate +
|
||||
(0.5 * zeta * geometry.thickness * result.shape.dr[i]) * normal;
|
||||
result.g2 = result.g2 + result.shape.ds[i] * coordinate +
|
||||
(0.5 * zeta * geometry.thickness * result.shape.ds[i]) * normal;
|
||||
result.g3 = result.g3 + (0.5 * geometry.thickness * result.shape.n[i]) * normal;
|
||||
}
|
||||
|
||||
result.jacobian = dot(cross(result.g1, result.g2), result.g3);
|
||||
if (!isFinite(result.jacobian) || std::fabs(result.jacobian) <= tolerance) {
|
||||
result.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 element Jacobian is near zero"));
|
||||
}
|
||||
|
||||
const auto e3 = normalizedIfValid(result.g3, tolerance);
|
||||
if (!e3) {
|
||||
result.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis normal is near zero"));
|
||||
return result;
|
||||
}
|
||||
auto e1 = normalizedIfValid(cross(result.g2, *e3), tolerance);
|
||||
if (!e1) {
|
||||
e1 = firstNormalizedCross({globalEY(), globalEZ(), globalEX()}, *e3, tolerance);
|
||||
}
|
||||
if (!e1) {
|
||||
result.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis tangent could not be constructed"));
|
||||
return result;
|
||||
}
|
||||
const auto e2 = normalizedIfValid(cross(*e3, *e1), tolerance);
|
||||
if (!e2) {
|
||||
result.diagnostics.push_back(
|
||||
mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis is not right-handed"));
|
||||
return result;
|
||||
}
|
||||
result.local = {*e1, *e2, *e3};
|
||||
return result;
|
||||
}
|
||||
|
||||
inline LocalBasis computeLocalBasis(const std::array<Vec3, 4>& coordinates) {
|
||||
Vec3 v1 = 0.5 * ((coordinates[1] - coordinates[0]) + (coordinates[2] - coordinates[3]));
|
||||
Vec3 v2 = 0.5 * ((coordinates[3] - coordinates[0]) + (coordinates[2] - coordinates[1]));
|
||||
Vec3 e1 = normalized(v1);
|
||||
v2 = v2 - dot(v2, e1) * e1;
|
||||
Vec3 e2 = normalized(v2);
|
||||
Vec3 e3 = normalized(cross(e1, e2));
|
||||
e2 = normalized(cross(e3, e1));
|
||||
return {e1, e2, e3};
|
||||
const MITC4Geometry geometry = buildMITC4Geometry(coordinates, 1.0);
|
||||
if (!geometry.ok()) {
|
||||
throw std::runtime_error("invalid MITC4 geometry");
|
||||
}
|
||||
const MITC4DirectorFrame& frame = geometry.nodal_frames[0];
|
||||
return {frame.v1, frame.v2, frame.vn};
|
||||
}
|
||||
|
||||
struct NaturalDerivatives {
|
||||
|
||||
Reference in New Issue
Block a user