feat(shell-surface-pressure-load): step 5 - integrate-nonlinear-follower-pressure
This commit is contained in:
KOKO\Mimi
@@ -0,0 +1,430 @@
|
||||
#include <gtest/gtest.h>
|
||||
#include <hdf5.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <atomic>
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iterator>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "fesa/app/fesa_application.h"
|
||||
#include "fesa/loads/shell_surface_pressure.h"
|
||||
|
||||
namespace {
|
||||
|
||||
class TempDirectory {
|
||||
public:
|
||||
explicit TempDirectory(const std::string& label) {
|
||||
static std::atomic<unsigned long long> sequence{0U};
|
||||
const auto tick =
|
||||
std::chrono::steady_clock::now().time_since_epoch().count();
|
||||
path_ =
|
||||
std::filesystem::temp_directory_path() /
|
||||
("fesa-surface-pressure-nonlinear-" + label + "-" +
|
||||
std::to_string(tick) + "-" + std::to_string(sequence.fetch_add(1U)));
|
||||
std::error_code error;
|
||||
if (!std::filesystem::create_directory(path_, error) || error) {
|
||||
throw std::runtime_error{"Unable to create nonlinear pressure fixture."};
|
||||
}
|
||||
}
|
||||
|
||||
TempDirectory(const TempDirectory&) = delete;
|
||||
TempDirectory& operator=(const TempDirectory&) = delete;
|
||||
|
||||
~TempDirectory() {
|
||||
std::error_code ignored;
|
||||
std::filesystem::remove_all(path_, ignored);
|
||||
}
|
||||
|
||||
const std::filesystem::path& Path() const noexcept { return path_; }
|
||||
|
||||
private:
|
||||
std::filesystem::path path_;
|
||||
};
|
||||
|
||||
class Hdf5File {
|
||||
public:
|
||||
explicit Hdf5File(const std::filesystem::path& path)
|
||||
: handle_{H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT)} {
|
||||
if (handle_ < 0) {
|
||||
throw std::runtime_error{"Unable to open nonlinear pressure results."};
|
||||
}
|
||||
}
|
||||
|
||||
Hdf5File(const Hdf5File&) = delete;
|
||||
Hdf5File& operator=(const Hdf5File&) = delete;
|
||||
|
||||
~Hdf5File() { (void)H5Fclose(handle_); }
|
||||
|
||||
hid_t Get() const noexcept { return handle_; }
|
||||
|
||||
private:
|
||||
hid_t handle_;
|
||||
};
|
||||
|
||||
void WriteText(const std::filesystem::path& path, const std::string& text) {
|
||||
std::ofstream stream{path, std::ios::binary | std::ios::trunc};
|
||||
stream.write(text.data(), static_cast<std::streamsize>(text.size()));
|
||||
if (!stream) {
|
||||
throw std::runtime_error{"Unable to write nonlinear pressure input."};
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<char> ReadBytes(const std::filesystem::path& path) {
|
||||
std::ifstream stream{path, std::ios::binary};
|
||||
return {std::istreambuf_iterator<char>{stream},
|
||||
std::istreambuf_iterator<char>{}};
|
||||
}
|
||||
|
||||
std::string NonlinearPressureDeck(const bool include_pressure,
|
||||
const bool collapsed) {
|
||||
const std::string surface = include_pressure
|
||||
? "*Surface, TYPE=ELEMENT, NAME=Pressure\n"
|
||||
"Shell-1.1, SPOS\n"
|
||||
: "";
|
||||
const std::string dsload =
|
||||
include_pressure ? "*Dsload\nPressure, P, 4.0\n" : "";
|
||||
const std::string node_three =
|
||||
collapsed ? "3, 1., 0., 0.\n" : "3, 1., 1., 0.\n";
|
||||
const std::string constraints =
|
||||
"N1, 1, 6\n"
|
||||
"N2, 1, 6\n"
|
||||
"N3, 2, 6\n"
|
||||
"N3, 1, 1, 1.0\n"
|
||||
"N4, 1, 6\n";
|
||||
return std::string{R"inp(*Part, name=ShellPart
|
||||
*Node
|
||||
1, 0., 0., 0.
|
||||
2, 1., 0., 0.
|
||||
)inp"} + node_three +
|
||||
R"inp(4, 0., 1., 0.
|
||||
*Element, type=S4
|
||||
1, 1, 2, 3, 4
|
||||
*Elset, elset=ShellSet
|
||||
1
|
||||
*Shell Section, elset=ShellSet, material=Steel
|
||||
0.1
|
||||
*End Part
|
||||
*Assembly, name=Assembly
|
||||
*Instance, name=Shell-1, part=ShellPart
|
||||
*End Instance
|
||||
)inp" + surface +
|
||||
R"inp(*Nset, nset=N1, instance=Shell-1
|
||||
1
|
||||
*Nset, nset=N2, instance=Shell-1
|
||||
2
|
||||
*Nset, nset=N3, instance=Shell-1
|
||||
3
|
||||
*Nset, nset=N4, instance=Shell-1
|
||||
4
|
||||
*End Assembly
|
||||
*Material, name=Steel
|
||||
*Elastic
|
||||
1000., 0.25
|
||||
*Boundary
|
||||
)inp" + constraints +
|
||||
R"inp(*Step, name=Step-1, nlgeom=YES
|
||||
*Static
|
||||
0.1, 1., 0.001, 0.1
|
||||
)inp" + dsload +
|
||||
"*End Step\n";
|
||||
}
|
||||
|
||||
std::vector<double> ReadFrameDataset(const Hdf5File& file,
|
||||
const std::uint32_t frame_index,
|
||||
const std::string& suffix) {
|
||||
const std::string path =
|
||||
"/steps/Step-1/frames/" + std::to_string(frame_index) + suffix;
|
||||
const hid_t dataset = H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT);
|
||||
if (dataset < 0) {
|
||||
throw std::runtime_error{"Unable to open nonlinear pressure dataset."};
|
||||
}
|
||||
std::vector<double> values(24U);
|
||||
const herr_t read = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
|
||||
H5P_DEFAULT, values.data());
|
||||
(void)H5Dclose(dataset);
|
||||
if (read < 0 ||
|
||||
!std::all_of(
|
||||
values.begin(), values.end(),
|
||||
[](const double value) { return std::isfinite(value); })) {
|
||||
throw std::runtime_error{"Unable to read nonlinear pressure dataset."};
|
||||
}
|
||||
return values;
|
||||
}
|
||||
|
||||
std::array<double, 6> ReadEquilibrium(const Hdf5File& file,
|
||||
const std::uint32_t frame_index) {
|
||||
const std::string path = "/steps/Step-1/frames/" +
|
||||
std::to_string(frame_index) + "/global/equilibrium";
|
||||
const hid_t dataset = H5Dopen2(file.Get(), path.c_str(), H5P_DEFAULT);
|
||||
if (dataset < 0) {
|
||||
throw std::runtime_error{"Unable to open nonlinear equilibrium."};
|
||||
}
|
||||
std::array<double, 6> values{};
|
||||
const herr_t read = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
|
||||
H5P_DEFAULT, values.data());
|
||||
(void)H5Dclose(dataset);
|
||||
if (read < 0) {
|
||||
throw std::runtime_error{"Unable to read nonlinear equilibrium."};
|
||||
}
|
||||
return values;
|
||||
}
|
||||
|
||||
double SumTranslation(const std::vector<double>& values,
|
||||
const std::size_t component) {
|
||||
double result = 0.0;
|
||||
for (std::size_t node = 0U; node < 4U; ++node) {
|
||||
result += values[node * 6U + component];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
double MaximumAbs(const std::array<double, 6>& values) {
|
||||
double result = 0.0;
|
||||
for (const double value : values) {
|
||||
result = (std::max)(result, std::abs(value));
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
fesa::ShellSurfacePressureInput CurrentPressureInput(
|
||||
const double node_three_x_displacement) {
|
||||
return {
|
||||
{{{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}}},
|
||||
{{{0.0, 0.0, 0.0},
|
||||
{0.0, 0.0, 0.0},
|
||||
{node_three_x_displacement, 0.0, 0.0},
|
||||
{0.0, 0.0, 0.0}}},
|
||||
fesa::ShellSurfaceSide::kSpos,
|
||||
4.0,
|
||||
1.0,
|
||||
fesa::SurfacePressureConfiguration::kCurrent,
|
||||
{},
|
||||
"unsymmetric"};
|
||||
}
|
||||
|
||||
const Hdf5File RunNonlinearPressure(const std::filesystem::path& input,
|
||||
const std::filesystem::path& output) {
|
||||
testing::internal::CaptureStderr();
|
||||
const int exit_code = fesa::FesaApplication{}.Run({input.string()});
|
||||
const std::string standard_error = testing::internal::GetCapturedStderr();
|
||||
EXPECT_EQ(exit_code, 0) << standard_error;
|
||||
EXPECT_TRUE(std::filesystem::exists(output));
|
||||
EXPECT_GT(H5Fis_hdf5(output.string().c_str()), 0);
|
||||
return Hdf5File{output};
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
AssemblesFollowerForceFromCurrentTrialGeometry) {
|
||||
TempDirectory directory{"current-force"};
|
||||
const auto input = directory.Path() / "current-force.inp";
|
||||
const auto output = directory.Path() / "current-force_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
|
||||
|
||||
EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
FormsResidualAndTangentWithApprovedExternalSigns) {
|
||||
TempDirectory directory{"external-sign"};
|
||||
const auto input = directory.Path() / "external-sign.inp";
|
||||
const auto output = directory.Path() / "external-sign_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
|
||||
const auto residual =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_residual");
|
||||
|
||||
EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
|
||||
EXPECT_NEAR(SumTranslation(residual, 2U), 6.0, 1.0e-12);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure, AppliesIncrementFactorExactlyOnce) {
|
||||
TempDirectory directory{"factor-once"};
|
||||
const auto input = directory.Path() / "factor-once.inp";
|
||||
const auto output = directory.Path() / "factor-once_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 0U, "/nodal/generalized_external_force");
|
||||
|
||||
EXPECT_NEAR(SumTranslation(external, 2U), -0.42, 1.0e-12);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure, UsesFullUnsymmetricTangentSolverPath) {
|
||||
TempDirectory directory{"unsymmetric"};
|
||||
const auto input = directory.Path() / "unsymmetric.inp";
|
||||
const auto output = directory.Path() / "unsymmetric_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
|
||||
const auto residual =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_residual");
|
||||
|
||||
EXPECT_NE(SumTranslation(external, 2U), 0.0);
|
||||
EXPECT_NE(SumTranslation(residual, 2U), 0.0);
|
||||
|
||||
const auto central =
|
||||
fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0));
|
||||
const auto forward =
|
||||
fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0 + 1.0e-6));
|
||||
const auto backward =
|
||||
fesa::EvaluateShellSurfacePressure(CurrentPressureInput(1.0 - 1.0e-6));
|
||||
ASSERT_TRUE(central.HasValue());
|
||||
ASSERT_TRUE(forward.HasValue());
|
||||
ASSERT_TRUE(backward.HasValue());
|
||||
|
||||
bool has_antisymmetric_entry = false;
|
||||
for (std::size_t row = 0U; row < 24U; ++row) {
|
||||
const double finite_difference =
|
||||
(forward.Value().force[row] - backward.Value().force[row]) / 2.0e-6;
|
||||
EXPECT_NEAR(finite_difference,
|
||||
central.Value().external_tangent[row * 24U + 12U], 1.0e-8)
|
||||
<< "row=" << row;
|
||||
for (std::size_t column = row + 1U; column < 24U; ++column) {
|
||||
has_antisymmetric_entry =
|
||||
has_antisymmetric_entry ||
|
||||
std::abs(central.Value().external_tangent[row * 24U + column] -
|
||||
central.Value().external_tangent[column * 24U + row]) >
|
||||
1.0e-12;
|
||||
}
|
||||
}
|
||||
EXPECT_TRUE(has_antisymmetric_entry);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
ConvergesOpenPatchAndPreservesAntisymmetricEvidence) {
|
||||
TempDirectory directory{"open-patch"};
|
||||
const auto input = directory.Path() / "open-patch.inp";
|
||||
const auto output = directory.Path() / "open-patch_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
|
||||
|
||||
EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
|
||||
EXPECT_GT(H5Lexists(file.Get(), "/steps/Step-1/accepted_frames", H5P_DEFAULT),
|
||||
0);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
SatisfiesAcceptedFreeResidualAndFactorOneGlobalBalance) {
|
||||
TempDirectory directory{"balance"};
|
||||
const auto input = directory.Path() / "balance.inp";
|
||||
const auto output = directory.Path() / "balance_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const auto equilibrium = ReadEquilibrium(file, 9U);
|
||||
|
||||
EXPECT_LE(MaximumAbs(equilibrium), 1.0e-8);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
IsObjectiveUnderRigidMotionAndTracksCurrentArea) {
|
||||
fesa::ShellSurfacePressureInput input{
|
||||
{{{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}}},
|
||||
{{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, -1.0, 1.0}, {0.0, -1.0, 1.0}}},
|
||||
fesa::ShellSurfaceSide::kSpos,
|
||||
4.0,
|
||||
1.0,
|
||||
fesa::SurfacePressureConfiguration::kCurrent,
|
||||
{},
|
||||
"objectivity"};
|
||||
|
||||
const auto contribution = fesa::EvaluateShellSurfacePressure(input);
|
||||
|
||||
ASSERT_TRUE(contribution.HasValue());
|
||||
double force_x = 0.0;
|
||||
double force_y = 0.0;
|
||||
double force_z = 0.0;
|
||||
for (std::size_t node = 0U; node < 4U; ++node) {
|
||||
force_x += contribution.Value().force[node * 6U];
|
||||
force_y += contribution.Value().force[node * 6U + 1U];
|
||||
force_z += contribution.Value().force[node * 6U + 2U];
|
||||
}
|
||||
EXPECT_NEAR(force_x, 0.0, 1.0e-12);
|
||||
EXPECT_NEAR(force_y, 4.0, 1.0e-12);
|
||||
EXPECT_NEAR(force_z, 0.0, 1.0e-12);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
ReducesSerialAndParallelContributionsDeterministically) {
|
||||
TempDirectory directory{"deterministic"};
|
||||
const auto first_input = directory.Path() / "first.inp";
|
||||
const auto first_output = directory.Path() / "first_results.h5";
|
||||
const auto second_input = directory.Path() / "second.inp";
|
||||
const auto second_output = directory.Path() / "second_results.h5";
|
||||
const std::string deck = NonlinearPressureDeck(true, false);
|
||||
WriteText(first_input, deck);
|
||||
WriteText(second_input, deck);
|
||||
|
||||
const Hdf5File first = RunNonlinearPressure(first_input, first_output);
|
||||
const Hdf5File second = RunNonlinearPressure(second_input, second_output);
|
||||
const auto first_force =
|
||||
ReadFrameDataset(first, 9U, "/nodal/generalized_external_force");
|
||||
const auto second_force =
|
||||
ReadFrameDataset(second, 9U, "/nodal/generalized_external_force");
|
||||
|
||||
EXPECT_EQ(first_force, second_force);
|
||||
EXPECT_NE(SumTranslation(first_force, 2U), 0.0);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
InvalidTrialRollsBackCommittedFrameAndFinalOutput) {
|
||||
TempDirectory directory{"rollback"};
|
||||
const auto input = directory.Path() / "rollback.inp";
|
||||
const auto output = directory.Path() / "rollback_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, true));
|
||||
const std::vector<char> previous{'p', 'r', 'i', 'o', 'r'};
|
||||
{
|
||||
std::ofstream stream{output, std::ios::binary | std::ios::trunc};
|
||||
stream.write(previous.data(),
|
||||
static_cast<std::streamsize>(previous.size()));
|
||||
}
|
||||
|
||||
const int exit_code = fesa::FesaApplication{}.Run({input.string()});
|
||||
|
||||
EXPECT_NE(exit_code, 0);
|
||||
EXPECT_EQ(ReadBytes(output), previous);
|
||||
}
|
||||
|
||||
TEST(NonlinearShellSurfacePressure,
|
||||
WritesSchemaOneExternalInternalResidualReactionAndEquilibrium) {
|
||||
TempDirectory directory{"schema"};
|
||||
const auto input = directory.Path() / "schema.inp";
|
||||
const auto output = directory.Path() / "schema_results.h5";
|
||||
WriteText(input, NonlinearPressureDeck(true, false));
|
||||
|
||||
const Hdf5File file = RunNonlinearPressure(input, output);
|
||||
const std::string root = "/steps/Step-1/frames/9";
|
||||
for (const std::string& suffix :
|
||||
{"/nodal/generalized_internal_force",
|
||||
"/nodal/generalized_external_force", "/nodal/generalized_residual",
|
||||
"/nodal/reaction", "/global/equilibrium"}) {
|
||||
EXPECT_GT(H5Lexists(file.Get(), (root + suffix).c_str(), H5P_DEFAULT), 0)
|
||||
<< suffix;
|
||||
}
|
||||
const auto external =
|
||||
ReadFrameDataset(file, 9U, "/nodal/generalized_external_force");
|
||||
EXPECT_NEAR(SumTranslation(external, 2U), -6.0, 1.0e-12);
|
||||
EXPECT_EQ(H5Lexists(file.Get(), "/surface_pressure", H5P_DEFAULT), 0);
|
||||
}
|
||||
@@ -101,6 +101,38 @@ All, 1, 6
|
||||
dsload + "*End Step\n";
|
||||
}
|
||||
|
||||
std::string NonlinearShellDeck() {
|
||||
return R"inp(*Part, name=ShellPart
|
||||
*Node
|
||||
1, 0., 0., 0.
|
||||
2, 1., 0., 0.
|
||||
3, 1., 1., 0.
|
||||
4, 0., 1., 0.
|
||||
*Element, type=S4
|
||||
1, 1, 2, 3, 4
|
||||
*Elset, elset=ShellSet
|
||||
1
|
||||
*Shell Section, elset=ShellSet, material=Steel
|
||||
0.1
|
||||
*End Part
|
||||
*Assembly, name=Assembly
|
||||
*Instance, name=Shell-1, part=ShellPart
|
||||
*End Instance
|
||||
*Nset, nset=All, instance=Shell-1
|
||||
1, 2, 3, 4
|
||||
*End Assembly
|
||||
*Material, name=Steel
|
||||
*Elastic
|
||||
1000., 0.25
|
||||
*Boundary
|
||||
All, 1, 6
|
||||
*Step, name=Step-1, nlgeom=YES
|
||||
*Static
|
||||
0.1, 1., 0.001, 0.1
|
||||
*End Step
|
||||
)inp";
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
TEST(ShellSurfacePressureApplication,
|
||||
@@ -143,3 +175,21 @@ TEST(ShellSurfacePressureApplication,
|
||||
0);
|
||||
EXPECT_GE(H5Fclose(file), 0);
|
||||
}
|
||||
|
||||
TEST(ShellSurfacePressureApplication,
|
||||
InputWithoutPressureRetainsExistingNonlinearBehavior) {
|
||||
TempDirectory directory{"nonlinear-no-pressure"};
|
||||
const auto input = directory.Path() / "nonlinear-no-pressure.inp";
|
||||
const auto output = directory.Path() / "nonlinear-no-pressure_results.h5";
|
||||
WriteText(input, NonlinearShellDeck());
|
||||
|
||||
EXPECT_EQ(fesa::FesaApplication{}.Run({input.string()}), 0);
|
||||
ASSERT_TRUE(std::filesystem::exists(output));
|
||||
ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0);
|
||||
const hid_t file =
|
||||
H5Fopen(output.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
|
||||
ASSERT_GE(file, 0);
|
||||
EXPECT_GT(H5Lexists(file, "/steps/Step-1/accepted_frames", H5P_DEFAULT), 0);
|
||||
EXPECT_EQ(H5Lexists(file, "/surface_pressure", H5P_DEFAULT), 0);
|
||||
EXPECT_GE(H5Fclose(file), 0);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user