Skip to content

fix(core): Explicit names for the elastic models #343

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

Closed
sanguinariojoe wants to merge 1 commit into from
Closed

fix(core): Explicit names for the elastic models #343

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
37 changes: 18 additions & 19 deletions source/Line.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ Line::setup(int number_in,
d = props->d;
A = pi / 4. * d * d;
rho = props->w / A;
ElasticMod = props->ElasticMod;
ElasticMod = (elastic_model)props->ElasticMod;
EA = props->EA;
EA_D = props->EA_D;
alphaMBL = props->alphaMBL;
Expand Down Expand Up @@ -242,7 +242,7 @@ Line::setup(int number_in,
// segments) (1 = fully submerged, 0 = out of water)

// viscoelastic things
if (ElasticMod > 1) {
if (ElasticMod != ELASTIC_CONSTANT) {
dl_1.assign(
N,
0.0); // segment stretch attributed to static stiffness portion [m]
Expand Down Expand Up @@ -615,7 +615,7 @@ Line::initialize()
// the segment. This is required here to initalize the state as non-zero,
// which avoids an initial transient where the segment goes from unstretched
// to stretched in one time step.
if (ElasticMod > 1) {
if (ElasticMod != ELASTIC_CONSTANT) {
for (unsigned int i = 0; i < N; i++) {
lstr[i] = unitvector(qs[i], r[i], r[i + 1]);
dl_1[i] = lstr[i] - l[i]; // delta l of the segment
Expand Down Expand Up @@ -790,32 +790,32 @@ Line::setState(const InstanceStateVarView state)

// Error check for number of columns (if VIV and Visco need row.size() = 8,
// if VIV xor Visco need row.size() = 7, if not VIV need row.size() = 6)
if ((state.row(0).size() != 8 && Cl > 0 && ElasticMod > 1) ||
(state.row(0).size() != 7 && ((Cl > 0) ^ (ElasticMod > 1))) ||
(state.row(0).size() != 6 && Cl == 0 && ElasticMod == 1)) {
if ((state.row(0).size() != 8 && Cl > 0 && ElasticMod != ELASTIC_CONSTANT) ||
(state.row(0).size() != 7 && ((Cl > 0) ^ (ElasticMod != ELASTIC_CONSTANT))) ||
(state.row(0).size() != 6 && Cl == 0 && ElasticMod == ELASTIC_CONSTANT)) {
LOGERR << "Invalid state.row size for Line " << number << endl;
throw moordyn::mem_error("Invalid state.row size");
}

// Error check for number of rows (if visco need N rows, if normal need N-1
// rows)
if ((state.rows() != N && ElasticMod > 1) ||
(state.rows() != N - 1 && ElasticMod == 1)) {
if ((state.rows() != N && ElasticMod != ELASTIC_CONSTANT) ||
(state.rows() != N - 1 && ElasticMod == ELASTIC_CONSTANT)) {
LOGERR << "Invalid number of rows in state matrix for Line " << number
<< endl;
throw moordyn::mem_error("Invalid number of rows in state matrix");
}

// If using the viscoelastic model, interate N rows, else iterate N-1 rows.
for (unsigned int i = 0; i < (ElasticMod > 1 ? N : N - 1); i++) {
for (unsigned int i = 0; i < (ElasticMod != ELASTIC_CONSTANT ? N : N - 1); i++) {
// node number is i+1
// segment number is i
if (i < N - 1) { // only assign the internal nodes
r[i + 1] = state.row(i).head<3>();
rd[i + 1] = state.row(i).segment<3>(3);
}

if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
dl_1[i] =
state.row(i)
.tail<1>()[0]; // [0] needed becasue tail<1> returns a one
Expand All @@ -826,7 +826,7 @@ Line::setState(const InstanceStateVarView state)
!(IC_gen)) { // not needed in IC_gen. Initializes as distribution on
// 0-2pi. State is initialized by init function in this
// code, which sets phi to range 0-2pi
if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
phi[i + 1] =
state.row(i)
.tail<2>()[0]; // if both VIV and viscoelastic second to
Expand Down Expand Up @@ -1010,7 +1010,7 @@ Line::getStateDeriv(InstanceStateVarView drdt)
// V[i] = A * l[i]; // volume attributed to segment

// Calculate segment stiffness
if (ElasticMod == 1) {
if (ElasticMod == ELASTIC_CONSTANT) {
// line tension
if (nEApoints > 0)
EA = getNonlinearEA(lstr[i], l[i]);
Expand All @@ -1028,18 +1028,17 @@ Line::getStateDeriv(InstanceStateVarView drdt)
BA = getNonlinearBA(ldstr[i], l[i]);
Td[i] = BA * ldstr[i] / l[i] * qs[i];

} else if (
ElasticMod >
1) { // viscoelastic model from
// https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018
} else {
// viscoelastic model from
// https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018
// note that dl_1[i] is the same as Line%dl_1 in MD-F. This is the
// deltaL of the first static spring k1.

if (ElasticMod == 2) {
if (ElasticMod == ELASTIC_VISCO_CTE) {
// constant dynamic stiffness
EA_2 = EA_D;

} else if (ElasticMod == 3) {
} else if (ElasticMod == ELASTIC_VISCO_MEAN) {
if (dl_1[i] >= 0.0) // spring k1 is in tension
// Mean load dependent dynamic stiffness: from combining
// eqn. 2 and eqn. 10 from original MD viscoelastic paper,
Expand Down Expand Up @@ -1468,7 +1467,7 @@ Line::getStateDeriv(InstanceStateVarView drdt)
// Update state derivative for VIV. i-1 indexing because this is
// only called for internal nodes (i.e. node 1 maps to row 0 in the
// state deriv matrix).
if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
drdt.row(i - 1).tail<2>()[0] =
phi_dot[i]; // second to last element if visco model
else
Expand Down
45 changes: 28 additions & 17 deletions source/Line.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,18 @@ class DECLDIR Line final
~Line();

private:
/// @brief Elasticity models
typedef enum
{
/// Constant EA
ELASTIC_CONSTANT = 1,
/// viscoelastic model with constant dynamic stiffness
ELASTIC_VISCO_CTE = 2,
/// mean load dependent dynamic stiffness
ELASTIC_VISCO_MEAN = 3,
} elastic_model;


/** @brief Get the non-linear stiffness. This is interpolated from a
* curve provided in the input file.
* @param l_stretched The actual length of the segment
Expand Down Expand Up @@ -142,9 +154,8 @@ class DECLDIR Line final
moordyn::real d;
/// line density (kg/m^3)
moordyn::real rho;
/// Elasticity model flag (1: constant EA, 2: viscoelastic model with
/// constant dynamic stiffness, 3: mean load depenedent dynamic stiffness)
unsigned int ElasticMod;
/// Elasticity model flag. See ::elastic_model
elastic_model ElasticMod;
/// line normal/static elasticity modulus * crosssectional area [N]
moordyn::real EA;
/// constant line dynamic stiffness modulus * area for viscoelastic stuff
Expand All @@ -158,8 +169,8 @@ class DECLDIR Line final
moordyn::real vbeta;
/// stiffness of spring 2 in viscoelastic model (dynamic stiffness). This is
/// the spring in series with the parallel spring-dashpot. if ElasticMod =
/// 2, EA_2 = EA_D. If ElasticMod = 3, EA_2 is load dependent dynamic
/// stiffness.
/// ELASTIC_VISCO_CTE, EA_2 = EA_D. If ElasticMod = ELASTIC_VISCO_MEAN, EA_2
/// is load dependent dynamic stiffness.
moordyn::real EA_2;
/// segment stretch attributed to static stiffness portion [m] (deltaL_1)
std::vector<moordyn::real> dl_1;
Expand Down Expand Up @@ -435,38 +446,38 @@ class DECLDIR Line final
// Line::setState but flipped) ------ Error check for number of columns
// (if VIV and Visco need row.size() = 8, if VIV xor Visco need
// row.size() = 7, if not VIV need row.size() = 6)
if ((state.row(0).size() != 8 && Cl > 0 && ElasticMod > 1) ||
(state.row(0).size() != 7 && ((Cl > 0) ^ (ElasticMod > 1))) ||
(state.row(0).size() != 6 && Cl == 0 && ElasticMod == 1)) {
if ((state.row(0).size() != 8 && Cl > 0 && ElasticMod != ELASTIC_CONSTANT) ||
(state.row(0).size() != 7 && ((Cl > 0) ^ (ElasticMod != ELASTIC_CONSTANT))) ||
(state.row(0).size() != 6 && Cl == 0 && ElasticMod == ELASTIC_CONSTANT)) {
LOGERR << "Invalid state.row size for Line " << number << endl;
throw moordyn::mem_error("Invalid state.row size");
}

// Error check for number of rows (if visco need N rows, if normal need
// N-1 rows)
if ((state.rows() != N && ElasticMod > 1) ||
(state.rows() != N - 1 && ElasticMod == 1)) {
if ((state.rows() != N && ElasticMod != ELASTIC_CONSTANT) ||
(state.rows() != N - 1 && ElasticMod == ELASTIC_CONSTANT)) {
LOGERR << "Invalid number of rows in state matrix for Line "
<< number << endl;
throw moordyn::mem_error("Invalid number of rows in state matrix");
}

// If using the viscoelastic model, iterate N rows, else iterate N-1
// rows.
for (unsigned int i = 0; i < (ElasticMod > 1 ? N : N - 1); i++) {
for (unsigned int i = 0; i < (ElasticMod != ELASTIC_CONSTANT ? N : N - 1); i++) {
// node number is i+1
// segment number is i
state.row(i).head<3>() = r[i + 1];
state.row(i).segment<3>(3) = rd[i + 1];

if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
state.row(i).tail<1>()[0] =
dl_1[i]; // [0] needed becasue tail<1> returns a one element
// vector. Viscoelastic state is always the last
// element in the row

if (Cl > 0) {
if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
state.row(i).tail<2>()[0] =
phi[i + 1]; // if both VIV and viscoelastic second to
// last element in the row
Expand Down Expand Up @@ -561,7 +572,7 @@ class DECLDIR Line final
*/
inline void setConstantEA(moordyn::real EA_in)
{
if (ElasticMod > 1) {
if (ElasticMod != ELASTIC_CONSTANT) {
LOGERR << "Cannot set constant EA for viscoelastic model" << endl;
throw moordyn::invalid_value_error(
"Cannot set constant EA for viscoelastic model");
Expand Down Expand Up @@ -1014,7 +1025,7 @@ class DECLDIR Line final
*/
inline const size_t stateN() const
{
if (ElasticMod > 1)
if (ElasticMod != ELASTIC_CONSTANT)
return getN(); // N rows for viscoelastic case
else
return getN() - 1; // N-1 rows for other cases
Expand All @@ -1032,10 +1043,10 @@ class DECLDIR Line final
*/
inline const size_t stateDims() const
{
if (Cl > 0 && ElasticMod > 1)
if (Cl > 0 && ElasticMod != ELASTIC_CONSTANT)
return 8; // 3 for position, 3 for velocity, 1 for VIV phase, 1 for
// viscoelasticity
else if ((Cl > 0) ^ (ElasticMod > 1))
else if ((Cl > 0) ^ (ElasticMod != ELASTIC_CONSTANT))
return 7; // 3 for position, 3 for velocity, 1 for VIV phase or
// viscoelasticity
else
Expand Down
Loading