Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[FluidBio][Bugfix] OSI-related indicators bugfix #11698

Merged
merged 7 commits into from
Oct 20, 2023
Merged

[FluidBio][Bugfix] OSI-related indicators bugfix #11698

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
cde34a8
Some corrections
rubenzorrilla Oct 15, 2023
c8c0eff
Merge branch 'master' into fluidbio/enhance-wss-calculation
rubenzorrilla Oct 16, 2023
7e893be
Merge branch 'master' into fluidbio/enhance-wss-calculation
rubenzorrilla Oct 18, 2023
fb0d54f
Merge branch 'master' into fluidbio/enhance-wss-calculation
rubenzorrilla Oct 19, 2023
cb583bd
Update test
rubenzorrilla Oct 19, 2023
3bb772b
Update applications/FluidDynamicsBiomedicalApplication/custom_utiliti…
rubenzorrilla Oct 19, 2023
671a936
Update wss_statistics_utilities.cpp
rubenzorrilla Oct 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 17 additions & 14 deletions ...ications/FluidDynamicsBiomedicalApplication/custom_utilities/wss_statistics_utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,33 +98,36 @@ void WssStatisticsUtilities::CalculateWSS(

void WssStatisticsUtilities::CalculateOSI(ModelPart &rModelPart)
{
KRATOS_ERROR_IF_NOT(rModelPart.GetProcessInfo().Has(STEP)) << "'STEP' is not present in '" << rModelPart.FullName() << "' ProcessInfo container." << std::endl;
const unsigned int step = rModelPart.GetProcessInfo()[STEP];
block_for_each(rModelPart.Nodes(), [&step](NodeType& rNode){
const double abs_tol = std::numeric_limits<double>::epsilon();
KRATOS_ERROR_IF_NOT(rModelPart.GetProcessInfo().Has(TIME)) << "'TIME' is not present in '" << rModelPart.FullName() << "' ProcessInfo container." << std::endl;
KRATOS_ERROR_IF_NOT(rModelPart.GetProcessInfo().Has(DELTA_TIME)) << "'DELTA_TIME' is not present in '" << rModelPart.FullName() << "' ProcessInfo container." << std::endl;
const double time = rModelPart.GetProcessInfo()[TIME];
const double delta_time = rModelPart.GetProcessInfo()[DELTA_TIME];
block_for_each(rModelPart.Nodes(), [&time, &delta_time, &abs_tol](NodeType& rNode){

// Accumulate the current step contribution to the Oscillatory Shear Index (OSI)
// Note that this requires the tangential WSS to be already computed (see CalculateWSS)
auto& r_temporal_osi = rNode.GetValue(TEMPORAL_OSI);
const auto& r_wss_tang_stress = rNode.GetValue(WSS_TANGENTIAL_STRESS);
r_temporal_osi += (r_wss_tang_stress - r_temporal_osi)/step;
r_temporal_osi += r_wss_tang_stress*delta_time;

// Calculates the sum of the WSS magnitudes for all time steps
double& r_twss = rNode.GetValue(TWSS);
double& r_tawss = rNode.GetValue(TAWSS);
r_twss += (norm_2(r_wss_tang_stress) - r_twss)/step;
r_tawss = norm_2(r_temporal_osi / step);
r_twss += norm_2(r_wss_tang_stress)*delta_time;
r_tawss = r_twss / time;

// Calculate OSI and OSI-dependent magnitudes
double& r_osi = rNode.GetValue(OSI);
r_osi = r_tawss / r_twss > 1.0 ? 0.0 : 0.5 * (1.0 - r_tawss/r_twss);
if (r_twss > 1.0e-12) {
if (std::abs(r_osi - 0.5) < 1.0e-12) {
rNode.GetValue(RRT) = 0.0;
} else {
rNode.GetValue(RRT) = 1.0 / ((1.0 - 2.0*r_osi) * r_twss);
const double temp_osi_norm = norm_2(r_temporal_osi);
r_osi = temp_osi_norm / r_twss > 1.0 ? 0.0 : 0.5 * (1.0 - temp_osi_norm/r_twss);
if (r_tawss > abs_tol) {
rNode.GetValue(ECAP) = r_osi / r_tawss;
if (std::abs(r_osi - 0.5) > abs_tol) {
rNode.GetValue(RRT) = 1.0 / ((1.0 - 2.0*r_osi) * r_tawss);
}
}
rNode.GetValue(ECAP) = r_osi / r_twss;
});
}

}
}
18 changes: 10 additions & 8 deletions ...ions/FluidDynamicsBiomedicalApplication/tests/cpp_tests/test_wss_statistics_utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,8 @@ void TetrahedraModelPartForWSSTests(ModelPart& rModelPart)
// Model part data
rModelPart.AddNodalSolutionStepVariable(REACTION);
rModelPart.GetProcessInfo().SetValue(DOMAIN_SIZE, 3);
rModelPart.GetProcessInfo().SetValue(STEP, 2); // Required as WSS checks step > buffer
rModelPart.GetProcessInfo().SetValue(TIME, 2.0); // Required for the OSI-related indicators
rModelPart.GetProcessInfo().SetValue(DELTA_TIME, 2.0); // Required for the OSI-related indicators

// Geometry creation
auto p_point_1 = Kratos::make_intrusive<Node>(1, 0.0, 0.0, 0.0);
Expand Down Expand Up @@ -109,23 +110,24 @@ KRATOS_TEST_CASE_IN_SUITE(WSSStatisticsUtilitiesOSI, FluidDynamicsBiomedicalAppl
array_1d<double,3> wss_tang_stress;
for (auto& r_node : r_test_skin_model_part.Nodes()) {
temporal_osi = r_node.Coordinates() / 2.0;
wss_tang_stress = r_node.Coordinates() * 2.0;
wss_tang_stress = r_node.Coordinates() * -2.0;
r_node.SetValue(TEMPORAL_OSI, temporal_osi);
r_node.SetValue(WSS_TANGENTIAL_STRESS, wss_tang_stress);
}

// Calculate OSI
WssStatisticsUtilities::CalculateOSI(r_test_skin_model_part);


// Check results
const double tolerance = 1.0e-8;
const auto &r_node = *(r_test_skin_model_part.NodesEnd() - 1);
std::vector<double> expected_temporal_osi = {1.25, 1.25, 1.25};
KRATOS_EXPECT_NEAR(r_node.GetValue(OSI), 0.1875, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(RRT), 0.923760430703, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(TWSS), 1.73205080757, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(ECAP), 0.108253175473, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(TAWSS), 1.08253175473, tolerance);
std::vector<double> expected_temporal_osi = {-3.5,-3.5,-3.5};
KRATOS_EXPECT_NEAR(r_node.GetValue(OSI), 0.0625, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(RRT), 0.329914439537, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(TWSS), 6.92820323028, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(ECAP), 0.0180421959122, tolerance);
KRATOS_EXPECT_NEAR(r_node.GetValue(TAWSS), 3.46410161514, tolerance);
KRATOS_EXPECT_VECTOR_NEAR(r_node.GetValue(TEMPORAL_OSI), expected_temporal_osi, tolerance);
}

Expand Down
Loading