Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
12 changes: 6 additions & 6 deletions doc/readthedocs/boundaryconditions.rst
Original file line number Diff line number Diff line change
Expand Up @@ -213,9 +213,9 @@ The input looks like this:

------------------------------DESIGN SURF ROBIN SPRING DASHPOT CONDITIONS
E <surfset> NUMDOF 3 ONOFF 0 0 0 \
STIFF <X_STIFF> <Y_STIFF> <Z_STIFF> TIMEFUNCTSTIFF 0 0 0 \
VISCO <X_VISCO> <Y_VISCO> <Z_VISCO> TIMEFUNCTVISCO 0 0 0 \
DISPLOFFSET 0.0 0.0 0.0 TIMEFUNCTDISPLOFFSET 0 0 0 FUNCTNONLINSTIFF 0 0 0 \
STIFF <X_STIFF> <Y_STIFF> <Z_STIFF> TIMEFUNCTSTIFF none none none \
VISCO <X_VISCO> <Y_VISCO> <Z_VISCO> TIMEFUNCTVISCO none none none \
DISPLOFFSET 0.0 0.0 0.0 TIMEFUNCTDISPLOFFSET none none none FUNCTNONLINSTIFF none none none \
DIRECTION xyz|refsurfnormal|cursurfnormal COUPLING none

- Commonly the Robin boundary condition couples the nodes at the surface to its original position.
Expand All @@ -236,9 +236,9 @@ The input looks like this:

------------------------------DESIGN SURF ROBIN SPRING DASHPOT CONDITIONS
E <surfset> NUMDOF 3 ONOFF 0 0 0 \
STIFF <X_STIFF> <Y_STIFF> <Z_STIFF> TIMEFUNCTSTIFF 0 0 0 \
VISCO <X_VISCO> <Y_VISCO> <Z_VISCO> TIMEFUNCTVISCO 0 0 0 \
DISPLOFFSET 0.0 0.0 0.0 TIMEFUNCTDISPLOFFSET 0 0 0 FUNCTNONLINSTIFF 0 0 0 \
STIFF <X_STIFF> <Y_STIFF> <Z_STIFF> TIMEFUNCTSTIFF none none none \
VISCO <X_VISCO> <Y_VISCO> <Z_VISCO> TIMEFUNCTVISCO none none none \
DISPLOFFSET 0.0 0.0 0.0 TIMEFUNCTDISPLOFFSET none none none FUNCTNONLINSTIFF none none none \
DIRECTION xyz|refsurfnormal|cursurfnormal COUPLING <couplingID>
//
------------------------------DESIGN SURF ROBIN SPRING DASHPOT COUPLING CONDITIONS
Expand Down
162 changes: 90 additions & 72 deletions src/constraint/4C_constraint_springdashpot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,13 +43,20 @@ namespace
Core::FE::CelltypeSequence<Core::FE::CellType::quad4, Core::FE::CellType::quad8,
Core::FE::CellType::quad9, Core::FE::CellType::tri3, Core::FE::CellType::tri6>;

double scale_by_optional_function(double baseline_value, int num_func, double total_time)
double scale_by_optional_function(
double baseline_value, Core::IO::Noneable<int> funct_id, double total_time)
{
return num_func ? Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(num_func)
.evaluate(total_time) *
baseline_value
: baseline_value;
if (funct_id.has_value())
{
return Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(funct_id.value())
.evaluate(total_time) *
baseline_value;
}
else
{
return baseline_value;
}
}

template <typename T>
Expand Down Expand Up @@ -260,12 +267,15 @@ CONSTRAINTS::SpringDashpot::SpringDashpot(std::shared_ptr<Core::FE::Discretizati

// safety checks of input
const auto* springstiff = &spring_->parameters().get<std::vector<double>>("STIFF");
const auto* numfuncstiff = &spring_->parameters().get<std::vector<int>>("TIMEFUNCTSTIFF");
const auto* numfuncnonlinstiff = &spring_->parameters().get<std::vector<int>>("FUNCTNONLINSTIFF");
const auto* function_id_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTSTIFF");
const auto* function_id_nonlinear_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("FUNCTNONLINSTIFF");

for (unsigned i = 0; i < (*numfuncnonlinstiff).size(); ++i)
for (unsigned i = 0; i < (*function_id_nonlinear_stiffness).size(); ++i)
{
if ((*numfuncnonlinstiff)[i] != 0 and ((*springstiff)[i] != 0 or (*numfuncstiff)[i] != 0))
if ((*function_id_nonlinear_stiffness)[i].has_value() and
((*springstiff)[i] != 0 or *function_id_stiffness)[i].has_value())
FOUR_C_THROW("Must not apply nonlinear stiffness and linear stiffness");
}

Expand Down Expand Up @@ -318,13 +328,16 @@ void CONSTRAINTS::SpringDashpot::evaluate_robin(std::shared_ptr<Core::LinAlg::Sp
// get values and switches from the condition
const auto& onoff = spring_->parameters().get<std::vector<int>>("ONOFF");
const auto& springstiff = spring_->parameters().get<std::vector<double>>("STIFF");
const auto& numfuncstiff = spring_->parameters().get<std::vector<int>>("TIMEFUNCTSTIFF");
const auto& function_id_stiffness =
spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTSTIFF");
const auto& dashpotvisc = spring_->parameters().get<std::vector<double>>("VISCO");
const auto& numfuncvisco = spring_->parameters().get<std::vector<int>>("TIMEFUNCTVISCO");
const auto& function_id_viscosity =
spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTVISCO");
const auto& constant_offset = spring_->parameters().get<std::vector<double>>("DISPLOFFSET");
const auto& numfuncdisploffset =
spring_->parameters().get<std::vector<int>>("TIMEFUNCTDISPLOFFSET");
const auto& numfuncnonlinstiff = spring_->parameters().get<std::vector<int>>("FUNCTNONLINSTIFF");
const auto& function_id_displacement_offset =
spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTDISPLOFFSET");
const auto& function_id_nonlinear_stiffness =
spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("FUNCTNONLINSTIFF");
const auto& direction = spring_->parameters().get<RobinSpringDashpotType>("DIRECTION");

// time-integration factor for stiffness contribution of dashpot, d(v_{n+1})/d(d_{n+1})
Expand Down Expand Up @@ -361,14 +374,14 @@ void CONSTRAINTS::SpringDashpot::evaluate_robin(std::shared_ptr<Core::LinAlg::Sp
if (!onoff_bool[i]) continue;

const double dashpot_viscosity =
scale_by_optional_function(dashpotvisc[i], numfuncvisco[i], total_time);
scale_by_optional_function(dashpotvisc[i], function_id_viscosity[i], total_time);

if (numfuncnonlinstiff[i] > 0)
if (function_id_nonlinear_stiffness[i].has_value())
{
// function is nonlinear
const auto& nonlinear_spring =
Global::Problem::instance()->function_by_id<Core::Utils::FunctionOfSpaceTime>(
numfuncnonlinstiff[i]);
function_id_nonlinear_stiffness[i].value());

spring_dashpot_evaluators[i] =
[=, &nonlinear_spring](double delta_displacement,
Expand All @@ -391,7 +404,7 @@ void CONSTRAINTS::SpringDashpot::evaluate_robin(std::shared_ptr<Core::LinAlg::Sp
{
// function is linear but might depend on time
const double spring_stiffness =
scale_by_optional_function(springstiff[i], numfuncstiff[i], total_time);
scale_by_optional_function(springstiff[i], function_id_stiffness[i], total_time);

spring_dashpot_evaluators[i] =
[=](double delta_displacement,
Expand Down Expand Up @@ -457,7 +470,7 @@ void CONSTRAINTS::SpringDashpot::evaluate_robin(std::shared_ptr<Core::LinAlg::Sp
{
scaled_constant_offset[axis_direction] =
scale_by_optional_function(constant_offset[axis_direction],
numfuncdisploffset[axis_direction], total_time);
function_id_displacement_offset[axis_direction], total_time);
}


Expand Down Expand Up @@ -557,48 +570,47 @@ void CONSTRAINTS::SpringDashpot::evaluate_robin(std::shared_ptr<Core::LinAlg::Sp

// compute stiffness, viscosity, and initial offset from functions
const double dof_stiffness =
(numfuncstiff)[dof] != 0
? springstiff[dof] *
Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(numfuncstiff[dof])
.evaluate(total_time)
function_id_stiffness[dof].has_value()
? springstiff[dof] * Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(
function_id_stiffness[dof].value())
.evaluate(total_time)
: springstiff[dof];
const double dof_viscosity =
numfuncvisco[dof] != 0
? dashpotvisc[dof] *
Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(numfuncvisco[dof])
.evaluate(total_time)
function_id_viscosity[dof].has_value()
? dashpotvisc[dof] * Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(
function_id_viscosity[dof].value())
.evaluate(total_time)
: dashpotvisc[dof];
const double dof_disploffset =
numfuncdisploffset[dof] != 0
? constant_offset[dof] *
Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(numfuncdisploffset[dof])
.evaluate(total_time)
function_id_displacement_offset[dof].has_value()
? constant_offset[dof] * Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(
function_id_displacement_offset[dof].value())
.evaluate(total_time)
: constant_offset[dof];

// displacement related forces and derivatives
double force_disp = 0.0;
double force_disp_deriv = 0.0;
if ((numfuncnonlinstiff)[dof] == 0)
{
force_disp = dof_stiffness * (dof_disp - dof_disploffset);
force_disp_deriv = dof_stiffness;
}
else
if (function_id_nonlinear_stiffness[dof].has_value())
{
std::array<double, 3> displ = {(*disp)[0], (*disp)[1], (*disp)[2]};
force_disp = Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfSpaceTime>(
(numfuncnonlinstiff)[dof] - 1)
function_id_nonlinear_stiffness[dof].value())
.evaluate(displ.data(), total_time, 0);

force_disp_deriv = (Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfSpaceTime>(
(numfuncnonlinstiff)[dof] - 1)
function_id_nonlinear_stiffness[dof].value())
.evaluate_spatial_derivative(displ.data(), total_time, 0))[dof];
}
else
{
force_disp = dof_stiffness * (dof_disp - dof_disploffset);
force_disp_deriv = dof_stiffness;
}

// velocity related forces and derivatives
const double force_vel = dof_viscosity * dof_vel;
Expand Down Expand Up @@ -677,36 +689,39 @@ void CONSTRAINTS::SpringDashpot::evaluate_force(Core::LinAlg::Vector<double>& fi
case RobinSpringDashpotType::cursurfnormal: // spring dashpot acts in curnormal direction

// safety checks
const auto* numfuncstiff = &spring_->parameters().get<std::vector<int>>("TIMEFUNCTSTIFF");
const auto* numfuncvisco = &spring_->parameters().get<std::vector<int>>("TIMEFUNCTVISCO");
const auto* numfuncdisploffset =
&spring_->parameters().get<std::vector<int>>("TIMEFUNCTDISPLOFFSET");
const auto* numfuncnonlinstiff =
&spring_->parameters().get<std::vector<int>>("FUNCTNONLINSTIFF");
for (int dof_numfuncstiff : *numfuncstiff)
const auto* function_id_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTSTIFF");
const auto* function_id_viscosity =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTVISCO");
const auto* function_id_displacement_offset =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>(
"TIMEFUNCTDISPLOFFSET");
const auto* function_id_nonlinear_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("FUNCTNONLINSTIFF");
for (auto dof_function_id_stiffness : *function_id_stiffness)
{
if (dof_numfuncstiff != 0)
if (dof_function_id_stiffness.has_value())
{
FOUR_C_THROW(
"temporal dependence of stiffness not implemented for current surface "
"evaluation");
}
}
for (int dof_numfuncvisco : *numfuncvisco)
for (auto dof_function_id_viscosity : *function_id_viscosity)
{
if (dof_numfuncvisco != 0)
if (dof_function_id_viscosity.has_value())
FOUR_C_THROW(
"temporal dependence of damping not implemented for current surface evaluation");
}
for (int dof_numfuncdisploffset : *numfuncdisploffset)
for (auto dof_function_id_displacement_offset : *function_id_displacement_offset)
{
if (dof_numfuncdisploffset != 0)
if (dof_function_id_displacement_offset.has_value())
FOUR_C_THROW(
"temporal dependence of offset not implemented for current surface evaluation");
}
for (int dof_numfuncnonlinstiff : *numfuncnonlinstiff)
for (auto dof_function_id_nonlinear_stiffness : *function_id_nonlinear_stiffness)
{
if (dof_numfuncnonlinstiff != 0)
if (dof_function_id_nonlinear_stiffness.has_value())
FOUR_C_THROW("Nonlinear spring not implemented for current surface evaluation");
}

Expand Down Expand Up @@ -800,36 +815,39 @@ void CONSTRAINTS::SpringDashpot::evaluate_force_stiff(Core::LinAlg::SparseMatrix
case RobinSpringDashpotType::cursurfnormal: // spring dashpot acts in curnormal direction

// safety checks
const auto* numfuncstiff = &spring_->parameters().get<std::vector<int>>("TIMEFUNCTSTIFF");
const auto* numfuncvisco = &spring_->parameters().get<std::vector<int>>("TIMEFUNCTVISCO");
const auto* numfuncdisploffset =
&spring_->parameters().get<std::vector<int>>("TIMEFUNCTDISPLOFFSET");
const auto* numfuncnonlinstiff =
&spring_->parameters().get<std::vector<int>>("FUNCTNONLINSTIFF");
for (int dof_numfuncstiff : *numfuncstiff)
const auto* function_id_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTSTIFF");
const auto* function_id_viscosity =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("TIMEFUNCTVISCO");
const auto* function_id_displacement_offset =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>(
"TIMEFUNCTDISPLOFFSET");
const auto* function_id_nonlinear_stiffness =
&spring_->parameters().get<std::vector<Core::IO::Noneable<int>>>("FUNCTNONLINSTIFF");
for (auto dof_function_id_stiffness : *function_id_stiffness)
{
if (dof_numfuncstiff != 0)
if (dof_function_id_stiffness.has_value())
{
FOUR_C_THROW(
"temporal dependence of stiffness not implemented for current surface "
"evaluation");
}
}
for (int dof_numfuncvisco : *numfuncvisco)
for (auto dof_function_id_viscosity : *function_id_viscosity)
{
if (dof_numfuncvisco != 0)
if (dof_function_id_viscosity.has_value())
FOUR_C_THROW(
"temporal dependence of damping not implemented for current surface evaluation");
}
for (int dof_numfuncdisploffset : *numfuncdisploffset)
for (auto dof_function_id_displacement_offset : *function_id_displacement_offset)
{
if (dof_numfuncdisploffset != 0)
if (dof_function_id_displacement_offset.has_value())
FOUR_C_THROW(
"temporal dependence of offset not implemented for current surface evaluation");
}
for (int dof_numfuncnonlinstiff : *numfuncnonlinstiff)
for (auto dof_function_id_nonlinear_stiffness : *function_id_nonlinear_stiffness)
{
if (dof_numfuncnonlinstiff != 0)
if (dof_function_id_nonlinear_stiffness.has_value())
FOUR_C_THROW("Nonlinear spring not implemented for current surface evaluation");
}

Expand Down
17 changes: 12 additions & 5 deletions src/fluid/4C_fluid_impedancecondition.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,7 @@ FLD::Utils::FluidImpedanceBc::FluidImpedanceBc(
r1_(impedancecond->parameters().get<double>("R1")),
r2_(impedancecond->parameters().get<double>("R2")),
c_(impedancecond->parameters().get<double>("C")),
functnum_(impedancecond->parameters().get<int>("FUNCT"))
functnum_(impedancecond->parameters().get<Core::IO::Noneable<int>>("FUNCT"))
{
if (myrank_ == 0)
{
Expand Down Expand Up @@ -356,10 +356,17 @@ void FLD::Utils::FluidImpedanceBc::calculate_impedance_tractions_and_update_resi
}
else if (treetype_ == "pressure_by_funct")
{
pressure = Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(functnum_)
.evaluate(time);
Q_np_fac = 0.0;
if (functnum_.has_value())
{
pressure = Global::Problem::instance()
->function_by_id<Core::Utils::FunctionOfTime>(functnum_.value())
.evaluate(time);
Q_np_fac = 0.0;
}
else
{
FOUR_C_THROW("Error with surface impendance condition: Need to define a positive FUNCT.");
}
}
else
{
Expand Down
4 changes: 2 additions & 2 deletions src/fluid/4C_fluid_impedancecondition.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,8 @@ namespace FLD
//! 'material' parameters required for artery tree
const double r1_, r2_, c_;

//! curve number
const int functnum_;
//! function number
const Core::IO::Noneable<int> functnum_;

//! traction vector for impedance bc
std::shared_ptr<Core::LinAlg::Vector<double>> impedancetbc_;
Expand Down
Loading
Loading